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