# [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.

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).
***************************************************************************
***

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
>
> _______________________________________________
> 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
```