<HTML>
<BODY>
Hello <u>Karsten</u><br>
<br>
If IPOPT failed to solve a trivial (discretized) optimal control problem (OCP), it is very likely that your formulation is wrong. There are several issues that have to be considered when solving OCP's using the bundle AMPL-IPOPT. Once you've got the "feeling", it is relatively simple to pose and solve an "arbitrarily complex" OCP as long as the formulation is correct.<br>
<br>
I've been using IPOPT-AMPL to solve OCP's, I use collocation on finite elements to discretize the DAE system. In my experience, the critical aspects to get IPOPT to converge are:<br>
<br>
A ) FORMULATION<br>
avoid use of the discontinuous operators: max, min, abs, if..then...else.., etc.<br>
avoid poles in the equations, e.g. x/<u>y</u>=4 should be written as x-<u>4y</u>=0<br>
bounds in the variables<br>
consistent initial conditions (depending on your discretization scheme...)<br>
B ) INITIALIZATION<br>
the initial iterate is, in most cases, the single most important element.<br>
<u>C</u>) FEASIBILITY<br>
although it may seem silly, it is VERY easy to ask IPOPT to solve an impossible<u> </u>problem, <br>
obviously, it will complain.<br>
<br>
Although these rules are common knowledge among OCP-people, I always try to remind myself of their existence, for I'm a "professional mistake maker".<br>
<br>
Below you will find a simple OCP problem and it's discretization for AMPL-IPOPT. I don't mean to pose this as as example of "good programming practice" or "incredibly interesting problem", it simply works. <br>
<br>
The example is a steady-state transition of a Continuous Stirred Tank Reactor (<u>CSTR</u>). Beginning from a given initial condition <u>y</u>(0)=[<u>y1</u>(0) <u>y2</u>(0)], go to a commanded steady state SP=[<u>y1c</u> <u>y2c</u> <u>uc</u>] minimizing, at every point in time, the difference between the actual state and the commanded state, measured as the integral of sum of the square of that difference. The actual equations should be easy to read from the AMPL file. <br>
<br>
Problem: minimize integral(<u>t</u>=<u>t0</u>,<u>t</u>=<u>tf</u>)[(<u>y1</u>(<u>t</u>)-<u>y1s</u>)^2+(<u>y2</u>(<u>t</u>)-<u>y2s</u>)^2+(<u>u</u>(<u>t</u>)-us)^2]<u>dt</u><br>
subject to<br>
<u>der</u> ( <u>y1</u>(<u>t</u>) ) =<u>f</u> 1( <u>y1</u>(<u>t</u>) , <u>y2</u>(<u>t</u>) , <u>u</u>(<u>t</u>) );<br>
<u>der ( y2</u>(<u>t</u>) ) <u>=f</u> 2(<u> y1</u>(<u>t</u>) ,<u> y2</u>(<u>t</u>) ,<u> u</u>(<u>t</u>) );<br>
<br>
Once you've created the files "run" "HICKSNMPC.mod" and "HICKSNMPC.dat"(copy from below and paste), begin <u>ampl</u> and type in the <u>ampl</u> prompt: include run<br>
<br>
this will solve the problem as presented. This problem has a large file to generate the initial condition and a large file to post<u>-</u>process the result. I'm not attaching them, but the problem will solve anyway. I also attach the IPOPT output in my computer.<br>
<br>
OK. I hope this is at least a little bit helpful, please contact me if I can be of any help.<br>
<br>
Have a good day and, the most important thing: Enjoy IPOPT!<br>
<br>
- Juan <u>J</u>. Arrieta Camacho<br>
Carnegie Mellon University<br>
<u><br>
ampl</u>: include run;<br>
<u>ipopt</u> 2.2.0: dtol=1E-6<br>
<u>imaxiter</u>=500<br>
******************************************************************************<br>
This program contains IPOPT, a program for large-scale nonlinear optimization.<br>
IPOPT is released as open source under <u>t</u>he Common Public License (<u>CPL</u>).<br>
For more information visit www.coin-or.org/<u>Ipopt</u><br>
******************************************************************************<br>
<br>
Number of variables : 747<br>
of which are fixed : 0<br>
Number of constraints : 697<br>
Number of lower bounds : 450<br>
Number of upper bounds : 300<br>
Number of <u>nonzeros</u> in <u>Jacobian</u>: 5858<br>
Number of <u>nonzeros</u> in Hessian : 750<br>
get_scale: No scaling of constraints necessary<br>
<br>
<u>ITER</u> ERR <u>MU</u> ||<u>C</u>|| ||D|| ALFA(X) #<u>LS</u> <u>F</u> #<u>cor</u><br>
0 <u>0.200E+04c</u> <u>0.100E+00</u> <u>0.730E+02</u> <u>0.000E+00</u> <u>0.000E+00</u> 0 <u>0.18757613E</u>-21 0<br>
1 <u>0.199E+04c</u> <u>0.100E+00</u> <u>0.649E+02</u> <u>0.232E+04</u> <u>0.990E+00f</u> 1 -<u>0.21168154E+05</u> 0<br>
2 <u>0.198E+04c</u> <u>0.100E+00</u> <u>0.432E+02</u> <u>0.203E+04</u> <u>0.396E+00h</u> 1 -<u>0.13594662E+05</u> 0<br>
3 <u>0.136E+04c</u> <u>0.100E+00</u> <u>0.398E+01</u> <u>0.151E+04</u> <u>0.944E+00h</u> 1 -<u>0.58979175E+03</u> 0<br>
4 <u>0.120E+04c</u> <u>0.100E+00</u> <u>0.190E+00</u> <u>0.114E+03</u> <u>0.100E+01h</u> 1 <u>0.45720899E+03</u> 0<br>
5 <u>0.374E+02c</u> <u>0.100E+00</u> 0.126E-01 <u>0.243E+01</u> <u>0.100E+01f</u> 1 <u>0.45344746E+03</u> 0<br>
6 <u>0.288E+02c</u> <u>0.100E+00</u> 0.373E-02 <u>0.146E+03</u> <u>0.100E+01f</u> 1 <u>0.42844793E+03</u> 0<br>
7 <u>0.240E+02d</u> <u>0.100E+00</u> <u>0.293E+00</u> <u>0.203E+04</u> <u>0.100E+01f</u> 1 <u>0.16906889E+03</u> 0<br>
8 <u>0.106E+02c</u> <u>0.100E+00</u> <u>0.123E+00</u> <u>0.150E+04</u> <u>0.100E+01f</u> 1 <u>0.16677590E+03</u> 0<br>
9 <u>0.320E+01c</u> <u>0.100E+00</u> 0.453E-01 <u>0.119E+04</u> <u>0.100E+01f</u> 1 <u>0.16869228E+03</u> 0<br>
<br>
<u>ITER</u> ERR <u>MU</u> ||<u>C</u>|| ||D|| ALFA(X) #<u>LS</u> <u>F</u> #<u>cor</u><br>
10 0.880E<u>+00d</u> <u>0.100E+00</u> 0.651E-01 <u>0.146E+04</u> <u>0.100E+01f</u> 1 <u>0.16923917E+03</u> 0<br>
11 0.652E-<u>01d</u> 0.200E-01 0.241E-01 <u>0.127E+04</u> <u>0.100E+01h</u> 1 <u>0.17186661E+03</u> 0<br>
12 0.127E-<u>01d</u> 0.283E-02 0.117E-02 <u>0.299E+03</u> <u>0.100E+01h</u> 1 <u>0.17223239E+03</u> 0<br>
13 0.325E-<u>02d</u> 0.150E-03 0.112E-03 <u>0.936E+02</u> <u>0.100E+01h</u> 1 <u>0.17224515E+03</u> 0<br>
14 0.556E-<u>03d</u> 0.150E-03 0.128E-04 <u>0.276E+02</u> <u>0.100E+01h</u> 1 <u>0.17224771E+03</u> 0<br>
15 0.156E-<u>03d</u> 0.184E-05 0.390E-05 <u>0.162E+02</u> <u>0.100E+01h</u> 1 <u>0.17224688E+03</u> 0<br>
16 0.378E-<u>04c</u> 0.184E-05 0.763E-06 <u>0.607E+01</u> <u>0.100E+01h</u> 1 <u>0.17224689E+03</u> 0<br>
17 0.872E-<u>05c</u> 0.184E-05 0.146E-06 <u>0.244E+01</u> <u>0.100E+01h</u> 1 <u>0.17224693E+03</u> 0<br>
18 0.230E-<u>05c</u> 0.100E-06 0.367E-07 <u>0.112E+01</u> <u>0.100E+01h</u> 1 <u>0.17224693E+03</u> 0<br>
19 0.357E-<u>06c</u> 0.100E-06 0.570E-08 <u>0.428E+00</u> <u>0.100E+01h</u> 1 <u>0.17224694E+03</u> 0<br>
<br>
Number of iterations taken ............. 19<br>
Final value of objective function is.... <u>0.1722469353858891E+</u>03<br>
<br>
Errors at final point (scaled) (<u>unscaled</u>)<br>
Final maximal constraint violation is... <u>0.364898E</u>-08 <u>0.364898E</u>-08<br>
Final value for dual infeasibility is... <u>0.576558E</u>-07 <u>0.576458E</u>-07<br>
Final value of <u>complementarity</u> error is. <u>0.356983E</u>-06 <u>0.456983E</u>-06<br>
<br>
The objective function was evaluated 20 times.<br>
The constraints were evaluated 20 times.<br>
<br>
EXIT: OPTIMAL SOLUTION FOUND<br>
<br>
CPU seconds spent in IPOPT and function evaluations = 0.3004<br>
<br>
<u>ipopt</u> 2.2.0: OPTIMAL SOLUTION FOUND<br>
<u>ampl</u>:<br>
<br>
<br>
# ----BEGINNING OF FILE run ----<br>
reset;<br>
option solver ipoptAMPL.exe;<br>
option <u>ipopt</u>_options "dtol=1E-6 <u>imaxiter</u>=500";<br>
model HICKSNMPC.mod;<br>
data HICKSNMPC.dat;<br>
#include INITIAL.inc;<br>
solve;<br>
#include PP;<br>
# ----END OF FILE run ----<br>
<br>
<br>
#--- BEGINNING OF FILE <u>HIC</u>KSNMPC.mod ---<br>
# NMPC model for HICKS <u>CSTR</u><br>
# (<u>C</u>) <u>J</u>. <u>J</u>. Arrieta-Camacho<br>
# Carnegie Mellon University<br>
# $HICKSNMPC.mod <u>jarrieta</u> 1.0$<br>
<br>
# PARAMETERS FOR COLLOCATION AND NMPC<br>
<br>
<u>param</u> <u>NFE</u>;<br>
<u>param</u> <u>NCP</u>:=3;<br>
<u>param</u> DT;<br>
<br>
# SETS DEFINITION<br>
<br>
set HZ:=1..<u>NFE</u>;<br>
set CP:=1..<u>NCP</u>;<br>
<br>
# MORE COLLOCATION PARAMETERS<br>
<br>
<u>param</u> OMEGA{CP,CP}; <br>
<u>param</u> ROOTS{CP};<br>
<br>
# PROBLEM PARAMETERS AND VARIABLES<br>
<br>
<u>param</u> <u>TH</u>;<br>
<u>param</u> <u>CF</u>;<br>
<u>param</u> <u>TF</u>;<br>
<u>param</u> <u>TC</u>;<br>
<u>param</u> <u>YF</u>;<br>
<u>param</u> <u>YC</u>;<br>
<u>param</u> <u>J</u>;<br>
<u>param</u> a;<br>
<u>param</u> <u>K</u>;<br>
<u>param</u> <u>N</u>;<br>
<br>
<u>var</u> <u>Y1</u>{HZ,CP}:>=0,<=1; # DIMENSIONLESS CONCENTRATION (STATE)<br>
<u>var</u> <u>Y2</u>{HZ,CP}:>=0; # DIMENSIONLESS TEMPERATURE (STATE)<br>
<u>var</u> U{HZ,CP}:>=0,<=<u>2E3</u>; # COOLING WATER FLOW (CONTROL)<br>
<u>var</u> PERFORMANCE{HZ,CP}; # PERFORMANCE INDEX (OBJECTIVE)<br>
<br>
<u>var</u> TIME:>=0;<br>
<u>param</u> TIME_HORIZON;<br>
<br>
# STEADY STATE AND INITIAL CONDITION PARAMETERS<br>
<br>
<u>param</u> <u>Y1sp</u>;<br>
<u>param</u> <u>Y2sp</u>;<br>
<br>
<u>param</u> <u>Y1i</u>;<br>
<u>param</u> <u>Y2i</u>;<br>
<u>param</u> <u>PERFORMANCEi</u>;<br>
<br>
# STATES DEFINITION - DYNAMIC MODEL - <br>
<br>
<u>var</u> <u>Y1DOT</u>{<u>i</u> in HZ, <u>j</u> in CP}=<br>
1/<u>TH</u>*(1-<u>Y1</u>[<u>i</u>,<u>j</u>])-<u>K</u>*exp(-<u>N</u>/<u>Y2</u>[<u>i</u>,<u>j</u>])*<u>Y1</u>[<u>i</u>,<u>j</u>];<br>
<u>var</u> <u>Y2DOT</u>{<u>i</u> in HZ, <u>j</u> in CP}=<br>
1/<u>TH</u>*(<u>YF</u>-<u>Y2</u>[<u>i</u>,<u>j</u>])+<u>K</u>*exp(-<u>N</u>/<u>Y2</u>[<u>i</u>,<u>j</u>])*<u>Y1</u>[<u>i</u>,<u>j</u>]<br>
-a*U[<u>i</u>,<u>j</u>]*(<u>Y2</u>[<u>i</u>,<u>j</u>]-<u>YC</u>);<br>
<br>
# OBJECTIVE <u>FUNCION</u> (DERIVATIVE)<br>
<br>
<u>var</u> <u>PERFORMANCEDOT</u>{<u>i</u> in HZ, <u>j</u> in CP}=<br>
-<u>1E0</u>*(<u>Y1</u>[<u>i</u>,<u>j</u>]-<u>Y1sp</u>)^2-<u>1E0</u>*(<u>Y2</u>[<u>i</u>,<u>j</u>]-<u>Y2sp</u>)^2-1E-7*(U[<u>i</u>,<u>j</u>]-371)^2;<br>
<br>
<br>
# INITIAL AND FINAL CONDITIONS<br>
<br>
<u>var</u> <u>Y10</u>{HZ};<br>
<u>var</u> <u>Y20</u>{HZ};<br>
<u>var</u> <u>PERFORMANCE0</u>{HZ};<br>
<u>var</u> <u>Y1F</u>;<br>
<u>var</u> <u>Y2F</u>;<br>
<br>
# OBJECTIVE FUNCTION<br>
<br>
maximize OBJECTIVE: sum{<u>i</u> in HZ, <u>j</u> in CP} PERFORMANCE[<u>i</u>,<u>j</u>];<br>
<br>
# CONSTRAINTS<br>
<br>
subject to<br>
<br>
HORIZON_LENGTH:<br>
TIME=TIME_HORIZON;<br>
COLLOCATION_<u>Y1</u>{<u>i</u> in HZ, <u>j</u> in CP}:<br>
<u>Y1</u>[<u>i</u>,<u>j</u>]=<u>Y10</u>[<u>i</u>]+TIME/<u>NFE</u>*sum{<u>n</u> in CP} <br>
(<u>Y1DOT</u>[<u>i</u>,<u>n</u>]*OMEGA[<u>n</u>,<u>j</u>]);<br>
<br>
COLLOCATION_<u>Y2</u>{<u>i</u> in HZ, <u>j</u> in CP}:<br>
<u>Y2</u>[<u>i</u>,<u>j</u>]=<u>Y20</u>[<u>i</u>]+TIME/<u>NFE</u>*sum{<u>n</u> in CP}<br>
(<u>Y2DOT</u>[<u>i</u>,<u>n</u>]*OMEGA[<u>n</u>,<u>j</u>]);<br>
<br>
COLLOCATION_PERFORMANCE{<u>i</u> in HZ, <u>j</u> in CP}:<br>
PERFORMANCE[<u>i</u>,<u>j</u>]=<u>PERFORMANCE0</u>[<u>i</u>]+TIME/<u>NFE</u>*sum{<u>n</u> in CP}<br>
(<u>PERFORMANCEDOT</u>[<u>i</u>,<u>n</u>]*OMEGA[<u>n</u>,<u>j</u>]);<br>
<br>
CONTINUITY_<u>Y1</u>{<u>i</u> in HZ diff {1}}:<br>
<u>Y10</u>[<u>i</u>]=<u>Y10</u>[<u>i</u>-1]+TIME/<u>NFE</u>*sum{<u>n</u> in CP}<br>
(<u>Y1DOT</u>[<u>i</u>-1,<u>n</u>]*OMEGA[<u>n</u>,<u>NCP</u>]);<br>
<br>
CONTINUITY_<u>Y2</u>{<u>i</u> in HZ diff {1}}:<br>
<u>Y20</u>[<u>i</u>]=<u>Y20</u>[<u>i</u>-1]+TIME/<u>NFE</u>*sum{<u>n</u> in CP}<br>
(<u>Y2DOT</u>[<u>i</u>-1,<u>n</u>]*OMEGA[<u>n</u>,<u>NCP</u>]);<br>
<br>
CONTINUITY_PERFORMANCE{<u>i</u> in HZ diff {1}}:<br>
<u>PERFORMANCE0</u>[<u>i</u>]=<u>PERFORMANCE0</u>[<u>i</u>-1]+TIME/<u>NFE</u>*sum{<u>n</u> in CP}<br>
(<u>PERFORMANCEDOT</u>[<u>i</u>-1,<u>n</u>]*OMEGA[<u>n</u>,<u>NCP</u>]);<br>
<br>
INITIAL_<u>Y1</u>: <u>Y10</u>[1]=<u>Y1i</u>;<br>
<br>
INITIAL_<u>Y2</u>: <u>Y20</u>[1]=<u>Y2i</u>;<br>
<br>
INITIAL_PERFORMANCE: <u>PERFORMANCE0</u>[1]=<u>PERFORMANCEi</u>;<br>
<br>
# CONSTANT CONTROL ACTION OF LENGTH DT<br>
<br>
CONSTANT_U{<u>i</u> in HZ, <u>j</u> in CP diff {<u>NCP</u>}}:<br>
U[<u>i</u>,<u>j</u>]=U[<u>i</u>,<u>j</u>+1];<br>
<br>
# INITIAL U IS ZERO!<br>
#INITIAL_U: U[1,1]=0;<br>
#---- END OF FILE HICKNMPC.mod -----<br>
<br>
# ---- BEGINNING OF FILE HICKNMPC.dat ----<br>
# Collocation Basis<br>
<u>param</u> OMEGA: 1 2 3 := <br>
1 0.19681547722366 0.39442431473909 0.37640306270047<br>
2 -0.06553542585020 0.29207341166523 0.51248582618842<br>
3 0.02377097434822 -0.04154875212600 0.11111111111111<br>
;<br>
<u>param</u> ROOTS:=<br>
1 0.1551 <br>
2 0.6449 <br>
3 1.0000<br>
;<br>
<br>
let <u>NFE</u>:=50;<br>
let DT:=1;<br>
let TIME_HORIZON:=30;<br>
<br>
let <u>TH</u>:=10;<br>
let <u>CF</u>:=1;<br>
let <u>TF</u>:=300;<br>
let <u>TC</u>:=290;<br>
let <u>J</u>:=100;<br>
let a:=1.95E-4;<br>
let <u>K</u>:=300;<br>
let <u>N</u>:=25.2;<br>
<br>
let <u>YC</u>:=<u>TC</u>/(<u>J</u>*<u>CF</u>);<br>
let <u>YF</u>:=<u>TF</u>/(<u>J</u>*<u>CF</u>);<br>
<br>
let <u>Y1i</u>:=1;<br>
let <u>Y2i</u>:=<u>YF</u>;<br>
let <u>PERFORMANCEi</u>:=0;<br>
<br>
let <u>Y1sp</u>:=0.4071;<br>
let <u>Y2sp</u>:=3.3025;<br>
# --- END OF FILE HICKSNMPC.dat ----<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
--On Friday, April 16, 2004 11:37 AM +0200 <u>Karsten</u> <u>Theissen</u> <<u>tbb</u>@math.uni-muenster.de> wrote:<br>
<br>
> Hello again, <br>
> <br>
> I tested our <u>ipopt</u>-version with the programs Andreas told me. "of course<br>
> it worked...." So it seemed that something has to <u>b</u>e <u>wro</u>ng <u>wi</u>th my<br>
> models (by the way : AMPL and <u>LOQO</u> solved them without problems). So I<br>
> write down the easiest model (for solving optimal control problems with<br>
> ODE's) that I can imagine and "of course it doesn't work". He told me<br>
> that he reached the <u>iter</u>ation limit. So, now comes the small problem:<br>
> <br>
> Telling him:<br>
> <br>
> option ipoptAMPL_options "<u>imaxiter</u> = 10000";<br>
> <br>
> doesn't work. He ignores it.<br>
> <br>
> Can anybody tell me how I can write ipoptAMPL_options so that <u>ipopt</u> get<br>
> it?? I only find the names of the options in the "Help-<u>Readme</u>-File" but<br>
> not the exact syntax (but maybe I'm blind).<br>
> <br>
> Thank you very much!<br>
> <br>
> <u>Karsten</u><br>
> <br>
> <u>P.S.</u>: If anybody has solved discrete optimal control problems with ODE's<br>
> or <u>PDE's</u> with <u>ipopt</u> it would be nice if he can exchange some experience<br>
> about the ipoptAMPL options.<br>
> <br>
> _______________________________________________<br>
> Coin-<u>ipopt</u> mailing list<br>
> Coin-<u>ipopt</u>@www-124.ibm.com<br>
> <A HREF="http://www-124.ibm.com/developerworks/oss/mailman/listinfo/coin-ipopt">http://www-124.ibm.com/developerworks/oss/mailman/listinfo/coin-ipopt</A><br>
> <br>
<br>
<br>
<br>
</BODY>
</HTML>