TomWiki - User contributions [en]

TomSym Chemical Equilibrium Problem

2011-11-09T12:48:18Z

Mbot:

{{Part Of Manual|title=the TomSym Manual|link=[[TomSym|TomSym Manual]]}}

TomSym implementation of GAMS Example (WALL,SEQ=76)

A Sample Nonlinear system to solve Chemical Equilibrium models.

Wall, T W, Greening, D, and Woolsey, R E D, Solving Complex Chemical Equilibria Using a Geometric-Programming Based Technique. OR 34, 3 (1987).

<source lang="matlab">
toms ba so4 baoh oh hso4 h

r1 = {ba*so4 == 1};
r2 = {baoh/ba/oh == 4.8};
r3 = {hso4/so4/h == .98};
r4 = {h*oh == 1};

b1 = {ba + 1e-7*baoh == so4 + 1e-5*hso4};
b2 = {2*ba + 1e-7*baoh + 1e-2*h == 2*so4 + 1e-5*hso4 + 1e-2*oh};

x0 = {1 == ba; 1 == so4; 1 == baoh
1 == oh; 1 == hso4; 1 == h };

solution = ezsolve(ba,{r1,r2,r3,r4,b1,b2},x0);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.13614 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Problem 1 f_k 1.000004657593424700
sum(|constr|) 0.000000000000846878
f(x_k) + sum(|constr|) 1.000004657594271600
f(x_0) 1.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 4 ConJacEv 4 Iter 3 MinorIter 5
Elapsed time: 0.008000 sec.

</pre>

TomSym Chemical Equilibrium Problem

2011-11-09T12:47:46Z

Mbot:

{{Part Of Manual|title=the TomSym Manual|link=[[TomSym|TomSym Manual]]}}

TomSym implementation of GAMS Example (CHEM,SEQ=21)

The problem is to find the equilibrium composition of a mixture of different atoms.

Bracken, J, and McCormick, G P, Chapter 5. In Selected Applications of Nonlinear Programming. John Wiley and Sons, New York, 1968, pp. 48-49.

c: compounds (H, H2, H2O, N, N2, NH, NO, O, O2, OH)

i: atoms (H hydrogen, N nitrogen, O oxygen)

<source lang="matlab">
% Parameters

% Atoms per compound
a = [1 2 2 0 0 1 0 0 0 1;
0 0 0 1 2 1 1 0 0 0;
0 0 1 0 0 0 1 1 2 1];

% Number of elements in mixture
h=2; n=1; o=1;
mix = [h n o]';

% Gibbs free energy at 3500 k and 750 psi
gibbs = [-10.021; -21.096; -37.986; -9.846; -28.653; ...
-18.918; -28.032; -14.640; -30.594; -26.11];

% Gibbs energy plus pressure
gplusp = gibbs + log(750*0.07031);

% Variables

% Number of mols in mixture
toms 10x1 x

% Component definition
eq1 = {};
for i=1:length(mix)
eq1 = {eq1; sum(a(i,:)*x) == mix(i)};
end

% Energy definition (xb = total number of mols in mixture)
xb = sum(x);

% Total free energy in mixture
energy = sum(x.*(gplusp + log((x)/xb)));

% Bounds on x and xb
cbnd = {0.001 <= x; 0.01 <= xb};

solution = ezsolve(energy,{eq1, cbnd});
</source>

<pre>
Problem type appears to be: con
Time for symbolic processing: 0.074995 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Problem 1 f_k -47.706514830836568000
sum(|constr|) 0.000000000000000666
f(x_k) + sum(|constr|) -47.706514830836568000
f(x_0) -0.209269531292347570

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 26 GradEv 24 Iter 15 MinorIter 37
CPU time: 0.015600 sec. Elapsed time: 0.014000 sec.

</pre>

PROPT Van der Pol Oscillator

2011-11-09T12:13:54Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Restricted second order information for the solution of optimal control problems using control vector parameterization. 2002, Eva Balsa Canto, Julio R. Banga, Antonio A. Alonso Vassilios S. Vassiliadis

Case Study 6.1: van der Pol oscillator

This case study has been studied by several authors, for example Morison, Gritsis, Vassiliadis and Tanartkit and Biegler.

==Problem Description==

The dynamic optimization problem is to minimize:

<math> J = x_3(t_f) </math>

subject to:

<math> \frac{dx_1}{dt} = (1-x_2^2)*x_1-x_2+u </math>

<math> \frac{dx_2}{dt} = x_1 </math>

<math> \frac{dx_3}{dt} = x_1^2+x_2^2+u^2 </math>

<math> -0.3 <= u <= 1.0 </math>

<math> x(t_0) = [0 \ 1 \ 0]' </math>

<math> t_f = 5 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 5, 60);
setPhase(p);

tomStates x1 x2 x3
tomControls u

% Initial guess
x0 = {icollocate({x1 == 0; x2 == 1; x3 == 0})
collocate(u == -0.01)};

% Box constraints
cbox = {-10 <= icollocate(x1) <= 10
-10 <= icollocate(x2) <= 10
-10 <= icollocate(x3) <= 10
-0.3 <= collocate(u) <= 1};

% Boundary constraints
cbnd = initial({x1 == 0; x2 == 1; x3 == 0});

% ODEs and path constraints
ceq = collocate({dot(x1) == (1-x2.^2).*x1-x2+u
dot(x2) == x1; dot(x3) == x1.^2+x2.^2+u.^2});

% Objective
objective = final(x3);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Van Der Pol';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
x3 = subs(collocate(x3),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.11775 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Van Der Pol f_k 2.867259538084708100
sum(|constr|) 0.000000020744545091
f(x_k) + sum(|constr|) 2.867259558829253300
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 26 ConJacEv 26 Iter 23 MinorIter 348
CPU time: 0.218401 sec. Elapsed time: 0.222000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-');
legend('x1','x2','x3');
title('vanDerPol state variables');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('vanDerPol control');
</source>

[[File:vanDerPol_01.png]]

PROPT Two Stage CSTR

2011-11-09T12:13:52Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

Section 6.3.1 Nonlinear two-stage CSTR system

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

==Problem Description==

The system consists of a series of two CSTRs, where there is a transportation delay tau = 0.1 from the first tank to the second. A truncated Taylor series expansion for the time delay.

Find u over t in [0; 2 ] to minimize

<math> J = x_5(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = f_1 </math>

<math> \frac{dx_2}{dt} = f_2 </math>

<math> \frac{dx_3}{dt} = x_1 - x_3 - tau*f_1 - R_2 + 0.25 </math>

<math> \frac{dx_4}{dt} = x_2 - 2*x_4 - u_2*(x_4 + 0.25) - tau*f_2 + R_2 - 0.25 </math>

<math> \frac{dx_5}{dt} = x_1^2+ x_2^2 + x_3^2 + x_4^2 + 0.1*(u_1^2 + u_2^2) </math>

<math> f_1 = 0.5 - x_1 - R_1 </math>

<math> f_2 = -2*(x_2 + 0.25) - u_1*(x_2 + 0.25) + R_1 </math>

<math> R_1 = (x_1 + 0.5)*exp(25*\frac{x_2}{x_2 + 2}) </math>

<math> R_2 = (x_3 + 0.25)*exp(25*\frac{x_4}{x_4 + 2}) </math>

The state variables x1 and x3 are normalized concentration variables in tanks 1 and 2, respectively, and x2 and x4 are normalized temperature variables in tanks 1 and 2, respectively. The variable x5 is introduced to provide the performance index to be minimized.

The initial condition are:

<math> x(0) = [0.15 \ -0.03 \ 0.10 \ 0 \ 0] </math>

<math> -0.5 <= u1 <= 0.5 </math>

<math> -0.5 <= u2 <= 0.5 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 2, 20);
setPhase(p);

tomStates x1 x2 x3 x4 x5
tomControls u1 u2

xi = [0.15;-0.03;0.10;0;0];

% Initial guess
x0 = {icollocate({x1 == xi(1); x2 == xi(2)
x3 == xi(3); x4 == xi(4); x5 == xi(5)})
collocate({u1 == 0; u2 == 0})};

% Box constraints
cbox = collocate({-0.5 <= u1 <= 0.5
-0.5 <= u2 <= 0.5});

% Boundary constraints
cbnd = initial({x1 == xi(1); x2 == xi(2)
x3 == xi(3); x4 == xi(4); x5 == xi(5)});

% ODEs and path constraints
R1 = (x1 + 0.5).*exp(25*x2./(x2 + 2));
R2 = (x3 + 0.25).*exp(25*x4./(x4 + 2));
f1 = 0.5 - x1 - R1;
f2 = -2*(x2 + 0.25) - u1.*(x2 + 0.25) + R1;
tau = 0.1;
ceq = collocate({
dot(x1) == f1; dot(x2) == f2
dot(x3) == x1-x3-tau*f1-R2+0.25
dot(x4) == x2-2*x4-u2.*(x4+0.25)-tau*f2+R2-0.25
dot(x5) == x1.^2+ x2.^2+x3.^2+x4.^2+0.1*(u1.^2+u2.^2)});

% Objective
objective = final(x5);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Two Stage CSTR';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
x3 = subs(collocate(x3),solution);
x4 = subs(collocate(x4),solution);
u1 = subs(collocate(u1),solution);
u2 = subs(collocate(u2),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.45209 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Two Stage CSTR f_k 0.023238023992802660
sum(|constr|) 0.000000172957099895
f(x_k) + sum(|constr|) 0.023238196949902555
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 26 ConJacEv 26 Iter 22 MinorIter 123
CPU time: 0.062400 sec. Elapsed time: 0.068000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Two Stage CSTR state variables');

subplot(2,1,2)
plot(t,u1,'+-',t,u2,'+-');
legend('u1','u2');
title('Two Stage CSTR control');
</source>

[[File:twoStageCSTR_01.png]]

PROPT Two-Phase Schwartz

2011-11-09T12:13:42Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Users Guide for dyn.Opt, Example 4

Schwartz, A. L., Theory and Implementation of Numerical Methods based on Runge-Kutta Integration for Solving Optimal Control Problems. Ph.D. Dissertation, University of California, Berkeley, 1989

==Problem Formulation==

Find u over t in [0; 2.9 ] to minimize

<math> J = 5*(x_1(t_F)^2 + x_2(t_F)^2) </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = u - 0.1*(1+2*x_1^2)*x_2 </math>

<math> x(0) = [1 \ 1] </math>

and path constraints for t<1:

<math> 1-9*(x_1-1)^2-(\frac{x_2-0.4}{0.3})^2 <=0 </math>

<math> -0.8-x_2 <= 0 \ --> \ -0.8 <= x_2</math>

<math> -1 <= u <= 1, (t<1) </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t1
p1 = tomPhase('p1', t1, 0, 1, 25);
toms t2
p2 = tomPhase('p2', t2, 1, 1.9, 25);

setPhase(p1);
tomStates x1p1 x2p1
tomControls up1

setPhase(p2);
tomStates x1p2 x2p2
tomControls up2

setPhase(p1);
% Initial guess
x01 = {icollocate({x1p1 == 1; x2p1 == 1})
collocate(up1==0)};

% Box constraints
cbox1 = {-0.8 <= icollocate(x2p1)
-1 <= collocate(up1) <= 1};

% Boundary constraints
cbnd1 = initial({x1p1 == 1; x2p1 == 1});

% ODEs and path constraints
ceq1 = collocate({
dot(x1p1) == x2p1
dot(x2p1) == up1 - 0.1*(1+2*x1p1.^2).*x2p1
1-9*(x1p1-1).^2-((x2p1-0.4)/0.3).^2 <= 0});

setPhase(p2);
% Initial guess
x02 = {icollocate({x1p2 == 1; x2p2 == 1})
collocate(up2==0)};

% Box constraints
cbox2 = {-50 <= collocate(up2) <= 50};

% ODEs and path constraints
ceq2 = collocate({
dot(x1p2) == x2p2
dot(x2p2) == up2-0.1*(1+2*x1p2.^2).*x2p2});

% Link phase
link = {final(p1,x1p1) == initial(p2,x1p2)
final(p1,x2p1) == initial(p2,x2p2)};

% Objective
objective = 5*(final(p2,x1p2)^2+final(p2,x2p2)^2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Two Phase Schwartz';
constr = {cbox1, cbnd1, ceq1, cbox2, ceq2, link};
solution = ezsolve(objective, constr, {x01, x02}, options);

t = subs(collocate(p1,t1),solution);
t = [t;subs(collocate(p2,t2),solution)];
x1 = subs(collocate(p1,x1p1),solution);
x1 = [x1;subs(collocate(p2,x1p2),solution)];
x2 = subs(collocate(p1,x2p1),solution);
x2 = [x2;subs(collocate(p2,x2p2),solution)];
u = subs(collocate(p1,up1),solution);
u = [u;subs(collocate(p2,up2),solution)];
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.20317 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Two Phase Schwartz f_k 0.000000000000002272
sum(|constr|) 0.000000000003288885
f(x_k) + sum(|constr|) 0.000000000003291157
f(x_0) 10.000000000000014000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 24 ConJacEv 24 Iter 18 MinorIter 361
CPU time: 0.046800 sec. Elapsed time: 0.056000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-');
legend('x1','x2');
title('Two Phase Schwartz state variables');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Two Phase Schwartz control');
</source>

[[File:twoPhaseSchwartz_01.png]]

PROPT Two-Link Robotic Arm

2011-11-09T12:13:39Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

12.4.2 Example 2: Two-link robotic arm

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

==Problem Formulation==

Find u over t in [0; t_F ] to minimize

<math> J = t_F </math>

subject to:

<math> \frac{dx_1}{dt} = \frac{sin(x_3)*(\frac{9}{4}*cos(x_3)*x_1^2+2*x_2^2) + \frac{4}{3}*(u_1-u_2) - \frac{3}{2}*cos(x_3)*u_2 } {\frac{31}{36} + \frac{9}{4}*sin(x_3)^2} </math>

<math> \frac{dx_2}{dt} = -\frac{sin(x_3)*(\frac{7}{2}*x_1^2+\frac{9}{4}*cos(x3)*x_2^2) - \frac{7}{3}*u_2 + \frac{3}{2}*cos(x_3)*(u_1-u_2) }{\frac{31}{36} + \frac{9}{4}*sin(x_3)^2} </math>

<math> \frac{dx_3}{dt} = x_2-x_1 </math>

<math> \frac{dx_4}{dt} = x_1 </math>

The initial condition are:

<math> x(0) = [0 \ 0 \ 0.5 \ 0] </math>

<math> x(t_F) = [0 \ 0 \ 0.5 \ 0.522] </math>

<math> -1 <= u(1:2) <= 1 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t t_f
p = tomPhase('p', t, 0, t_f, 30);
setPhase(p);

tomStates x1 x2 x3 x4
tomControls u1 u2

% Initial guess
x0 = {t_f == 3
icollocate({x1 == 0; x2 == 0
x3 == 0.5; x4 == 0.522})
collocate({u1 == 1-2*t/t_f
u2 == 1-2*t/t_f})};

% Box constraints
cbox = {2.6 <= t_f <= 100
-1 <= collocate(u1) <= 1
-1 <= collocate(u2) <= 1};

% Boundary constraints
cbnd = {initial({x1 == 0; x2 == 0
x3 == 0.5; x4 == 0})
final({x1 == 0; x2 == 0
x3 == 0.5; x4 == 0.522})};

% ODEs and path constraints
ceq = collocate({
dot(x1) == ( sin(x3).*(9/4*cos(x3).*x1.^2+2*x2.^2) ...
+4/3*(u1-u2) - 3/2*cos(x3).*u2 )./ (31/36 + 9/4*sin(x3).^2)
dot(x2) == -( sin(x3).*(7/2*x1.^2 + 9/4*cos(x3).*x2.^2) ...
- 7/3*u2 + 3/2*cos(x3).*(u1-u2) )./ (31/36 + 9/4*sin(x3).^2)
dot(x3) == x2-x1
dot(x4) == x1});

% Objective
objective = t_f;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Two Link Robotic Arm';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
x3 = subs(collocate(x3),solution);
x4 = subs(collocate(x4),solution);
u1 = subs(collocate(u1),solution);
u2 = subs(collocate(u2),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.41192 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Two Link Robotic Arm f_k 2.983364855223867600
sum(|constr|) 0.000000154455553673
f(x_k) + sum(|constr|) 2.983365009679421300
f(x_0) 3.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 20 ConJacEv 20 Iter 16 MinorIter 278
CPU time: 0.062400 sec. Elapsed time: 0.062000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Two Link Robotic Arm state variables');

subplot(2,1,2)
plot(t,u1,'+-',t,u2,'+-');
legend('u1','u2');
title('Two Link Robotic Arm control');
</source>

[[File:twoLinkRoboticArm_01.png]]

PROPT Turbo Generator

2011-11-09T12:13:36Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

OCCAL - A Mixed symbolic-numeric Optimal Control CALcluator

Section 4 Example 1

==Problem Formulation==

Find u over t in [0; t ] to minimize

<math> J = \int_0^{t} (alpha_1*( (x_1-x_1^{s})^2 + (x_4-x_4^{s})^2) + </math>

<math> alpha_2*x_2^{2} + alpha_3*(x_3-x_3^{s})^{2} + </math>

<math> beta_1*(u_1-u_1^{s})^{2} + beta_2*(u_2-u_2^{s})^{2} ) \mathrm{d}t </math>

subject to:

<math> \frac{dx_1}{dt} = x_2*x_4 </math>

<math> \frac{dx_2}{dt} = \frac{1}{M}*(u_1-s_4*x_1*x_4-s_5*x_1*x_3-kappa_d*x_2</math>

<math> \frac{dx_3}{dt} = u_2-A*x_3+c*x_4 </math>

<math> \frac{dx_4}{dt} = -x_1*x_2 </math>

The initial condition are:

<math> x(0) = [x_1^{s} \ x_2^{s} \ x_3^{s} \ x_4^{s}] </math>

<math> x_{1:4}^{s} = [0.60295 \ 0.0 \ 1.87243 \ 0.79778] </math>

<math> alpha = [2.5 \ 1.0 \ 0.1] </math>

<math> beta = [1.0 \ 1.0] </math>

<math> M = 0.04225 </math>

<math> s_{4:5} = [0.0 \ 0.0] </math>

<math> c = 0 </math>

<math> A = 0.17 </math>

<math> u_{1:2}^{s} = [0.80 \ 0.73962] </math>

<math> kappa_d = 0.02535 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 20, 30);
setPhase(p);

tomStates x1 x2 x3 x4
tomControls u1 u2

% Initial guess
x0i = [0.60295;0;1.87243;0.79778];
x0 = {icollocate({x1 == x0i(1); x2 == x0i(2)
x3 == x0i(3); x4 == x0i(4)})
collocate({u1 == 0; u2 == 0})};

% Boundary constraints
cbnd = initial({x1 == x0i(1); x2 == x0i(2)
x3 == x0i(3); x4 == x0i(4)});

% ODEs and path constraints
u1s = 0.80; u2s = 0.73962;
A = 0.17; c = 0;
s4 = 0; s5 = 0;
M = 0.04225; alpha1 = 2.5;
alpha2 = 1.0; alpha3 = 0.1;
beta1 = 1.0; beta2 = 1.0;
kappa_d = 0.02535;

ceq = collocate({dot(x1) == x2.*x4
dot(x2) == 1/M.*(u1-s4*x1.*x4-s5*x1.*x3-kappa_d*x2)
dot(x3) == u2-A*x3+c*x4; dot(x4) == -x1.*x2});

% Objective
objective = integrate(alpha1*( (x1-x0i(1)).^2 + ...
(x4-x0i(4)).^2) + alpha2*x2.^2 + alpha3*(x3-x0i(3)).^2 + ...
beta1*(u1-u1s).^2 + beta2*(u2-u2s).^2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Turbo Generator';
solution = ezsolve(objective, {cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
x3 = subs(collocate(x3),solution);
x4 = subs(collocate(x4),solution);
u1 = subs(collocate(u1),solution);
u2 = subs(collocate(u2),solution);
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.17864 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Turbo Generator f_k 15.019841547670794000
sum(|constr|) 0.000000000046608960
f(x_k) + sum(|constr|) 15.019841547717403000
f(x_0) -57.012069754799789000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 35 ConJacEv 35 Iter 23 MinorIter 119
CPU time: 0.046800 sec. Elapsed time: 0.048000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Turbo Generator state variables');

subplot(2,1,2)
plot(t,u1,'+-',t,u2,'+-');
legend('u1','u2');
title('Turbo Generator control');
</source>

[[File:turboGenerator_01.png]]

PROPT Tubular Reactor

2011-11-09T12:13:32Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Global Optimization of Chemical Processes using Stochastic Algorithms 1996, Julio R. Banga, Warren D. Seider

Case Study V: Global optimization of a bifunctional catalyst blend in a tubular reactor

==Problem Description==

Luus et al and Luus and Bojkov consider the optimization of a tubular reactor in which methylcyclopentane is converted into benzene. The blend of two catalysts, for hydrogenation and isomerization is described by the mass fraction u of the hydrogenation catalyst. The optimal control problem is to find the catalyst blend along the length of the reactor defined using a characteristic reaction time t in the interval 0 <= t <= t_f where t_f = 2000 gh/mol corresponding to the reactor exit such that the benzene concentration at the exit of the reactor is maximized.

Find u(t) to maximize:

<math> J = x_7(t_f) </math>

Subject to:

<math> \frac{dx_1}{dt} = -k_1*x_1 </math>

<math> \frac{dx_2}{dt} = k_1*x_1-(k_2+k_3)*x_2+k_4*x_5 </math>

<math> \frac{dx_3}{dt} = k_2*x_2 </math>

<math> \frac{dx_4}{dt} = -k_6*x_4+k_5*x_5 </math>

<math> \frac{dx_5}{dt} = k_3*x_2+k_6*x_4-(k_4+k_5+k_8+k_9)*x_5+k_7*x_6+k_{10}*x_7 </math>

<math> \frac{dx_6}{dt} = k_8*x_5-k_7*x_6 </math>

<math> \frac{dx_7}{dt} = k_9*x_5-k_{10}*x_7 </math>

where xi, i = 1,..,7 are the mole fractions of the chemical species (i = 1 for methylcyclopentane, i = 7 for benzene), the rate constants are functions of the catalyst blend u(t):

<math> k_i = c(i,1)+c(i,2)*u+c(i,3)*u^2+c(i,4)*u^3 </math>

and the coefficients cij are given by Luus et al. The upper and lower bounds on the mass fraction of the hydrogenation catalyst are.

<math> 0.6 <= u <= 0.9 </math>

The initial vector of mole fractions is:

<math> x(t_0) = [1 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0]' </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.

% Data
c = [0.2918487e-2, -0.8045787e-2, 0.6749947e-2, -0.1416647e-2;...
0.9509977e1, -0.3500994e2, 0.4283329e2, -0.1733333e2; ...
0.2682093e2, -0.9556079e2, 0.1130398e3, -0.4429997e2; ...
0.2087241e3, -0.7198052e3, 0.8277466e3, -0.3166655e3; ...
0.1350005e1, -0.6850027e1, 0.1216671e2, -0.6666689e1; ...
0.1921995e-1, -0.7945320e-1, 0.1105666, -0.5033333e-1;...
0.1323596, -0.4696255, 0.5539323, -0.2166664; ...
0.7339981e1, -0.2527328e2, 0.2993329e2, -0.1199999e2; ...
-0.3950534, 0.1679353e1, -0.1777829e1, 0.4974987; ...
-0.2504665e-4, 0.1005854e-1, -0.1986696e-1, 0.9833470e-2]';

% Independent variable
toms t

% Initial guess
x1opt = 1; x2opt = 0;
x3opt = 0; x4opt = 0;
x5opt = 0; x6opt = 0;
x7opt = 0; uopt = 0.9;

% Loop over number of collocation points
for n=[10 80]
</source>

==Problem setup==

<source lang="matlab">
p = tomPhase('p', t, 0, 2000, n);
setPhase(p);

tomStates x1 x2 x3 x4 x5 x6 x7
tomControls u

% Collocate initial guess
x0 = {icollocate({x1 == x1opt; x2 == x2opt
x3 == x3opt; x4 == x4opt; x5 == x5opt
x6 == x6opt; x7 == x7opt})
collocate(u == uopt)};

% Box constraints
cbox = {icollocate({
0 <= x1 <= 1; 0 <= x2 <= 1
0 <= x3 <= 1; 0 <= x4 <= 1
0 <= x5 <= 1; 0 <= x6 <= 1
0 <= x7 <= 1})
0.6 <= collocate(u) <= 0.9};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 0
x3 == 0; x4 == 0; x5 == 0
x6 == 0; x7 == 0});

% ODEs and path constraints
uvec = [1, u, u.^2, u.^3];
k = (c'*uvec')';
ceq = collocate({
dot(x1) == -k(1)*x1
dot(x2) == k(1)*x1-(k(2)+k(3))*x2+k(4)*x5
dot(x3) == k(2)*x2
dot(x4) == -k(6)*x4+k(5)*x5
dot(x5) == k(3)*x2+k(6)*x4-(k(4)+k(5)+k(8)+k(9))*x5+...
k(7)*x6+k(10).*x7
dot(x6) == k(8)*x5-k(7)*x6
dot(x7) == k(9)*x5-k(10)*x7});

% Objective
objective = -final(x7);
options = struct;
options.name = 'Tubular Reactor';
options.d2c = false;
Prob = sym2prob('con',objective,{cbox, cbnd, ceq}, x0, options);
Prob.xInit = 35; % Use 35 random starting points.
Prob.GO.localSolver = 'snopt';
if n<=10
Result = tomRun('multiMin', Prob, 1);
else
Result = tomRun('snopt', Prob, 1);
end
solution = getSolution(Result);

% Store optimum for use in initial guess
x1opt = subs(x1,solution);
x2opt = subs(x2,solution);
x3opt = subs(x3,solution);
x4opt = subs(x4,solution);
x5opt = subs(x5,solution);
x6opt = subs(x6,solution);
x7opt = subs(x7,solution);
uopt = subs(u,solution);
</source>

<pre>
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Tubular Reactor - Trial 10 f_k -0.010093483541388937
sum(|constr|) 0.000000004530098742
f(x_k) + sum(|constr|) -0.010093479011290195

Solver: multiMin with local solver snopt. EXIT=0. INFORM=0.
Find local optima using multistart local search
Did 35 local tries. Found 1 global, 31 minima. TotFuncEv 1488. TotConstrEv 1418

FuncEv 1488 GradEv 1418 ConstrEv 1418 ConJacEv 23 Iter 907
CPU time: 2.480416 sec. Elapsed time: 2.489000 sec.

</pre>

<pre>
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Tubular Reactor f_k -0.010082374603677158
sum(|constr|) 0.000000000037585022
f(x_k) + sum(|constr|) -0.010082374566092137
f(x_0) -0.010093483541388937

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 46 GradEv 44 ConstrEv 44 ConJacEv 44 Iter 31 MinorIter 531
CPU time: 2.059213 sec. Elapsed time: 2.055000 sec.

</pre>

<source lang="matlab">
end
</source>

==Plot result==

<source lang="matlab">
t = collocate(t);
x7 = collocate(x7opt);
u = collocate(uopt);

subplot(2,1,1)
plot(t,x7,'*-');
legend('x7');
title('Tubular Reactor state variable x7');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Tubular Reactor control');
</source>

[[File:tubularReactor_01.png]]

PROPT Transfer Min Swing

2011-11-09T12:13:29Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Example 7.1: K.L. Teo, K. K. Leong, G.J. Goh

==Problem Formulation==

Find u over t in [0; 1 ] to minimize

<math> J = \int_0^{1} 4.5 (x_3^2 + x_6^2) \mathrm{d}t </math>

subject to:

<math> \frac{dx_1}{dt} = 9*x_4 </math>

<math> \frac{dx_2}{dt} = 9*x_5 </math>

<math> \frac{dx_3}{dt} = 9*x_6 </math>

<math> \frac{dx_4}{dt} = 9*(x_7+17.2656*x_3) </math>

<math> \frac{dx_5}{dt} = 9*x_8 </math>

<math> \frac{dx_6}{dt} = -9*\frac{x_7+27.0756*x_3+2*x_5*x_6}{x_2} </math>

<math> \frac{dx_7}{dt} = 9*u_1 </math>

<math> \frac{dx_8}{dt} = 9*u_2 </math>

<math> x(0) = [0 \ 22 \ 0 \ 0 \ -1 \ 0 \ NaN \ NaN] </math>

<math> x(1) = [10 \ 14 \ 0 \ 2.5 \ 0 \ 0 \ NaN \ NaN] </math>

<math> |x_4| <= 2.5 </math>

<math> |x_5| <= 1.0 </math>

<math> |x_7| <= 2.83374 </math>

<math> -0.80865 <= x_8 <= 0.71265 </math>

<math> |u| <= 10 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t phi1 phi2

% Starting guess
speed = 5;
xopt = 1.2*t;
yopt = 1.6*t;
thetaopt = pi/4;

phi1opt = 1;
phi2opt = 1;
x1opt = 10*t;
x2opt = 22-8*t;
x3opt = 0;
x4opt = 2.5*t;
x5opt = -1+t;
x6opt = 0;
x7opt = 0;
x8opt = 0;
u1opt = 0;
u2opt = 0;
</source>

==Solve the problem, using a successively larger number collocation points==

<source lang="matlab">
for n=[20 40]
</source>

<source lang="matlab">
% Create a new phase and states, using n collocation points
p = tomPhase('p', t, 0, 1, n);
setPhase(p);
tomStates x1 x2 x3 x4 x5 x6 x7 x8
tomControls u1 u2

% Initial guess
x0 = {phi1 == phi1opt; phi2 == phi2opt
icollocate({
x1 == x1opt; x2 == x2opt
x3 == x3opt; x4 == x4opt
x5 == x5opt; x6 == x6opt
x7 == x7opt; x8 == x8opt})
collocate({
u1 == u1opt; u2 == u2opt})};

% Box constraints
cbox = {-10 <= phi1 <= 10
-10 <= phi2 <= 10
-2.5 <= icollocate(x4) <= 2.5
-1 <= icollocate(x5) <= 1
-2.83374 <= icollocate(x7) <= 2.83374
-0.80865 <= icollocate(x8) <= 0.71265
-10 <= collocate(u1) <= 10
-10 <= collocate(u2) <= 10};

% Boundary constraints
cbnd = {initial({x1 == 0
x2 == 22; x3 == 0
x4 == 0; x5 == -1
x6 == 0; x7 == phi1
x8 == phi2
})
final({x1 == 10
x2 == 14; x3 == 0
x4 == 2.5; x5 == 0
x6 == 0})};

% ODEs and path constraints
ceq = collocate({dot(x1) == 9*x4
dot(x2) == 9*x5; dot(x3) == 9*x6
dot(x4) == 9*(x7+17.2656*x3)
dot(x5) == 9*x8
dot(x6) == -9*(x7+27.0756*x3+2*x5.*x6)./x2
dot(x7) == 9*u1; dot(x8) == 9*u2});

% Objective
objective = integrate(4.5*(x3.^2 + x6.^2));
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Transfer Min Swing';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

% Optimal x and u to use as starting guess in the next iteration
x1opt = subs(x1, solution);
x2opt = subs(x2, solution);
x3opt = subs(x3, solution);
x4opt = subs(x4, solution);
x5opt = subs(x5, solution);
x6opt = subs(x6, solution);
x7opt = subs(x7, solution);
x8opt = subs(x8, solution);
u1opt = subs(u1, solution);
u2opt = subs(u2, solution);
phi1opt = subs(phi1, solution);
phi2opt = subs(phi2, solution);
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.25702 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Transfer Min Swing f_k 0.005155677076381541
sum(|constr|) 0.000000000220310318
f(x_k) + sum(|constr|) 0.005155677296691859
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 107 ConJacEv 107 Iter 100 MinorIter 374
CPU time: 0.187201 sec. Elapsed time: 0.186000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.25108 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Transfer Min Swing f_k 0.005157874701172190
sum(|constr|) 0.000000001277942393
f(x_k) + sum(|constr|) 0.005157875979114583
f(x_0) 0.005155729096774491

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 126 ConJacEv 126 Iter 123 MinorIter 594
CPU time: 0.951606 sec. Elapsed time: 0.965000 sec.

</pre>

<source lang="matlab">
end

t = subs(collocate(t),solution);
x1 = collocate(x1opt);
x2 = collocate(x2opt);
x3 = collocate(x3opt);
x7 = collocate(x7opt);
x8 = collocate(x8opt);
</source>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-');
legend('x1','x2','x3');
title('Transfer Min Swing state variables');

subplot(2,1,2)
plot(t,x7,'+-',t,x8,'+-');
legend('x7','x8');
title('Transfer Min Swing control');
</source>

[[File:transferMinSwing_01.png]]

PROPT Time Delay 2 (Approximate)

2011-11-09T12:13:26Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

8.3.2 Example 2

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

Linear time-delay system considered by Palanisamy et al.

==Problem Formulation==

Find u over t in [0; 2 ] to minimize

<math> J = x_2(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = t*x_1 + x_1(t-tau) + u </math>

<math> \frac{dx_2}{dt} = x_1^2 + u^2 </math>

<math> tau = 1 </math>

A Taylor series expansion gives:

<math> \frac{dx_1}{dt} \approx \frac{(t+1)*x_1 + u }{1 + tau} </math>

The initial condition are:

<math> x(0) = [1 \ 0] </math>

<math> -inf <= u <= inf </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 2, 50);
setPhase(p);

tomStates x1 x2
tomControls u

% Initial guess
x0 = {icollocate({x1 == 1; x2 == 0})
collocate(u == 0)};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 0});

% ODEs and path constraints
tau = 1;
ceq = collocate({
dot(x1) == ((t+1).*x1+u)/(1+tau)
dot(x2) == x1.^2+u.^2});

% Objective
objective = final(x2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Time Delay 2 Appr.';
solution = ezsolve(objective, {cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.070509 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Time Delay 2 Appr. f_k 5.340734691399960700
sum(|constr|) 0.000000236393505091
f(x_k) + sum(|constr|) 5.340734927793465500
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 31 ConJacEv 31 Iter 26 MinorIter 138
CPU time: 0.062400 sec. Elapsed time: 0.068000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,u,'+-');
legend('u');
title('Time Delay 2 control');
</source>

[[File:timeDelay2_approximate_01.png]]

PROPT Time Delay 2

2011-11-09T12:13:25Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

8.3.2 Example 2

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

Linear time-delay system considered by Palanisamy et al.

==Problem Formulation==

Find u over t in [0; 2 ] to minimize

<math> J = x_2(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = t*x_1 + x_1(t-tau) + u </math>

<math> \frac{dx_2}{dt} = x_1^2 + u^2 </math>

<math> tau = 1 </math>

The initial condition are:

<math> x(t<=0) = [1 \ 0] </math>

<math> -inf <= u <= inf </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p1 = tomPhase('p1', t, 0, 2, 50);
setPhase(p1);

tomStates x1 x2
tomControls u

% Initial guess
x0 = {icollocate({x1 == 1; x2 == 0})
collocate(u == 0)};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 0});

% Expression for x1(t-tau)
tau = 1;
x1delayed = ifThenElse(t<tau, 1, subs(x1,t,t-tau));

% ODEs and path constraints
ceq = collocate({
dot(x1) == t.*x1 + x1delayed + u
dot(x2) == x1.^2 + u.^2});

% Objective
objective = final(x2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Time Delay 2';
solution = ezsolve(objective, {cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.07238 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Time Delay 2 f_k 4.796108536142885900
sum(|constr|) 0.000000305572513302
f(x_k) + sum(|constr|) 4.796108841715399000
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 32 ConJacEv 32 Iter 27 MinorIter 137
CPU time: 0.062400 sec. Elapsed time: 0.071000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,u,'+-');
legend('u');
title('Time Delay 2 control');
</source>

[[File:timeDelay2_01.png]]

PROPT Time Delay 1 (Approximate)

2011-11-09T12:13:11Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

8.3.1 Example 1

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

Linear time-delay system used for optimal control studies by Chan and Perkins

==Problem Formulation==

Find u over t in [0; 5 ] to minimize

<math> J = x_3(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = -10*x_1-5*x_2-2*x_1(t-tau)-x_2(t-tau)+u </math>

<math> \frac{dx_3}{dt} = 0.5*(10*x_1^2+x_2^2+u^2) </math>

<math> tau = 0.25 </math>

A Taylor series expansion gives:

<math> \frac{dx_2}{dt} \approx (-12*x_1+(2*tau-6)*x_2 + u)/(1-tau) </math>

The initial condition are:

<math> x(0) = [1 \ 1 \ 0] </math>

<math> -inf <= u <= inf </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 5, 50);
setPhase(p);

tomStates x1 x2 x3
tomControls u

% Initial guess
x0 = {icollocate({x1 == 1; x2 == 1; x3 == 0})
collocate(u == 0)};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 1; x3 == 0});

% ODEs and path constraints
tau = 0.25;
ceq = collocate({dot(x1) == x2
dot(x2) == (-12*x1+(2*tau-6)*x2 + u)/(1-tau)
dot(x3) == 0.5*(10*x1.^2+x2.^2+u.^2)});

% Objective
objective = final(x3);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Time Delay 1 Appr.';
solution = ezsolve(objective, {cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.093181 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Time Delay 1 Appr. f_k 2.387051416916647900
sum(|constr|) 0.000000035059748131
f(x_k) + sum(|constr|) 2.387051451976395900
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 67 ConJacEv 67 Iter 48 MinorIter 213
CPU time: 0.249602 sec. Elapsed time: 0.256000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,u,'+-');
legend('u');
title('Time Delay 1 control');
</source>

[[File:timeDelay1_approximate_01.png]]

PROPT Time Delay 1

2011-11-09T12:13:08Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

8.3.1 Example 1

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

Linear time-delay system used for optimal control studies by Chan and Perkins

==Problem Formulation==

Find u over t in [0; 5 ] to minimize

<math> J = x_3(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = -10*x_1-5*x_2-2*x_1(t-tau)-x_2(t-tau)+u </math>

<math> \frac{dx_3}{dt} = 0.5*(10*x_1^2+x_2^2+u^2) </math>

<math> tau = 0.25 </math>

The initial condition are:

<math> x(t<=0) = [1 \ 1 \ 0] </math>

<math> -inf <= u <= inf </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p1 = tomPhase('p1', t, 0, 5, 50);
setPhase(p1);

tomStates x1 x2 x3
tomControls u

% Initial guess
x0 = {icollocate({x1 == 1
x2 == 1; x3 == 0})
collocate(u == 0)};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 1; x3 == 0});

% Expressions for x1(t-tau) and x2(t-tau)
tau = 0.25;
x1delayed = ifThenElse(t<tau, 1, subs(x1,t,t-tau));
x2delayed = ifThenElse(t<tau, 1, subs(x2,t,t-tau));

% ODEs and path constraints
ceq = collocate({dot(x1) == x2
dot(x2) == -10*x1 - 5*x2 - 2*x1delayed - x2delayed + u
dot(x3) == 0.5*(10*x1.^2+x2.^2+u.^2)});

% Objective
objective = final(x3);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Time Delay 1';
solution = ezsolve(objective, {cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.10007 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Time Delay 1 f_k 2.525970860473680400
sum(|constr|) 0.000000011181939804
f(x_k) + sum(|constr|) 2.525970871655620200
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 70 ConJacEv 70 Iter 49 MinorIter 204
CPU time: 0.249602 sec. Elapsed time: 0.251000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,u,'+-');
legend('u');
title('Time Delay 1 control');
</source>

[[File:timeDelay1_01.png]]

PROPT A Simple Terminal Constraint Problem

2011-11-09T12:13:00Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Problem 1: Miser3 manual

==Problem Description==

Find u(t) over t in [0; 1 ] to minimize

<math> J = \int_0^{1} (x^2 + u^2) \mathrm{d}t </math>

subject to:

<math> \frac{dx}{dt} = u </math>

<math> x(0) = 1 </math>

<math> x(1) = 0.75 </math>

<math> x(0.75) = 0.9 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t1
p1 = tomPhase('p1', t1, 0, 0.75, 20);
toms t2
p2 = tomPhase('p2', t2, 0.75, 0.25, 20);

setPhase(p1);
tomStates x1p1
tomControls up1
setPhase(p2);
tomStates x1p2
tomControls up2
setPhase(p1);

% Initial guess
x01 = {icollocate({x1p1 == 1-0.1*t1/0.75})
collocate(up1==0.9*t1/0.75)};

% Box constraints
cbox1 = {-10 <= icollocate(p1,x1p1) <= 10
-10 <= collocate(p1,up1) <= 10};

% Boundary constraints
cbnd1 = initial(x1p1 == 1);

% ODEs and path constraints
ceq1 = collocate(dot(x1p1) == up1);

% Objective
objective1 = integrate(x1p1.^2+up1.^2);

setPhase(p2);
% Initial guess
x02 = {icollocate({x1p2 == 1-0.1*t2})
collocate(up2==0.9+0.1*t2)};

% Box constraints
cbox2 = {-10 <= icollocate(p2,x1p2) <= 10
-10 <= collocate(p2,up2) <= 10};

% Boundary constraints
cbnd2 = {initial(x1p2 == 0.9)
final(x1p2 == 0.75)};

% ODEs and path constraints
ceq2 = collocate(dot(x1p2) == up2);

% Objective
objective2 = integrate(x1p2.^2+up2.^2);

% Objective
objective = objective1 + objective2;

% Link phase
link = {final(p1,x1p1) == initial(p2,x1p2)};
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Terminal Constraint 2';
constr = {cbox1, cbnd1, ceq1, cbox2, cbnd2, ceq2, link};
solution = ezsolve(objective, constr, {x01, x02}, options);

t = subs(collocate(p1,t1),solution);
t = [t;subs(collocate(p2,t2),solution)];
x1 = subs(collocate(p1,x1p1),solution);
x1 = [x1;subs(collocate(p2,x1p2),solution)];
u = subs(collocate(p1,up1),solution);
u = [u;subs(collocate(p2,up2),solution)];
</source>

<pre>
Problem type appears to be: qp
Time for symbolic processing: 0.098918 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: 1: Terminal Constraint 2 f_k 0.920531441471499900
sum(|constr|) 0.000000000474625644
f(x_k) + sum(|constr|) 0.920531441946125570
f(x_0) 0.000000000000000000

Solver: CPLEX. EXIT=0. INFORM=1.
CPLEX Barrier QP solver
Optimal solution found

FuncEv 8 GradEv 8 ConstrEv 8 Iter 8
CPU time: 0.015600 sec. Elapsed time: 0.004000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-');
legend('x1');
title('Terminal Constraint state variable');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Terminal Constraint control');
</source>

[[File:terminalConstraint2_01.png]]

PROPT Temperature Control

2011-11-09T12:12:37Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Optimal Control CY3H2, Lecture notes by Victor M. Becerra, School of Systems Engineering, University of Reading

Heating a room using the least possible energy.

==Problem Description==

Find u over t in [0; 1] to minimize:

<math> J = \frac{1}{2} \int_0^{1} u^2 \mathrm{d}t </math>

subject to:

<math> \frac{dx}{dt} = -2*x + u </math>

<math> x(0) = 0, </math>

<math> x(1) = 10 </math>

<source lang="matlab">
% Copyright (c) 2009-2009 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 1, 20);
setPhase(p);

tomStates x
tomControls u

% Initial guess
x0 = {icollocate(x == 10*t)
collocate(u == 1)};

% Box constraints
cbox = collocate(0 <= u);

% Boundary constraints
cbnd = {initial(x == 0)
final(x == 10)};

% ODEs and path constraints
ceq = collocate(dot(x) == -2*x+u);

% Objective
objective = 0.5*integrate(u^2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Temperature Control';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x = subs(collocate(x),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: qp
Time for symbolic processing: 0.030845 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: 1: Temperature Control f_k 203.731472072683000000
sum(|constr|) 0.000000000059135659
f(x_k) + sum(|constr|) 203.731472072742150000
f(x_0) 0.000000000000000000

Solver: CPLEX. EXIT=0. INFORM=1.
CPLEX Barrier QP solver
Optimal solution found

FuncEv 9 GradEv 9 ConstrEv 9 Iter 9
Elapsed time: 0.003000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1);
subplot(2,1,1)
plot(t,x,'*-');
legend('Temperature');
subplot(2,1,2)
plot(t,u,'*-');
legend('Energy');
</source>

[[File:temperatureControl_01.png]]

PROPT Stirred Tank

2011-11-09T12:12:34Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Users Guide for dyn.Opt, Example 5a, 5b, 5c

Stirred-Tank Chemical Reactor - Kirk, D. E., Optimal control theory: An introduction, Prentice-Hall, 1970.

5a - unconstrained with terminal penalty 5b - unconstrained 5c - control constraint

==Problem Description==

Find u over t in [0; 0.78 ] to minimize

Does not say u^2 in text

5a:

<math> J = \int_0^{0.78} (x_1^2+x_2^2+0.1*u^2)/2 \mathrm{d}t + x_1(t_F)^2 + x_2(t_F)^2 </math>

5b:

<math> J = \int_0^{0.78} (x_1^2+x_2^2+0.1*u^2)/2 \mathrm{d}t </math>

5c:

<math> J = \int_0^{0.78} (x_1^2+x_2^2)/2 \mathrm{d}t </math>

subject to:

<math> a_1 = x_1 + 0.25 </math>

<math> a_2 = x_2 + 0.5 </math>

<math> a_3 = x_1 + 2.0 </math>

<math> a_4 = a_2*exp(25.0*\frac{x_1}{a_3}) </math>

<math> \frac{dx_1}{dt} = -2.0*a_1 + a_4 - a_1*u </math>

<math> \frac{dx_2}{dt} = 0.5 - x_2 - a_4 </math>

<math> x(0) = [0.05 \ 0] </math>

5b, 5c - x(t_F) = [0 0];

5c - <= 1

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t

for i=1:3
</source>

<source lang="matlab">
p = tomPhase('p', t, 0, 0.78, 40);
setPhase(p);

tomStates x1 x2
tomControls u

% Initial guess
x0 = {icollocate({x1 == 0.05; x2 == 0})
collocate(u == 0)};

% Box constraints
cbox = {-1.99 <= icollocate(x1) <= 100
-100 <= icollocate(x2) <= 100
-1000 <= collocate(u) <= 1000};
% x1 cannot be equal to -2, setting to greater
% to avoid singularity in a2*exp(25.0*x1/a3)

% Boundary constraints
cbnd = initial({x1 == 0.05; x2 == 0});

% ODEs and path constraints
a1 = x1 + 0.25; a2 = x2 + 0.5;
a3 = x1 + 2.0; a4 = a2.*exp(25.0*x1./a3);

ceq = collocate({
dot(x1) == -2.0*a1 + a4 - a1.*u
dot(x2) == 0.5 - x2 - a4});
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
if i==1
objective = final(x1)^2+final(x2)^2+...
integrate((x1.^2+x2.^2+0.1*u.^2)/2);
options.name = 'Stirred Tank 5a';
solution1 = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t1 = subs(collocate(t),solution1);
x11 = subs(collocate(x1),solution1);
x21 = subs(collocate(x2),solution1);
u1 = subs(collocate(u),solution1);
elseif i == 2
cbnd = {cbnd; final({x1 == 0; x2 == 0})};
objective = integrate((x1.^2+x2.^2+0.1*u.^2)/2);
options.name = 'Stirred Tank 5b';
solution2 = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
else
cbnd = {cbnd; final({x1 == 0; x2 == 0})};
cbox = {-1.99 <= icollocate(x1) <= 100
-100 <= icollocate(x2) <= 100
-1 <= collocate(u) <= 1};
objective = integrate((x1.^2+x2.^2)/2);
options.name = 'Stirred Tank 5c';
solution3 = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
end
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.20654 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Stirred Tank 5a f_k 0.014213969120012288
sum(|constr|) 0.000000005238862716
f(x_k) + sum(|constr|) 0.014213974358875005
f(x_0) 0.003474999999999964

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 30 ConJacEv 30 Iter 27 MinorIter 113
CPU time: 0.062400 sec. Elapsed time: 0.054000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.19774 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Stirred Tank 5b f_k 0.016702811155814214
sum(|constr|) 0.000000899223600415
f(x_k) + sum(|constr|) 0.016703710379414628
f(x_0) 0.000974999999999998

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 18 ConJacEv 18 Iter 16 MinorIter 121
CPU time: 0.031200 sec. Elapsed time: 0.037000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.18987 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Stirred Tank 5c f_k 0.000989922252663804
sum(|constr|) 0.000000035597654545
f(x_k) + sum(|constr|) 0.000989957850318349
f(x_0) 0.000974999999999998

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 14 ConJacEv 13 Iter 10 MinorIter 139
CPU time: 0.015600 sec. Elapsed time: 0.028000 sec.

</pre>

<source lang="matlab">
end
</source>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t1,x11,'*-',t1,x21,'*-');
legend('x1','x2');
title('Stirred Tank state variables');

subplot(2,1,2)
plot(t1,u1,'+-');
legend('u');
title('Stirred Tank control');
</source>

[[File:stirredTank_01.png]]

PROPT Spring Mass Damper (2 Degree Freedom)

2011-11-09T12:12:31Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

The Direct Approach of General Dynamic Optimal Control: Application on General Software

Tawiwat Veeraklaew, Ph.D. and Settapong Malisuwan, Ph.D. Chulachomklao Royal Military Academy Nakhon-Nayok, Thailand

==Problem Formulation==

Find u over t in [0; 2 ] to minimize

<math> J = u_1+u_2 </math>

subject to:

<math> \frac{dx}{dt} = A*x + B*u </math>

<source lang="matlab">
% Copyright (c) 2007-2009 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 2, 60);
setPhase(p);

tomStates x1 x2 x3 x4
tomControls u1 u2
x = [x1;x2;x3;x4];
u = [u1;u2];

m1 = 1.0; m2 = 1.0; c1 = 1.0; c3 = 1.0;
c2 = 2.0; k1 = 3.0; k2 = 3.0; k3 = 3.0;

B = [0 0; 1/m1 0;
0 0; 0 1/m2];

A = [0 1 0 0;...
1/m1*[-(k1+k2) -(c1+c2) k2 c2];...
0 0 0 1;...
1/m2*[ k2 c2 -(k2+k3) -(c2+c3)]];

x0i = [5; 0; 10; 0];
xfi = [0; 0; 0; 0];

% Box constraints
cbox = {0 <= collocate(u) <= 9};

% Boundary constraints
cbnd = {initial(x == x0i)
final(x == xfi)};

% ODEs and path constraints
ceq = collocate(dot(x) == A*x+B*u);

% Objective
objective = integrate(u1+u2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Spring Mass Damper';
solution = ezsolve(objective, {cbox, cbnd, ceq}, [], options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
x3 = subs(collocate(x3),solution);
x4 = subs(collocate(x4),solution);
u1 = subs(collocate(u1),solution);
u2 = subs(collocate(u2),solution);
</source>

<pre>
Problem type appears to be: lp
Time for symbolic processing: 0.068338 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Spring Mass Damper f_k 16.485256203083797000
sum(|constr|) 0.000000005591243026
f(x_k) + sum(|constr|) 16.485256208675040000
f(x_0) 0.000000000000000000

Solver: CPLEX. EXIT=0. INFORM=1.
CPLEX Dual Simplex LP solver
Optimal solution found

FuncEv 243 Iter 243
CPU time: 0.078000 sec. Elapsed time: 0.041000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Spring Mass Damper state variables');

subplot(2,1,2)
plot(t,u1,'+-',t,u2,'+-');
legend('u1','u2');
title('Spring Mass Damper control');
</source>

[[File:springSystem_01.png]]

PROPT Singular Control 6

2011-11-09T12:12:29Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Viscocity Solutions of Hamilton-Jacobi Equations and Optimal Control Problems. Alberto Bressan, S.I.S.S.A, Trieste, Italy.

A singular control example.

==Problem Description==

Find u over t in [0; 10 ] to maximize:

<math> J = x_3(t_f) </math>

subject to:

<math> \frac{dx_1}{dt} = u </math>

<math> \frac{dx_2}{dt} = -x_1 </math>

<math> \frac{dx_2}{dt} = x_2-x_1^2 </math>

<math> x(t_0) = [0 \ 0 \ 0] </math>

<math> |u| <= 1 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
t_F = 10;
p = tomPhase('p', t, 0, t_F, 80);
setPhase(p);

tomStates x1 x2 x3
tomControls u

x = [x1; x2; x3];

% Initial guess
x0 = {icollocate({x1 == 0, x2 == 0, x3 == 0})
collocate(u==0)};

% Box constraints
cbox = {-1 <= collocate(u) <= 1};

% Boundary constraints
cbnd = initial(x == [0;0;0]);

% ODEs and path constraints
ceq = collocate({dot(x1) == u; dot(x2) == -x(1)
dot(x3) == x(2)-x(1).^2});

% Objective
objective = -final(x(3));
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Singular Control 6';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

% Extract optimal states and controls from solution
t = collocate(subs(t,solution));
u = collocate(subs(u,solution));
x1 = collocate(subs(x1,solution));
x2 = collocate(subs(x2,solution));
x3 = collocate(subs(x3,solution));
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.076014 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Singular Control 6 f_k -55.555568442322517000
sum(|constr|) 0.000000011402048583
f(x_k) + sum(|constr|) -55.555568430920466000
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 58 ConJacEv 58 Iter 49 MinorIter 678
CPU time: 0.733205 sec. Elapsed time: 0.738000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(3,1,1)
plot(t,x1,'*-',t,x2/10,'*-',t,x3/50,'*-');
legend('x1','x2/10','x3/50');
title('Singular Control 6 state variables');

subplot(3,1,2)
plot(t,u,'+-');
legend('u');
title('Singular Control 6 control');

subplot(3,1,3)
plot(t,-1*(t<t_F/3)+1/2*(t>=t_F/3),'*-');
legend('Known u');
title('Singular Control 6 known solution');
</source>

[[File:singularControl6_01.png]]

PROPT Singular Control 5

2011-11-09T12:12:27Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

10.3 Yeo's singular control problem

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

==Problem Formulation==

Find u over t in [0; 1 ] to minimize:

<math> J = x_4(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = -x_3*u + 16*x_5 - 8 </math>

<math> \frac{dx_3}{dt} = u </math>

<math> \frac{dx_4}{dt} = x_1^2 + x_2^2 + 0.0005*(x_2+16*x_5-8-0.1*x_3*u^2)^2 </math>

<math> \frac{dx_5}{dt} = 1 </math>

The initial condition are:

<math> x(0) = [0 \ -1 \ -sqrt(5) \ 0 \ 0] </math>

<math> -4 <= u <= 10 </math>

The state x4 is implemented as a cost directly. x4 in the implementation is x5. u has a low limit of 9 in the code.

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 1, 80);
setPhase(p)

tomStates x1 x2 x3 x4
tomControls u

% Initial guess
x0 = {icollocate({x1 == 0; x2 == -1
x3 == -sqrt(5); x4 == 0})
collocate(u == 3)};

% Box constraints
cbox = {0 <= collocate(u) <= 10};

% Boundary constraints
cbnd = initial({x1 == 0; x2 == -1
x3 == -sqrt(5); x4 == 0});

% ODEs and path constraints
ceq = collocate({dot(x1) == x2
dot(x2) == -x3.*u + 16*x4 - 8
dot(x3) == u; dot(x4) == 1});

% Objective
objective = integrate(x1.^2 + x2.^2 + ...
0.0005*(x2+16*x4-8-0.1*x3.*u.^2).^2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Singular Control 5';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: con
Time for symbolic processing: 0.16034 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Singular Control 5 f_k 0.119318613152343460
sum(|constr|) 0.000000483936371958
f(x_k) + sum(|constr|) 0.119319097088715420
f(x_0) 1.024412849382206200

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 374 GradEv 372 ConstrEv 372 ConJacEv 372 Iter 345 MinorIter 938
CPU time: 4.180827 sec. Elapsed time: 4.185000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,u,'+-');
legend('u');
title('Singular Control 5 control');
</source>

[[File:singularControl5_01.png]]

PROPT Singular Control 4

2011-11-09T12:12:23Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

10.2.3 Example 4

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

==Problem Formulation==

Find u over t in [0; 5 ] to minimize

<math> J = x_4(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = x_3 </math>

<math> \frac{dx_3}{dt} = u </math>

<math> \frac{dx_4}{dt} = x_1^2 </math>

The initial condition are:

<math> x(0) = [1 \ 0 \ 0 \ 0] </math>

<math> -1 <= u <= 1 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 5, 100);
setPhase(p)

tomStates x1 x2 x3 x4
tomControls u

% Initial guess
x0 = {icollocate({x1 == 1; x2 == 0
x3 == 0; x4 == 0})
collocate(u == 0)};

% Box constraints
cbox = {-1 <= collocate(u) <= 1};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 0
x3 == 0; x4 == 0});

% ODEs and path constraints
ceq = collocate({dot(x1) == x2; dot(x2) == x3
dot(x3) == u; dot(x4) == x1.^2});

% Objective
objective = final(x4);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Singular Control 4';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.089315 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Singular Control 4 f_k 1.252389645383044100
sum(|constr|) 0.000000063932037643
f(x_k) + sum(|constr|) 1.252389709315081800
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 92 ConJacEv 92 Iter 89 MinorIter 652
CPU time: 9.313260 sec. Elapsed time: 2.553000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,u,'+-');
legend('u');
title('Singular Control 4 control');
</source>

[[File:singularControl4_01.png]]

PROPT Singular Control 3

2011-11-09T12:12:21Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

10.2.3 Example 3

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

==Problem Formulation==

Find u over t in [0; 5 ] to minimize

<math> J = x_3(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = u </math>

<math> \frac{dx_3}{dt} = x_1^2 + x_2^2 </math>

The initial condition are:

<math> x(0) = [0 \ 1 \ 0] </math>

<math> -1 <= u <= 1 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 5, 60);
setPhase(p);

tomStates x1 x2 x3
tomControls u

% Initial guess
x0 = {icollocate({x1 == 0; x2 == 1; x3 == 0})
collocate(u == 0)};

% Box constraints
cbox = {-1 <= collocate(u) <= 1};

% Boundary constraints
cbnd = initial({x1 == 0; x2 == 1; x3 == 0});

% ODEs and path constraints
ceq = collocate({dot(x1) == x2
dot(x2) == u; dot(x3) == x1.^2 + x2.^2});

% Objective
objective = final(x3);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Singular Control 3';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.08177 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Singular Control 3 f_k 0.753994561590099140
sum(|constr|) 0.000000015978036398
f(x_k) + sum(|constr|) 0.753994577568135590
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 43 ConJacEv 43 Iter 34 MinorIter 366
CPU time: 0.265202 sec. Elapsed time: 0.264000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,u,'+-');
legend('u');
title('Singular Control 3 control');
</source>

[[File:singularControl3_01.png]]

PROPT Singular Control 2

2011-11-09T12:12:18Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

10.2.2 Example 2

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

==Problem Formulation==

Find u over t in [0; 5 ] to minimize

<math> J = x_3(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = u </math>

<math> \frac{dx_3}{dt} = x_1^2 </math>

The initial condition are:

<math> x(0) = [0 \ 1 \ 0] </math>

<math> -1 <= u <= 1 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 5, 50);
setPhase(p);

tomStates x1 x2 x3
tomControls u

% Initial guess
x0 = {icollocate({x1 == 0; x2 == 1; x3 == 0})
collocate(u == 0)};

% Box constraints
cbox = {-1 <= collocate(u) <= 1};

% Boundary constraints
cbnd = initial({x1 == 0; x2 == 1; x3 == 0});

% ODEs and path constraints
ceq = collocate({dot(x1) == x2
dot(x2) == u; dot(x3) == x1.^2});

% Objective
objective = final(x3);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Singular Control 2';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.07418 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Singular Control 2 f_k 0.268336059478540940
sum(|constr|) 0.000000004483166969
f(x_k) + sum(|constr|) 0.268336063961707890
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 107 ConJacEv 107 Iter 100 MinorIter 353
CPU time: 0.343202 sec. Elapsed time: 0.333000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,u,'+-');
legend('u');
title('Singular Control 2 control');
</source>

[[File:singularControl2_01.png]]

PROPT Singular Control 1

2011-11-09T12:12:14Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

10.2.1 Example 1

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

==Problem Formulation==

Find u over t in [0; 2 ] to minimize

<math> J = x_2(t_F) </math>

subject to:

<math> \frac{dx_1}{dt} = u </math>

<math> \frac{dx_2}{dt} = 0.5*x_1^2 </math>

The initial condition are:

<math> x(0) = [1 \ 0] </math>

<math> -1 <= u <= 1 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 2, 50);
setPhase(p);

tomStates x1 x2
tomControls u

% Initial guess
x0 = {icollocate({x1 == 1; x2 == 0})
collocate(u == 0)};

% Box constraints
cbox = {-1 <= collocate(u) <= 1};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 0});

% ODEs and path constraints
ceq = collocate({dot(x1) == u
dot(x2) == 0.5*x1.^2});

% Objective
objective = final(x2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Singular Control 1';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.064118 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Singular Control 1 f_k 0.166665695130345510
sum(|constr|) 0.000000330654862346
f(x_k) + sum(|constr|) 0.166666025785207870
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 41 ConJacEv 41 Iter 39 MinorIter 164
CPU time: 0.078001 sec. Elapsed time: 0.079000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,u,'+-');
legend('u');
title('Singular Control 1 control');
</source>

[[File:singularControl1_01.png]]

PROPT Singular CSTR

2011-11-09T12:12:13Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

10.4 Nonlinear two-stage CSTR problem

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

==Problem Formulation==

Find u over t in [0; t_F ] to minimize:

<math> J = x(t_F)'*x(t_F) + t_F </math>

(the state variables are moved to bounds)

subject to:

<math> \frac{dx_1}{dt} = -3*x_1+g_1 </math>

<math> \frac{dx_2}{dt} = -11.1558*x_2+g_1-8.1558*(x_2+0.1592)*u_1 </math>

<math> \frac{dx_3}{dt} = 1.5*(0.5*x_1-x_3)+g_2 </math>

<math> \frac{dx_4}{dt} = 0.75*x_2-4.9385*x_4+g_2-3.4385*(x_4+0.122)*u_2 </math>

<math> g_1 = 1.5e7*(0.5251-x_1)*exp(-\frac{10}{x_2+0.6932})- </math>

<math> 1.5e10*(0.4748+x_1)*exp(-\frac{15}{x_2+0.6932}) - 1.4280 </math>

<math> g_2 = 1.5e7*(0.4236-x_2)*exp(-\frac{10}{x_4+0.6560})- </math>

<math> 1.5e10*(0.5764+x_3)*exp(-\frac{15}{x4+0.6560}) - 0.5086 </math>

The initial condition are:

<math> x(0) = [0.1962 \ -0.0372 \ 0.0946 \ 0] </math>

<math> -1 <= u(1:2) <= 1 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t t_f
p = tomPhase('p', t, 0, t_f, 30);
setPhase(p)

tomStates x1 x2 x3 x4
tomControls u1 u2

% Initial guess
x0 = {t_f == 0.3
icollocate({x1 == 0.1962; x2 == -0.0372
x3 == 0.0946; x4 == 0})
collocate({u1 == 0; u2 == 0})};

% Box constraints
cbox = {0.1 <= t_f <= 100
-1 <= collocate(u1) <= 1
-1 <= collocate(u2) <= 1};

% Boundary constraints
cbnd = {initial({x1 == 0.1962; x2 == -0.0372
x3 == 0.0946; x4 == 0})
final({x1 == 0; x2 == 0
x3 == 0; x4 == 0})};

% ODEs and path constraints
g1 = 1.5e7*(0.5251-x1).*exp(-10./(x2+0.6932)) ...
- 1.5e10*(0.4748+x1).*exp(-15./(x2+0.6932)) - 1.4280;
g2 = 1.5e7*(0.4236-x2).*exp(-10./(x4+0.6560)) ...
- 1.5e10*(0.5764+x3).*exp(-15./(x4+0.6560)) - 0.5086;

ceq = collocate({
dot(x1) == -3*x1+g1
dot(x2) == -11.1558*x2+g1-8.1558*(x2+0.1592).*u1
dot(x3) == 1.5*(0.5*x1-x3)+g2
dot(x4) == 0.75*x2-4.9385*x4+g2-3.4385*(x4+0.122).*u2});

% Objective
objective = t_f;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Singular CSTR';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
x3 = subs(collocate(x3),solution);
x4 = subs(collocate(x4),solution);
u1 = subs(collocate(u1),solution);
u2 = subs(collocate(u2),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.51719 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Singular CSTR f_k 0.324402684069356740
sum(|constr|) 0.000000010809237995
f(x_k) + sum(|constr|) 0.324402694878594740
f(x_0) 0.299999999999999990

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 75 ConJacEv 75 Iter 42 MinorIter 427
CPU time: 0.187201 sec. Elapsed time: 0.192000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Singular CSTR state variables');

subplot(2,1,2)
plot(t,u1,'+-',t,u2,'+-');
legend('u1','u2');
title('Singular CSTR control');
</source>

[[File:singularCSTR_01.png]]

PROPT Singular Arc Problem

2011-11-09T12:12:10Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Problem 3: Miser3 manual

==Problem Formulation==

Find u(t) over t in [0; t_f ] to minimize

<math> J = t_f </math>

subject to:

<math> \frac{dx_1}{dt} = u </math>

<math> \frac{dx_2}{dt} = cos(x_1) </math>

<math> \frac{dx_3}{dt} = sin(x_1) </math>

<math> x_2(t_f) = x_3(t_f) = 0 </math>

<math> |u| <= 2 </math>

<math> x(0) = [\frac{pi}{2} \ 4 \ 0] </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
toms t_f
p = tomPhase('p', t, 0, t_f, 60);
setPhase(p);

tomStates x1 x2 x3
tomControls u

% Initial guess
x0 = {t_f == 20
icollocate({
x1 == pi/2+pi/2*t/t_f
x2 == 4-4*t/t_f; x3 == 0})
collocate(u == 0)};

% Box constraints
cbox = {2 <= t_f <= 1000
-2 <= collocate(u) <= 2};

% Boundary constraints
cbnd = {initial({x1 == pi/2; x2 == 4; x3 == 0})
final({x2 == 0; x3 == 0})};

% ODEs and path constraints
ceq = collocate({dot(x1) == u
dot(x2) == cos(x1); dot(x3) == sin(x1)});

% Objective
objective = t_f;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Singular Arc';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
x3 = subs(collocate(x3),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.1042 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Singular Arc f_k 4.321198387073171600
sum(|constr|) 0.000000179336713690
f(x_k) + sum(|constr|) 4.321198566409885100
f(x_0) 20.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 77 ConJacEv 77 Iter 70 MinorIter 377
CPU time: 0.483603 sec. Elapsed time: 0.473000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-');
legend('x1','x2','x3');
title('Singular Arc state variables');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Singular Arc control');
</source>

[[File:singularArc_01.png]]

PROPT Simple Bang Bang Problem

2011-11-09T12:12:08Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Function Space Complementarity Methods for Optimal Control Problems, Dissertation, Martin Weiser

==Problem Description==

Find u over t in [-0.5; 0.5 ] to minimize:

<math> J = \int_{-\frac{1}{2}}^{\frac{1}{2}} t*u \mathrm{d}t </math>

subject to:

<math> |u| <= 1 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, -0.5, 1, 20);
setPhase(p);
tomStates x
tomControls u

% Initial guess
x0 = {collocate(u == 1-2*(t+0.5))
icollocate(x == 1-2*(t+0.5))};

% Box constraints
cbox = {-1 <= icollocate(x) <= 1
-1 <= collocate(u) <= 1};

% ODEs and path constraints
ceq = collocate(dot(x) == 0);

% Objective
objective = integrate(t.*u);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Simple Bang Bang Problem';
solution = ezsolve(objective, {cbox, ceq}, x0, options);
t = subs(collocate(t),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lp
Time for symbolic processing: 0.013436 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Simple Bang Bang Problem f_k -0.250490325030180820
sum(|constr|) 0.000000000000870996
f(x_k) + sum(|constr|) -0.250490325029309850
f(x_0) 0.000000000000000000

Solver: CPLEX. EXIT=0. INFORM=1.
CPLEX Dual Simplex LP solver
Optimal solution found

Elapsed time: 0.002000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1);
plot(t,u,'*-');
legend('u');
ylim([-1.1,1.1]);
</source>

[[File:simpleBangBang_01.png]]

PROPT Space Shuttle Reentry

2011-11-09T12:12:03Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

SOCS 6.5.0 Manual

7.4.4 Maximum Crossrange Space Shuttle Reentry Problem.

==Problem Formulation==

Find u over t in [0; t ] to maximize

<math> J = lat </math>

subject to:

The equations given in the code below.

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t t_f

% Scaled time
p1 = tomPhase('p1', t, 0, t_f, 30);
setPhase(p1);

tomStates alt long lat vel ggamma azi
tomControls aalpha bbeta

% Constants
tGuess = 2000;
tMax = 4000;
tMin = 100;
cr2d = 180/pi;
betalim = 90;
weight = 203000;
cm2w = 32.174;
cea = 20902900;
mmu = 0.14076539e17;
rho0 = 0.002378;
href = 23800;
cl0 = -0.20704;
cl1 = 0.029244;
cd0 = 0.07854;
cd1 = -6.1592e-3;
cd2 = 6.21408e-4;
sref = 2690;

alt0 = 260000;
altf = 80000;

vel0 = 25600;
velf = 2500;

% Initial guess
x0 = {
t_f == 1000
icollocate({
alt == alt0-(alt0-altf)*t/t_f
long == -0.5*90/cr2d
lat == -89/cr2d
vel == vel0-(vel0-velf)*t/t_f
ggamma == -1/cr2d-4/cr2d*t/t_f
azi == pi/2-pi*t/t_f
})
collocate({
aalpha == 0
bbeta == 1/cr2d
})
};

% Boundary constraints
cbnd = {
initial({
alt == alt0
long == -0.5*75.3153/cr2d
lat == 0
vel == 25600
ggamma == -1/cr2d
azi == 90/cr2d
aalpha == 17/cr2d
bbeta == -betalim/cr2d
})
final({
alt == altf
vel == velf
ggamma == -5/cr2d
})};

% Box constraints
cbox = {
100 <= t_f <= 5000
0 <= icollocate(alt) <= 300000
-0.5*90/cr2d <= icollocate(long) <= 0.5*90/cr2d
-89/cr2d <= icollocate(lat) <= 89/cr2d
1000 <= icollocate(vel) <= 40000
-89/cr2d <= icollocate(ggamma) <= 89/cr2d
-pi <= icollocate(azi) <= pi
-89/cr2d <= collocate(aalpha) <= 89/cr2d
-betalim/cr2d <= collocate(bbeta) <= 1/cr2d
};

mass = weight/cm2w;
alphad = cr2d*aalpha;
radius = cea+alt;
grav = mmu./radius.^2;

rhodns = rho0*exp(-alt/href);
dynp = 0.5*rhodns.*vel.^2;
subl = cl0+cl1*alphad;
subd = cd0+cd1+cd2*alphad.*alphad;
drag = dynp.*subd*sref;
lift = dynp.*subl*sref;
vrelg = vel./radius-grav./vel;

% ODEs and path constraints
ceq = collocate({
dot(alt) == vel.*sin(ggamma)
dot(long) == vel.*cos(ggamma).*sin(azi)./(radius.*cos(lat))
dot(lat) == vel.*cos(ggamma).*cos(azi)./radius
dot(vel) == -drag./mass-grav.*sin(ggamma)
dot(ggamma) == lift.*cos(bbeta)./(mass.*vel)+cos(ggamma).*vrelg
dot(azi) == lift.*sin(bbeta)./(mass.*vel.*cos(ggamma))+...
vel.*cos(ggamma).*sin(azi).*sin(lat)./(radius.*cos(lat))
});

% Objective
objective = -final(lat)*180/pi;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Shuttle Reentry';
options.Prob.SOL.optPar(30) = 100000;
options.scale = 'auto';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
</source>

<pre>
Problem type appears to be: lpcon
Auto-scaling
Time for symbolic processing: 0.90962 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Shuttle Reentry f_k -15.142919504944096000
sum(|constr|) 0.000000004261906835
f(x_k) + sum(|constr|) -15.142919500682190000
f(x_0) -0.000000000000795808

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 3700 ConJacEv 3698 Iter 424 MinorIter 5688
CPU time: 40.139057 sec. Elapsed time: 10.405000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,3,1)
ezplot(alt)
legend('alt');
title('Shuttle Reentry altitude');

subplot(2,3,2)
ezplot(vel)
legend('vel');
title('Shuttle Reentry velocity');

subplot(2,3,3)
ezplot(long)
legend('long');
title('Shuttle Reentry longitude');

subplot(2,3,4)
ezplot(lat)
legend('lat');
title('Shuttle Reentry latitude');

subplot(2,3,5)
ezplot(aalpha)
legend('aalpha');
title('Shuttle Reentry aalpha');

subplot(2,3,6)
ezplot(bbeta)
legend('bbeta');
title('Shuttle Reentry bbeta');
</source>

[[File:shuttleEntry_01.png]]

PROPT Second Order System

2011-11-09T12:11:59Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Users Guide for dyn.Opt, Example 1

Optimal control of a second order system

End time says 1 in problem text.

==Problem Formulation==

Find u over t in [0; 2 ] to minimize

<math> J = \int_0^{2} u^2/2 \mathrm{d}t </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = u </math>

<math> x_1(0) = 1 </math>

<math> x_1(2) = 0 </math>

<math> x_2(0) = 1 </math>

<math> x_2(2) = 0 </math>

<math> -100 <= u <= 100 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 2, 30);
setPhase(p);

tomStates x1 x2
tomControls u

% Initial guess
x0 = {icollocate({x1 == 1-t/2; x2 == -1+t/2})
collocate(u == -3.5+6*t/2)};

% Box constraints
cbox = {-100 <= icollocate(x1) <= 100
-100 <= icollocate(x2) <= 100
-100 <= collocate(u) <= 100};

% Boundary constraints
cbnd = {initial({x1 == 1; x2 == 1})
final({x1 == 0; x2 == 0})};

% ODEs and path constraints
ceq = collocate({dot(x1) == x2; dot(x2) == u});

% Objective
objective = integrate(u.^2/2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Second Order System';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: qp
Time for symbolic processing: 0.050283 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: 1: Second Order System f_k 3.249999999994048800
sum(|constr|) 0.000000000461823724
f(x_k) + sum(|constr|) 3.250000000455872700
f(x_0) 0.000000000000000000

Solver: CPLEX. EXIT=0. INFORM=1.
CPLEX Barrier QP solver
Optimal solution found

FuncEv 6 GradEv 6 ConstrEv 6 Iter 6
Elapsed time: 0.004000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-');
legend('x1','x2');
title('Second Order System state variables');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Second Order System control');
</source>

[[File:secondOrderSystem_01.png]]

PROPT Satellite Control

2011-11-09T12:11:57Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Users Guide for dyn.Opt, Example 7a, 7b, 7c

==A satellite control problem==

Find T (controls) over t in [0; 100 ] to minimize

7c is free end time

7a:

<math> J = \int_0^{100} ((T1^2+T2^2+T3^2)*0.5) \mathrm{d}t + </math>

<math> w_1^2 + w_2^2 + w_3^2 + (e_1 - 0.70106)^2 + (e_2 - 0.0923)^2 + (e_3 - 0.56098)^2 + (e_4 - 0.43047)^2 </math>

7b:

<math> J = \int_0^{100} ((T1^2+T2^2+T3^2)*0.5) \mathrm{d}t </math>

7c:

<math> J = t_F </math>

subject to:

<math> I1 = 1.0e6 </math>

<math> I2 = 833333.0 </math>

<math> I3 = 916667.0 </math>

<math> T1S = 550 </math>

<math> T2S = 50 </math>

<math> T3S = 550 </math>

<math> \frac{de_1}{dt} = 0.5*(w_1*e_4-w_2*e_3+w_3*e_2) </math>

<math> \frac{de_2}{dt} = 0.5*(w_1*e_3+w_2*e_4-w_3*e_1) </math>

<math> \frac{de_3}{dt} = 0.5*(-w_1*e_2+w_2*e_1+w_3*e_4) </math>

<math> \frac{de_4}{dt} = -0.5*(w_1*e_1+w_2*e_2+w_3*e_3) </math>

<math> \frac{dw_1}{dt} = \frac{(I2-I3)*w_2*w_3+T1*T1S}{I1} </math>

<math> \frac{dw_2}{dt} = \frac{(I3-I1)*w3*w1+T2*T2S}{I2} </math>

<math> \frac{dw_3}{dt} = \frac{(I1-I2)*w1*w2+T3*T3S}{I3} </math>

<math> e(0) = [0 \ 0 \ 0 \ 1] </math>

<math> w(0) = [0.01 \ 0.005 \ 0.001] </math>

7b, 7c - x(100) = [0.70106 0.0923 0.56098 NaN 0 0 0]; 7c - free time 7c - -1 <= T <= 1

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t

% Starting guess
e1opt = 0;
e2opt = 0;
e3opt = 0;
e4opt = 0;
w1opt = 0;
w2opt = 0;
w3opt = 0;
T1opt = 0;
T2opt = 0;
T3opt = 0;

% Final times
tfs = zeros(3,1);

for i=1:3
if i == 3
toms t_f
runs = [10 20 101];
else
runs = [10 40];
end
if i == 2
e1opt = 0; e2opt = 0;
e3opt = 0; e4opt = 0;
w1opt = 0; w2opt = 0;
w3opt = 0; T1opt = 0;
T2opt = 0; T3opt = 0;
end
for n=runs
</source>

<source lang="matlab">
if i == 3
p = tomPhase('p', t, 0, t_f, n);
else
p = tomPhase('p', t, 0, 100, n);
end

setPhase(p);
tomStates e1 e2 e3 e4 w1 w2 w3
tomControls T1 T2 T3

if i == 3
x0 = {t_f == 100};
else
x0 = {};
end
% Initial guess
x0 = {x0; collocate({T1 == T1opt
T2 == T2opt; T3 == T3opt})
icollocate({
e1 == e1opt; e2 == e2opt
e3 == e3opt; e4 == e4opt
w1 == w1opt; w2 == w2opt
w3 == w3opt})};

if i == 3
cbox = {
1 <= t_f <= 1000
-1 <= collocate(T1) <= 1
-1 <= collocate(T2) <= 1
-1 <= collocate(T3) <= 1};
else
cbox = {};
end

% Boundary constraints
cbnd = initial({e1 == 0
e2 == 0; e3 == 0
e4 == 1; w1 == 0.01
w2 == 0.005; w3 == 0.001});
</source>

==Problem 7b and 7c modifications==

<source lang="matlab">
if i ~= 1
cbnd = {cbnd
final({
e1 == 0.70106
e2 == 0.0923
e3 == 0.56098
w1 == 0
w2 == 0
w3 == 0
})};
end

% ODEs and path constraints
I1 = 1.0e6;
I2 = 833333.0;
I3 = 916667.0;
T1Sc = 550;
T2Sc = 50;
T3Sc = 550;
ceq = collocate({
dot(e1) == 0.5*(w1.*e4-w2.*e3+w3.*e2)
dot(e2) == 0.5*(w1.*e3+w2.*e4-w3.*e1)
dot(e3) == 0.5*(-w1.*e2+w2.*e1+w3.*e4)
dot(e4) == -0.5*(w1.*e1+w2.*e2+w3.*e3)
dot(w1) == ((I2-I3)*w2.*w3+T1*T1Sc)/I1
dot(w2) == ((I3-I1)*w3.*w1+T2*T2Sc)/I2
dot(w3) == ((I1-I2)*w1.*w2+T3*T3Sc)/I3});

% Objective
if i == 1
objective = final(w1)^2 + final(w2)^2 + final(w3)^2 +...
(final(e1) - 0.70106)^2 + (final(e2) - 0.0923)^2 + ...
(final(e3) - 0.56098)^2 + (final(e4) - 0.43047)^2 ...
+ integrate((T1.^2+T2.^2+T3.^2)*0.5);
elseif i == 2
objective = integrate((T1.^2+T2.^2+T3.^2)*0.5);
else
objective = t_f;
end
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
if i == 1
options.name = 'Satellite Control 7a';
elseif i == 2
options.name = 'Satellite Control 7b';
else
options.name = 'Satellite Control 7c';
end
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

e1opt = subs(e1, solution);
e2opt = subs(e2, solution);
e3opt = subs(e3, solution);
e4opt = subs(e4, solution);
w1opt = subs(w1, solution);
w2opt = subs(w2, solution);
w3opt = subs(w3, solution);
T1opt = subs(T1, solution);
T2opt = subs(T2, solution);
T3opt = subs(T3, solution);
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.41532 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Satellite Control 7a f_k 0.463944669252504660
sum(|constr|) 0.000000135115456869
f(x_k) + sum(|constr|) 0.463944804367961540
f(x_0) 0.139185999999999230

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 27 ConJacEv 27 Iter 15 MinorIter 105
CPU time: 0.031200 sec. Elapsed time: 0.040000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.4063 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Satellite Control 7a f_k 0.463944366974705650
sum(|constr|) 0.000000000342695654
f(x_k) + sum(|constr|) 0.463944367317401300
f(x_0) -0.536077645550791850

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 3 ConJacEv 3 Iter 2 MinorIter 255
CPU time: 0.078000 sec. Elapsed time: 0.082000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.37013 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Satellite Control 7b f_k 71.411903807059332000
sum(|constr|) 0.000000066784327540
f(x_k) + sum(|constr|) 71.411903873843656000
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 20 ConJacEv 20 Iter 17 MinorIter 205
CPU time: 0.031200 sec. Elapsed time: 0.041000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.37064 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Satellite Control 7b f_k 71.411593978635054000
sum(|constr|) 0.000000015292228483
f(x_k) + sum(|constr|) 71.411593993927283000
f(x_0) 71.236662445283315000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 5 ConJacEv 5 Iter 3 MinorIter 346
CPU time: 0.156001 sec. Elapsed time: 0.155000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.44253 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Satellite Control 7c f_k 99.185331597188878000
sum(|constr|) 0.000010903612748681
f(x_k) + sum(|constr|) 99.185342500801625000
f(x_0) 100.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 8 ConJacEv 8 Iter 6 MinorIter 83
CPU time: 0.031200 sec. Elapsed time: 0.024000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.4436 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Satellite Control 7c f_k 98.945461990937730000
sum(|constr|) 0.000017743233813344
f(x_k) + sum(|constr|) 98.945479734171542000
f(x_0) 100.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 5 ConJacEv 5 Iter 4 MinorIter 154
CPU time: 0.046800 sec. Elapsed time: 0.041000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.44408 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Satellite Control 7c f_k 98.834204968061897000
sum(|constr|) 0.000012503591248200
f(x_k) + sum(|constr|) 98.834217471653147000
f(x_0) 100.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 7 ConJacEv 7 Iter 6 MinorIter 972
CPU time: 5.756437 sec. Elapsed time: 2.873000 sec.

</pre>

<source lang="matlab">
end
tfs(i) = subs(final(t),solution);
% We only want to plot the solution from the first problem
if i == 1
tp = subs(collocate(t),solution);
e1p = subs(collocate(e1),solution);
e2p = subs(collocate(e2),solution);
e3p = subs(collocate(e3),solution);
e4p = subs(collocate(e4),solution);
w1p = subs(collocate(w1),solution);
w2p = subs(collocate(w2),solution);
w3p = subs(collocate(w3),solution);
T1p = subs(collocate(T1),solution);
T2p = subs(collocate(T2),solution);
T3p = subs(collocate(T3),solution);
end
end

disp(sprintf('\nFinal time for 7a = %1.4g',tfs(1)));
disp(sprintf('\nFinal time for 7b = %1.4g',tfs(2)));
disp(sprintf('\nFinal time for 7c = %1.4g',tfs(3)));
</source>

<pre>

Final time for 7a = 100

Final time for 7b = 100

Final time for 7c = 98.83

</pre>

==Plot result==

<source lang="matlab">
subplot(3,1,1)
plot(tp,e1p,'*-',tp,e2p,'*-',tp,e3p,'*-',tp,e4p,'*-');
legend('e1','e2','e3','e4');
title('Satellite Control state variables (e)');

subplot(3,1,2)
plot(tp,w1p,'*-',tp,w2p,'*-',tp,w3p);
legend('w1','w2','w3');
title('Satellite Control state variables (w)');

subplot(3,1,3)
plot(tp,T1p,'+-',tp,T2p,'+-',tp,T3p,'+-');
legend('T1','T2','T3');
title('Satellite Control controls');
</source>

[[File:satelliteControl_01.png]]

PROPT Time-optimal Trajectories for Robot Manipulators

2011-11-09T12:11:51Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Users Guide for dyn.Opt, Example 2

Dissanayake, M., Goh, C. J., & Phan-Thien, N., Time-optimal Trajectories for Robot Manipulators, Robotica, Vol. 9, pp. 131-138, 1991.

==Problem Formulation==

Find u over t in [0; t_F ] to minimize

<math> J = t_F </math>

subject to:

<math> x(0) = [0 \ -2 \ 0 \ 0] </math>

<math> x(t_F) = [1 \ -1 \ 0 \ 0] </math>

<math> L_1 = 0.4; </math>

<math> L_2 = 0.4; </math>

<math> m_1 = 0.5; </math>

<math> m_2 = 0.5; </math>

<math> Eye_1 = 0.1; </math>

<math> Eye_2 = 0.1; </math>

<math> el_1 = 0.2; </math>

<math> el_2 = 0.2; </math>

<math> cos(x2) = cos(x_2); </math>

<math> H_{11} = Eye_1 + Eye_2 + m_1*el_1^2 + m_2*(L_1^2+el_2^2+2.0*L_1*el_2*cos(x2)); </math>

<math> H_{12} = Eye_2 + m_2*el_2^2 + m_2*L_1*el_2*cos(x2); </math>

<math> H_{22} = Eye_2 + m_2*el_2^2; </math>

<math> h = m_2*L_1*el_2*sin(x_2); </math>

<math> delta = \frac{1.0}{H_{11}*H_{22}-H_{12}*H_{12}}; </math>

<math> \frac{dx_1}{dt} = x_3 </math>

<math> \frac{dx_2}{dt} = x_4 </math>

<math> \frac{dx_3}{dt} = delta*(2.0*h*H_{22}*x_3*x_4+h*H_{22}*x_4^2 + h*H_{12}*x_3^2+H_{22}*u_1-H_{12}*u_2); </math>

<math> \frac{dx_4}{dt} = delta*(-2.0*h*H_{12}*x_3*x_4-h*H_{11}*x_3^2- h*H_{12}*x_4^2+H_{11}*u_2-H_{12}*u_1); </math>

<math> -10 <= u <= 10 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
toms t_f

tfopt = 7;
x1opt = 1*t/t_f;
x2opt = -2+1*t/t_f;
x3opt = 2;
x4opt = 4;
u1opt = 10-20*t/t_f;
u2opt = -10+20*t/t_f;
</source>

==Solve the problem, using a successively larger number collocation points==

<source lang="matlab">
for n=[30 60]
</source>

<source lang="matlab">
p = tomPhase('p', t, 0, t_f, n);
setPhase(p);

tomStates x1 x2 x3 x4
tomControls u1 u2

% Initial guess
x0 = {t_f == tfopt
icollocate({
x1 == x1opt
x2 == x2opt
x3 == x3opt
x4 == x4opt})
collocate({
u1 == u1opt
u2 == u2opt})};

% Box constraints
cbox = {
0.1 <= t_f <= 50
-10 <= collocate(u1) <= 10
-10 <= collocate(u2) <= 10};

% Boundary constraints
cbnd = {initial({x1 == 0; x2 == -2
x3 == 0; x4 == 0})
final({x1 == 1; x2 == -1
x3 == 0; x4 == 0})};

% ODEs and path constraints
L_1 = 0.4; L_2 = 0.4;
m_1 = 0.5; m_2 = 0.5;
Eye_1 = 0.1; Eye_2 = 0.1;
el_1 = 0.2; el_2 = 0.2;

H_11 = Eye_1 + Eye_2 + m_1*el_1^2+ ...
m_2*(L_1^2+el_2^2+2.0*L_1*el_2*cos(x2));
H_12 = Eye_2 + m_2*el_2^2 + m_2*L_1*el_2*cos(x2);
H_22 = Eye_2 + m_2*el_2^2;
h = m_2*L_1*el_2*sin(x2);
delta = 1.0./(H_11.*H_22-H_12.^2);

ceq = collocate({
dot(x1) == x3
dot(x2) == x4
dot(x3) == delta.*(2.0*h.*H_22.*x3.*x4 ...
+h.*H_22.*x4.^2+h.*H_12.*x3.^2+H_22.*u1-H_12.*u2)
dot(x4) == delta.*(-2.0*h.*H_12.*x3.*x4 ...
-h.*H_11.*x3.^2-h.*H_12.*x4.^2+H_11.*u2-H_12.*u1)});

% Objective
objective = t_f;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Robot Manipulators';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

% Optimal x, y, and speed, to use as starting guess
% in the next iteration
tfopt = subs(final(t), solution);
x1opt = subs(x1, solution);
x2opt = subs(x2, solution);
x3opt = subs(x3, solution);
x4opt = subs(x4, solution);
u1opt = subs(u1, solution);
u2opt = subs(u2, solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.59011 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Robot Manipulators f_k 0.391698237386178090
sum(|constr|) 0.000063099633743277
f(x_k) + sum(|constr|) 0.391761337019921390
f(x_0) 7.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 68 ConJacEv 68 Iter 26 MinorIter 437
CPU time: 0.124801 sec. Elapsed time: 0.134000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.58064 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Robot Manipulators f_k 0.391820155673057110
sum(|constr|) 0.000000000018899679
f(x_k) + sum(|constr|) 0.391820155691956770
f(x_0) 0.391698237386178090

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 16 ConJacEv 16 Iter 12 MinorIter 440
CPU time: 0.187201 sec. Elapsed time: 0.191000 sec.

</pre>

<source lang="matlab">
end

t = subs(collocate(t),solution);
x1 = subs(collocate(x1opt),solution);
x2 = subs(collocate(x2opt),solution);
x3 = subs(collocate(x3opt),solution);
x4 = subs(collocate(x4opt),solution);
u1 = subs(collocate(u1opt),solution);
u2 = subs(collocate(u2opt),solution);
</source>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Robot Manipulators state variables');

subplot(2,1,2)
plot(t,u1,'+-',t,u2,'+-');
legend('u1','u2');
title('Robot Manipulators control');
</source>

[[File:robotManipulators_01.png]]

PROPT Robot Arm Movement

2011-11-09T12:11:49Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Benchmarking Optimization Software with COPS Elizabeth D. Dolan and Jorge J. More ARGONNE NATIONAL LABORATORY

==Problem Formulation==

Find u(t) over t in [0; 1 ] to minimize

<math> J = t_f </math>

subject to:

<math> L*\frac{d^{2}rho}{dt^{2}} = u_1 </math>

<math> I_1*\frac{d^{2}theta}{dt^{2}} = u_2 </math>

<math> I_2*\frac{d^{2}phi}{dt^{2}} = u_3 </math>

<math> 0 <= rho <= L </math>

<math> |theta| <= pi </math>

<math> 0 <= phi <= pi </math>

<math> |u_{1:3}| <= 1 </math>

<math> I_1 = \frac{((L-rho)^3+rho^3)}{3}*sin(phi)^2 </math>

<math> I_2 = \frac{((L-rho)^3+rho^3)}{3} </math>

The boundary conditions are:

<math> [rho_0 \ theta_0 \ phi_0] = [4.5 \ 0 \ \frac{pi}{4}] </math>

<math> [rho_1 \ theta_1 \ phi_1] = [4.5 \ \frac{2*pi}{3} \ \frac{pi}{4}] </math>

<math> L = 5 </math>

All first order derivatives are 0 at boundaries.

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
toms t_f

% Initial guess
tfopt = 1;
x1opt = 4.5;
x2opt = 0;
x3opt = 2*pi/3*t.^2;
x4opt = 0;
x5opt = pi/4;
x6opt = 0;
u1opt = 0;
u2opt = 0;
u3opt = 0;

for n=[20 100]
</source>

<source lang="matlab">
%rho d(rho)dt theta d(theta)dt phi d(phi)dt
p = tomPhase('p', t, 0, t_f, n);
setPhase(p);
tomStates x1 x2 x3 x4 x5 x6
tomControls u1 u2 u3

% Initial guess
x0 = {t_f == tfopt
icollocate({x1 == x1opt
x2 == x2opt; x3 == x3opt
x4 == x4opt; x5 == x5opt
x6 == x6opt})
collocate({u1 == u1opt
u2 == u2opt; u3 == u3opt})};

% Box constraints
L = 5;
cbox = {
0.1 <= t_f <= 10
0 <= icollocate(x1) <= L
-pi <= icollocate(x3) <= pi
0 <= icollocate(x5) <= pi
-1 <= collocate(u1) <= 1
-1 <= collocate(u2) <= 1
-1 <= collocate(u3) <= 1};

% Boundary constraints
cbnd = {initial({x1 == 4.5; x2 == 0
x3 == 0; x4 == 0
x5 == pi/4; x6 == 0})
final({x1 == 4.5; x2 == 0
x3 == 2*pi/3
x4 == 0
x5 == pi/4
x6 == 0
})};

I1 = ((L-x1).^3+x1.^3)./3.*sin(x5).^2;
I2 = ((L-x1).^3+x1.^3)/3;

% ODEs and path constraints
ceq = collocate({dot(x1) == x2
dot(x2) == u1/L; dot(x3) == x4
dot(x4) == u2./I1; dot(x5) == x6
dot(x6) == u3./I2});

% Objective
objective = t_f;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Robot Arm';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

% Optimal x, y, and speed, to use as starting guess in the next iteration
x1opt = subs(x1, solution);
x2opt = subs(x2, solution);
x3opt = subs(x3, solution);
x4opt = subs(x4, solution);
x5opt = subs(x5, solution);
u1opt = subs(u1, solution);
u2opt = subs(u2, solution);
u3opt = subs(u3, solution);
tfopt = subs(final(t), solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.3121 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Robot Arm f_k 9.146367546080650600
sum(|constr|) 0.000000278037569148
f(x_k) + sum(|constr|) 9.146367824118220500
f(x_0) 1.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 17 ConJacEv 17 Iter 11 MinorIter 225
CPU time: 0.031200 sec. Elapsed time: 0.037000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.3096 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Robot Arm f_k 9.140854009735493300
sum(|constr|) 0.000002639179211847
f(x_k) + sum(|constr|) 9.140856648914704500
f(x_0) 9.146367546080650600

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 7 ConJacEv 7 Iter 4 MinorIter 729
CPU time: 1.606810 sec. Elapsed time: 0.634000 sec.

</pre>

<source lang="matlab">
end

t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x3 = subs(collocate(x3),solution);
x5 = subs(collocate(x5),solution);
u1 = subs(collocate(u1),solution);
u2 = subs(collocate(u2),solution);
u3 = subs(collocate(u3),solution);
</source>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x3,'*-',t,x5,'*-');
legend('rho','theta','phi');
title('Robot Arm state variables');

subplot(2,1,2)
plot(t,u1,'*-',t,u2,'*-',t,u3,'*-');
legend('u1','u2','u3');
title('Robot Arm control variables');
</source>

[[File:robotArmMovement_01.png]]

PROPT Rigid Body Rotation

2011-11-09T12:11:48Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

On smooth optimal control determination, Ilya Ioslovich and Per-Olof Gutman, Technion, Israel Institute of Technology.

Example 1: Rigid body rotation

==Problem Description==

Find u over t in [0; 1 ] to minimize:

<math> J = \frac{1}{4}*\int_0^{1} (u_1^2+u_2^2)^2 \mathrm{d}t </math>

subject to:

<math> \frac{dx}{dt} = a*y+u_1 </math>

<math> \frac{dy}{dt} = -a*x+u_2 </math>

<math> \frac{du_1}{dt} = a*u_2 </math>

<math> \frac{du_2}{dt} = -a*u_1 </math>

<math> x(t_0) = [0.9 \ 0.75] </math>

<math> x(t_f) = [0 \ 0] </math>

<math> a = 2 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 1, 20);
setPhase(p);

tomStates x y u1 u2

% Boundary constraints
cbnd = {initial({x == 0.9; y == 0.75})
final({x == 0; y == 0})};

% ODEs and path constraints
a = 2;
ceq = collocate({dot(x) == a*y+u1; dot(y) == -a*x+u2
dot(u1) == a*u2; dot(u2) == -a*u1});

% Objective
objective = 0.25*integrate((u1.^2+u2.^2).^2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Rigid Body Rotation';
solution = ezsolve(objective, {cbnd, ceq}, [], options);
t = subs(collocate(t),solution);
x = subs(collocate(x),solution);
y = subs(collocate(y),solution);
u1 = subs(collocate(u1),solution);
u2 = subs(collocate(u2),solution);
</source>

<pre>
Problem type appears to be: con
Time for symbolic processing: 0.11175 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Rigid Body Rotation f_k 0.470939062500258190
sum(|constr|) 0.000000000003070916
f(x_k) + sum(|constr|) 0.470939062503329120
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 3 GradEv 1 MinorIter 39
Elapsed time: 0.010000 sec.

</pre>

==Plot result==

<source lang="matlab">
figure(1);
subplot(2,1,1);
plot(t,x,'*-',t,y,'*-');
legend('x','y');

subplot(2,1,2);
plot(t,u1,'*-',t,u2,'*-');
legend('u1','u2');
</source>

[[File:rigidBodyRotation_01.png]]

PROPT Rayleigh Unconstrained

2011-11-09T12:11:44Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Lecture Notes for ECE/MAE 7360, Robust and Optimal Control (part 2) Fall 2003, Jinsong Liang, Nov. 20, 2003 Utah State University at Logan

==Problem Formulation==

Find u over t in [0; t_F ] to minimize

<math> J = \int_0^{2.5} x_1^2 + u^2 \mathrm{d}t </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = -x_1+(1.4-0.14*x_2^2)*x_2+4*u </math>

<math> x_1(0) = -5 </math>

<math> x_2(0) = -5 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 2.5, 50);
setPhase(p);

tomStates x1 x2
tomControls u

% Initial guess
x0 = {icollocate({x1 == -5; x2 == -5})
collocate(u == 0)};

% Box constraints
cbox = {-100 <= icollocate(x1) <= 100
-100 <= icollocate(x2) <= 100
-100 <= collocate(u) <= 100};

% Boundary constraints
cbnd = initial({x1 == -5; x2 == -5});

% ODEs and path constraints
ceq = collocate({dot(x1) == x2
dot(x2) == -x1+(1.4-0.14*x2.^2).*x2+4*u});

% Objective
objective = integrate(x1.^2+u.^2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Rayleigh Unconstrained';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.10859 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Rayleigh Unconstrained f_k 29.376079656023276000
sum(|constr|) 0.000000003739617241
f(x_k) + sum(|constr|) 29.376079659762894000
f(x_0) 62.499999999999517000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 47 ConJacEv 47 Iter 36 MinorIter 144
CPU time: 0.093601 sec. Elapsed time: 0.091000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-');
legend('x1','x2');
title('Rayleigh Unconstrained state variables');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Rayleigh Unconstrained System control');
</source>

[[File:rayleighUnconstrained_01.png]]

PROPT Quadruple Integral

2011-11-09T12:11:31Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

ITERATIVE DYNAMIC PROGRAMMING, REIN LUUS

12.4.5 Example 5

CHAPMAN & HALL/CRC Monographs and Surveys in Pure and Applied Mathematics

==Problem Formulation==

Find u over t in [0; t_F ] to minimize

<math> J = t_F </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = x_3 </math>

<math> \frac{dx_3}{dt} = x_4 </math>

<math> \frac{dx_4}{dt} = u </math>

The initial condition are:

<math> x(0) = [0.1 \ 0.2 \ 0.3 \ 0] </math>

<math> x(t_F) = [0 \ 0 \ 0 \ 0] </math>

<math> -1 <= u <= 1 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
toms t_f
p = tomPhase('p', t, 0, t_f, 30);
setPhase(p);

tomStates x1 x2 x3 x4
tomControls u

% Initial guess
x0 = {t_f == 5
icollocate({x1 == 0.1-0.1*t/t_f
x2 == 0.2-0.2*t/t_f; x3 == 0.3-0.3*t/t_f})
collocate(u == -1)};

% Box constraints
cbox = {0.1 <= t_f <= 100
-1 <= collocate(u) <= 1};

% Boundary constraints
cbnd = {initial({x1 == 0.1; x2 == 0.2; x3 == 0.3; x4 == 0})
final({x1 == 0; x2 == 0; x3 == 0; x4 == 0})};

% ODEs and path constraints
ceq = collocate({
dot(x1) == x2; dot(x2) == x3
dot(x3) == x4; dot(x4) == u});

% Objective
objective = t_f;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Quadruple Integral';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
x3 = subs(collocate(x3),solution);
x4 = subs(collocate(x4),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.12111 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Quadruple Integral f_k 4.849598187644263100
sum(|constr|) 0.000000078353877264
f(x_k) + sum(|constr|) 4.849598265998140300
f(x_0) 5.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 58 ConJacEv 58 Iter 31 MinorIter 474
CPU time: 0.109201 sec. Elapsed time: 0.090000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Quadruple Integral state variables');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Quadruple Integral control');
</source>

[[File:quadrupleIntegral_01.png]]

PROPT Quadratic constraint problem

2011-11-09T12:11:30Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Paper: LINEAR-QUADRATIC OPTIMAL CONTROL WITH INTEGRAL QUADRATIC CONSTRAINTS. OPTIMAL CONTROL APPLICATIONS AND METHODS Optim. Control Appl. Meth., 20, 79-92 (1999)

E. B. LIM(1), Y. Q. LIU(2), K. L. TEO(2) AND J. B. MOORE(1)

(1) Department of Systems Engineering, Research School of Information Sciences and Engineering, Australian National University, Canberra ACT 0200, Australia

(2) School of Mathematics and Statistics, Curtin University of Technology, Perth, WA 6845, Australia

==Problem Formulation==

Find u(t) over t in [0; 1 ] to minimize

<math> J = 0.5*x_1(1)^2 + 0.5*\int_0^1 (x_1^2 + u_1^2 + u_2^2) \mathrm{d}t </math>

subject to:

<math> \frac{dx_1}{dt} = 3*x_1+x_2 + u_1 </math>

<math> \frac{dx_2}{dt} = -x_1+2*x_2 + u_2 </math>

<math> x_1(0) = 4 </math>

<math> x_2(0) = -4 </math>

<math> 0.5*x_2(1)^2 + 0.5*\int_0^1 (x_1^2 + u_1^2 + u_2^2) <= 80 </math>

Introduce a new variable to remove integral in constraint:

<math> \frac{dx_3}{dt} = 0.5 * (x_1.^2 + u_1.^2 + u_2.^2) </math>

resulting in event constraint:

<math> 0.5*x_2(1)^2 + x_3(1) <= 8 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 1, 50);
setPhase(p);

tomStates x1 x2 x3
tomControls u1 u2

% Initial guess
x0 = {icollocate({
x1 == 4-5*t
x2 == -4-1*t
x3 == 50*t
})
collocate({
u1 == -10+10*t
u2 == 14-12*t})};

% Boundary constraints
cbnd = {
initial({
x1 == 4
x2 == -4
x3 == 0
})
final(x2)^2/2+final(x3) <= 80};

% ODEs and path constraints
ceq = collocate({
dot(x1) == 3*x1+x2 + u1
dot(x2) == -x1+2*x2 + u2
dot(x3) == 1/2 * (x2.^2 + u1.^2 + u2.^2)
});

% Objective
objective = final(x1)^2/2 + final(x3);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Quadratic Constraint';
solution = ezsolve(objective, {cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
x3 = subs(collocate(x3),solution);
u1 = subs(collocate(u1),solution);
u2 = subs(collocate(u2),solution);
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.12268 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Quadratic Constraint f_k 67.888740121887480000
sum(|constr|) 0.000000192408232753
f(x_k) + sum(|constr|) 67.888740314295717000
f(x_0) 50.499999999999915000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 32 ConJacEv 32 Iter 31 MinorIter 279
CPU time: 0.218401 sec. Elapsed time: 0.217000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3/10,'*-');
legend('x1','x2','x3/10');
title('Quadratic Constraint state variables');

subplot(2,1,2)
plot(t,u1,'+-',t,u2,'+-');
legend('u1','u2');
title('Quadratic Constraint control');
</source>

[[File:quadraticConstraint_01.png]]

PROPT Plug-Flow Tubular Reactor

2011-11-09T12:11:26Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

A HYBRID METHOD FOR THE OPTIMAL CONTROL OF CHEMICAL PROCESSES 1998, E F Carrasco, J R Banga

Case Study II: Plug-Flow Tubular Reactor

==Problem description==

This case study considers a plug-flow reactor as studied by Reddy and Husain, Luus and Mekarapiruk and Luus. The objective is to maximize the normalized concentration of the desired product.

Find u(t) to maximize

<math> J = x_1(t_f) </math>

subject to:

<math> \frac{dx_1}{dt} = (1-x_1)*k_1-x_1*k_2 </math>

<math> \frac{dx_2}{dt} = 300*((1-x_1)*k_1-x_1*k_2)-u*(x_2-290) </math>

where x1 denotes the normalized concentration of de desired product, and x2 is the temperature. The initial conditions are:

<math> x(t_0) = [0 \ 380]' </math>

The rate constants are given by:

<math> k_1 = 1.7536e5*exp(-\frac{1.1374e4}{1.9872*x_2}) </math>

<math> k_2 = 2.4885e10*exp(-\frac{2.2748e4}{1.9872*x_2}) </math>

where the final time t_f = 5 min. The constraint on the control variable (the normalized coolant flow rate) is:

<math> 0 <= u <= 0.5 </math>

In addition, there is an upper path constraint on the temperature:

<math> x_2(t) <= 460 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 5, 30);
setPhase(p);

tomStates x1 x2
tomControls u

% Initial guess
x0 = {icollocate({x1 == 0.6*t/5
x2 == 380})
collocate(u == 0.25)};

% Box constraints
cbox = {0 <= icollocate(x1) <= 10
100 <= icollocate(x2) <= 460
0 <= collocate(u) <= 0.5};

% Boundary constraints
cbnd = initial({x1 == 0; x2 == 380});

% ODEs and path constraints
k1 = 1.7536e5*exp(-1.1374e4/1.9872./x2);
k2 = 2.4885e10*exp(-2.2748e4/1.9872./x2);
ceq = collocate({
dot(x1) == (1-x1).*k1-x1.*k2
dot(x2) == 300*((1-x1).*k1-x1.*k2)-u.*(x2-290)});

% Objective
objective = -final(x1);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Plug Flow Reactor';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.18266 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Plug Flow Reactor f_k -0.677125737913608190
sum(|constr|) 0.000031723901907782
f(x_k) + sum(|constr|) -0.677094014011700440
f(x_0) -0.599999999999999200

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 346 ConJacEv 346 Iter 112 MinorIter 325
CPU time: 0.234001 sec. Elapsed time: 0.227000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2/100,'*-');
legend('x1','x2/100');
title('Plug Flow Reactor state variables');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Plug Flow Reactor control');
</source>

[[File:plugFlowReactor_01.png]]

PROPT Penicillin Plant

2011-11-09T12:11:21Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Fed-batch Fermentor Control: Dynamic Optimization of Batch Processes II. Role of Measurements in Handling Uncertainty 2001, B. Srinivasan, D. Bonvin, E. Visser, S. Palanki

Illustrative example: Nominal Optimization of a Fed-Batch Fermentor for Penicillin Production.

==Problem description==

This particular example was featured in the work of B. Srinivasan et al. 2001. The optimal trajectories for the problem was provided in the work.

In this problem, the objective is to maximize the concentration of penicillin, P, produced in a fed-batch bioreactor, given a finite amount of time.

<pre>
Reactions: S -> X, S -> P
Conditions: Fed-batch, isothermal.
Objective: Maximize the concentration of product P at a given final time.
Manipulated variable: Feed rate of S.
Constraints: Input bounds; upper limit on the biomass concentration,
which is motivated by oxygen-transfer limitation typically occurring at
large biomass concentration.
</pre>

<math> J = -P(t_f) </math>

subject to:

<math> \frac{X}{dt} = my*X - \frac{u}{V}*X </math>

<math> \frac{S}{dt} = -\frac{my*X}{Yx} - \frac{v*X}{Yp} + \frac{u}{V}*(Sin-S) </math>

<math> \frac{P}{dt} = v*X - \frac{u}{V}*P </math>

<math> \frac{V}{dt} = u </math>

<math> my = \frac{um*S}{Km+S+S^{2}/Ki} </math>

Programmer: Wee Kiat Lim (Nanyang Technological University)

<source lang="matlab">
% Copyright (c) 2009-2009 by Tomlab Optimization Inc.
</source>

==Problem setup==

Penalty for variations in u

<source lang="matlab">
penalty_constant = 0.001;

% Various constants
miu_m = 0.02; Km = 0.05; Ki = 5;
Yx = 0.5; Yp = 1.2; v = 0.004;
Sin = 200; umin = 0; umax = 1;
Xmin = 0; Xmax = 3.7; Smin = 0;

% no. of collocation points to use
narr = [20 80];
for n=narr
toms t1
toms tcut
p1 = tomPhase('p1', t1, 0, tcut, n);

setPhase(p1);
tomStates X1 S1 P1 V1 %Vs %Scaling is disabled here
tomControls u1

% Initial guess
if n == narr(1)
x01 = {tcut == 75
icollocate({X1 == 1+2.7*t1/tcut; S1 == 0.5;
P1 == 0.6*t1/tcut; V1 == 150})
collocate(u1 == 0.03+0.06*t1/tcut)};
else
x01 = {tcut == tcutg
icollocate({X1 == Xg1; S1 == Sg1; P1 == Pg1; V1 == Vg1})
collocate(u1 == ug1)};
end

% Box constraints
cbox1 = {75 <= tcut <= 85
0 <= icollocate(X1) <= Xmax
Smin <= icollocate(S1) <= 100
0 <= icollocate(P1) <= 5
1 <= icollocate(V1) <= 300
umin <= collocate(u1) <= umax};

% Boundary constraints
cbnd1 = initial({X1 == 1; S1 == 0.5;
P1 == 0; V1 == 150});

miu1 = (miu_m*S1)/(Km + S1 + S1^2/Ki);

% ODEs and path constraints
temp11 = miu1*X1;
temp21 = u1/V1;
temp31 = v*X1;
ceq1 = collocate ({
dot(X1) == temp11 - u1/V1*X1
dot(S1) == -temp11/Yx - temp31/Yp + temp21*(Sin - S1)
dot(P1) == temp31 - temp21*P1
dot(V1) == u1});

if n == narr(1)
% No objective in first phase
objective = 0;
else
% Variation penalty
objective = penalty_constant*integrate(dot(u1)^2);
end

toms t2
p2 = tomPhase('p2', t2, tcut, 150-tcut, n);

setPhase(p2);
tomStates X2 S2 P2 V2 %Vs %Scaling is disabled here
tomControls u2

% Initial guess
if n == narr(1)
x02 = {
icollocate({X2 == Xmax; S2 == 0; P2 == 0.6+t2/150; V2 == 150});
collocate(u2 == 0.01);
};
else
x02 = {
icollocate({X2 == Xg2; S2 == Sg2; P2 == Pg2; V2 == Vg2})
collocate(u2 == ug2)
};
end

% Box constraints
umax2 = 0.03;
cbox2 = {0 <= icollocate(X2) <= Xmax
Smin <= icollocate(S2) <= 100
0 <= icollocate(P2) <= 5
1 <= icollocate(V2) <= 300
umin <= collocate(u2) <= umax2
initial(S2) <= 0.2};

miu2 = (miu_m*S2)/(Km + S2 + S2^2/Ki);

% ODEs and path constraints
temp12 = miu2*X2;
temp22 = u2/V2;
temp32 = v*X2;
ceq2 = collocate ({
dot(X2) == temp12 - u2/V2*X2
dot(S2) == -temp12/Yx - temp32/Yp + temp22*(Sin - S2)
dot(P2) == temp32 - temp22*P2
dot(V2) == u2});

% Phase links
links = {initial(X2) == final(p1,X1)
initial(S2) == final(p1,S1)
initial(P2) == final(p1,P1)
initial(V2) == final(p1,V1)};

if n == narr(1)
% Objective (Negative sign is added to 'maximize' P)
objective = -final(P2);
ptype = 'lpcon';
solver = 'snopt';
else
objective = objective-final(P2)+penalty_constant*integrate(dot(u2)^2);
ptype = 'con';
solver = 'snopt';
end

% Solve the problem
options = struct;
options.name = 'Penicillin Plant';
Prob = sym2prob(ptype, objective, {cbox1, cbnd1, ceq1, cbox2, ceq2, links}, {x01, x02}, options);
Result = tomRun(solver, Prob, 1);
solution = getSolution(Result);

ug1 = subs(u1,solution);
Xg1 = subs(X1,solution);
Sg1 = subs(S1,solution);
Pg1 = subs(P1,solution);
Vg1 = subs(V1,solution);
ug2 = subs(u2,solution);
Xg2 = subs(X2,solution);
Sg2 = subs(S2,solution);
Pg2 = subs(P2,solution);
Vg2 = subs(V2,solution);
tcutg = solution.tcut;
end
</source>

<pre>
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Penicillin Plant f_k -1.682729742163946000
sum(|constr|) 0.000005098522916677
f(x_k) + sum(|constr|) -1.682724643641029200
f(x_0) -1.599999999999999600

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 652 ConJacEv 652 Iter 210 MinorIter 2952
CPU time: 1.326009 sec. Elapsed time: 1.312000 sec.
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Penicillin Plant f_k -1.682693890525556800
sum(|constr|) 0.000001881327092638
f(x_k) + sum(|constr|) -1.682692009198464300
f(x_0) -1.682727201466558900

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 18 GradEv 16 ConstrEv 16 ConJacEv 16 Iter 15 MinorIter 976
CPU time: 1.419609 sec. Elapsed time: 1.419000 sec.

</pre>

==Plot result==

Optimal states and control trajectories

<source lang="matlab">
uopt = subs([collocate(p1,u1);collocate(p2,u2)],solution);
Xopt = subs([collocate(p1,X1);collocate(p2,X2)],solution);
Sopt = subs([collocate(p1,S1);collocate(p2,S2)],solution);
Popt = subs([collocate(p1,P1);collocate(p2,P2)],solution);
Vopt = subs([collocate(p1,V1);collocate(p2,V2)],solution);

t = subs([collocate(p1,t1);collocate(p2,t2)],solution);
np = length(t);

Pfinal=subs(final(p2,P2),solution);

% Plots of the trajectories
figure(1)
subplot(3,1,1);
plot(t,Popt,'*-');
title(['Final Penicillin concentration is ',num2str(Pfinal),' g/L.'])
ylabel('Penicillin Conc')
xlabel('Time (hrs)')
subplot(3,1,2);
plot(t,Xopt,'*-');
ylabel('Cell Mass Conc')
xlabel('Time (hrs)')
subplot(3,1,3);
plot(t,Sopt,'*-');
ylabel('Substrate Conc')
xlabel('Time (hrs)')

figure(2)
subplot(2,1,1);
plot(t,Vopt,'*-');
title(['Final Penicillin concentration is ',num2str(Pfinal),' g/L.'])
ylabel('Volume of medium')
xlabel('Time (hrs)')
subplot(2,1,2);
plot(t,uopt,'*-');
ylabel('Feed flowrate')
xlabel('Time (hrs)')

fprintf('\n')
fprintf('Optimization completed... \n')
fprintf('Final Penicillin concentration is %5.4f g/L.\n',Pfinal)
</source>

<pre>

Optimization completed...
Final Penicillin concentration is 1.6827 g/L.

</pre>

[[File:penicillinPlant_01.png]]

[[File:penicillinPlant_02.png]]

PROPT Pendulum Gravity Estimation

2011-11-09T12:11:18Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

User's Guide for DIRCOL

Problem 2.5 Pendulum

==Problem Formulation==

Find u over t in [0; 2 ] to minimize

<math> J = \frac{1}{2}*((x_{t_f}-x_{f_{meas}})^2 + (y_{t_f}-y_{f_{meas}})^2) </math>

subject to:

<math> \frac{dx}{dt} = u </math>

<math> \frac{dy}{dt} = v </math>

<math> \frac{du}{dt} = lambda*x/m </math>

<math> \frac{dv}{dt} = lambda*y/m-g </math>

<math> x^2 + y^2 - L^2 = 0 </math>

<math> [x_0 \ y_0 \ u_0 \ v_0] = [0.4 \ -0.3 \ 0 \ 0] </math>

<math> [x_{f_{meas}} \ y_{f_{meas}}] = [-0.231625 \ -0.443109] </math>

<math> L = 0.5 </math>

<math> m = 0.3 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t g

% Initial guess
gopt = 20;
xopt = 0.4-(0.4+0.231625)*t/2;
yopt = -0.3-(-0.3+0.443109)*t/2;
uopt = 0;
vopt = 0;
lambdaopt = -5;

for n=[20 51]
</source>

<source lang="matlab">
p = tomPhase('p', t, 0, 2, n);
setPhase(p);
tomStates x y u v
tomControls lambda

% Initial guess
x0 = {g == gopt
icollocate({
x == xopt
y == yopt
u == uopt
v == vopt})
collocate(lambda == lambdaopt)};

% Box constraints
cbox = {1 <= g <= 100};

% Boundary constraints
cbnd = initial({x == 0.4; y == -0.3
u == 0; v == 0});

L = 0.5;
m = 0.3;
xmeas = -0.231625;
ymeas = -0.443109;

% ODEs and path constraints
ceq = collocate({
dot(x) == u
dot(y) == v
dot(u) == lambda.*x/m
dot(v) == lambda.*y/m-g
x.^2 + y.^2 - L^2 == 0});

% Objective
objective = 1/2*((final(x)-xmeas)^2+(final(y)-ymeas)^2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Pendulum Gravity';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
gopt = subs(g, solution);
xopt = subs(x, solution);
yopt = subs(y, solution);
uopt = subs(u, solution);
vopt = subs(v, solution);
lambdaopt = subs(lambda, solution);
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.14939 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Pendulum Gravity f_k 0.000000017178935982
sum(|constr|) 0.000001384358823592
f(x_k) + sum(|constr|) 0.000001401537759574
f(x_0) -0.124997863253000020

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 98 ConJacEv 98 Iter 43 MinorIter 96
CPU time: 0.078000 sec. Elapsed time: 0.088000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.15159 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Pendulum Gravity f_k 0.000000000017699578
sum(|constr|) 0.000000032138058518
f(x_k) + sum(|constr|) 0.000000032155758096
f(x_0) -0.124997846074064280

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 13 ConJacEv 13 Iter 12 MinorIter 190
CPU time: 0.202801 sec. Elapsed time: 0.200000 sec.

</pre>

<source lang="matlab">
end
</source>

==Show result==

<source lang="matlab">
disp(sprintf('Gravity estimated to %g',gopt));
</source>

<pre>
Gravity estimated to 9.82655

</pre>

PROPT Path Tracking Robot (Two-Phase)

2011-11-09T12:11:15Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

User's Guide for DIRCOL

2.7 Optimal path tracking for a simple robot. A robot with two rotational joints and simplified equations of motion has to move along a prescribed path with constant velocity.

==Problem Formulation==

Find u over t in [0; 2 ] to minimize

<math> J = \int_0^{2} (\sum_{i=1}^{2} (w_i*(q_i(t) - q_{i,ref})^2) + \sum_{i=1}^{2} (w_{2+i}*(\frac{dq}{dt}_i(t) - \frac{dq}{dt}_{i,ref})^2) \mathrm{d}t </math>

subject to:

<math> \frac{d^2q_1}{dt^2} = u_1 </math>

<math> \frac{d^2q_2}{dt^2} = u_2 </math>

A transformation gives:

<math> \frac{dx_1}{dt} = x_3 </math>

<math> \frac{dx_2}{dt} = x_4 </math>

<math> \frac{dx_3}{dt} = u_1 </math>

<math> \frac{dx_4}{dt} = u_2 </math>

<math> x_{1:4}(0) = [0 \ 0 \ 0.5 \ 0] </math>

<math> x_{1:4}(2) = [0.5 \ 0.5 \ 0 \ 0.5] </math>

<math> w_{1:4} = [100 \ 100 \ 500 \ 500] </math>

<math> x_1{1,ref} = \frac{t}{2} \ (0<t<1), \ \frac{1}{2} \ (1<t<2) </math>

<math> x_2{1,ref} = 0 \ (0<t<1), \ \frac{t-1}{2} \ (1<t<2) </math>

<math> x_3{1,ref} = \frac{1}{2} \ (0<t<1), \ 0 \ (1<t<2) </math>

<math> x_4{1,ref} = 0 \ (0<t<1), \ \frac{1}{2} \ (1<t<2) </math>

<math> |u| < 10 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
t_F = 2;
toms t1
p1 = tomPhase('p1', t1, 0, t_F/2, 25);
toms t2
p2 = tomPhase('p2', t2, t_F/2, t_F/2, 25);
setPhase(p1);

tomStates x1p1 x2p1 x3p1 x4p1
tomControls u1p1 u2p1

% Box constraints
cbox1 = {-10 <= collocate(u1p1) <= 10
-10 <= collocate(u2p1) <= 10};

% Boundary constraints
w1 = 100; w2 = 100;
w3 = 500; w4 = 500;
err1p1 = w1*(x1p1-t1/2).^2;
err2p1 = w2*(x2p1).^2;
err3p1 = w3*(x3p1-1/2).^2;
err4p1 = w4*(x4p1).^2;

cbnd1 = initial({x1p1 == 0; x2p1 == 0
x3p1 == 0.5; x4p1 == 0});

% ODEs and path constraints
ceq1 = collocate({dot(x1p1) == x3p1
dot(x2p1) == x4p1; dot(x3p1) == u1p1
dot(x4p1) == u2p1});

% Objective
objective1 = integrate(err1p1+err2p1+err3p1+err4p1);

% Phase 2
setPhase(p2);
tomStates x1p2 x2p2 x3p2 x4p2
tomControls u1p2 u2p2

% Box constraints
cbox2 = {-10 <= collocate(u1p2) <= 10
-10 <= collocate(u2p2) <= 10};

% Boundary constraints
err1p2 = w1*(x1p2-1/2).^2;
err2p2 = w2*(x2p2-(t2-1)/2).^2;
err3p2 = w3*(x3p2).^2;
err4p2 = w4*(x4p2-1/2).^2;

cbnd2 = final({x1p2 == 0.5
x2p2 == 0.5
x3p2 == 0
x4p2 == 0.5});

% ODEs and path constraints
ceq2 = collocate({
dot(x1p2) == x3p2
dot(x2p2) == x4p2
dot(x3p2) == u1p2
dot(x4p2) == u2p2});

% Objective
objective2 = integrate(err1p2+err2p2+err3p2+err4p2);

% Objective
objective = objective1 + objective2;

% Link phase
link = {final(p1,x1p1) == initial(p2,x1p2)
final(p1,x2p1) == initial(p2,x2p2)
final(p1,x3p1) == initial(p2,x3p2)
final(p1,x4p1) == initial(p2,x4p2)};
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Path Tracking Robot (Two-Phase)';
options.solver = 'sqopt7';
constr = {cbox1, cbnd1, ceq1, cbox2, cbnd2, ceq2, link};
solution = ezsolve(objective, constr, [], options);
t = subs(collocate(p1,t1),solution);
t = [t;subs(collocate(p2,t2),solution)];
x1 = subs(collocate(p1,x1p1),solution);
x1 = [x1;subs(collocate(p2,x1p2),solution)];
x2 = subs(collocate(p1,x2p1),solution);
x2 = [x2;subs(collocate(p2,x2p2),solution)];
x3 = subs(collocate(p1,x3p1),solution);
x3 = [x3;subs(collocate(p2,x3p2),solution)];
x4 = subs(collocate(p1,x4p1),solution);
x4 = [x4;subs(collocate(p2,x4p2),solution)];
u1 = subs(collocate(p1,u1p1),solution);
u1 = [u1;subs(collocate(p2,u1p2),solution)];
u2 = subs(collocate(p1,u2p1),solution);
u2 = [u2;subs(collocate(p2,u2p2),solution)];
</source>

<pre>
Problem type appears to be: qp
Time for symbolic processing: 0.25294 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: 1: Path Tracking Robot (Two-Phase) f_k 1.049952991376699200
sum(|constr|) 0.000000000261665098
f(x_k) + sum(|constr|) 1.049952991638364300
f(x_0) 0.000000000000000000

Solver: SQOPT. EXIT=0. INFORM=1.
SQOPT 7.2-5 QP solver
Optimality conditions satisfied

Iter 297
CPU time: 0.046800 sec. Elapsed time: 0.123000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1);
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Path Tracking Robot state variables');

subplot(2,1,2);
plot(t,u1,'*-',t,u2,'*-');
legend('u1','u2');
title('Path Tracking Robot control variables');
</source>

[[File:pathTrackingRobot_twophase_01.png]]

PROPT Path Tracking Robot

2011-11-09T12:11:10Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

User's Guide for DIRCOL

2.7 Optimal path tracking for a simple robot. A robot with two rotational joints and simplified equations of motion has to move along a prescribed path with constant velocity.

==Problem Formulation==

Find u over t in [0; 2 ] to minimize

<math> J = \int_0^{2} (\sum_{i=1}^{2} (w_i*(q_i(t) - q_{i,ref})^2) + \sum_{i=1}^{2} (w_{2+i}*(\frac{dq}{dt}_i(t) - \frac{dq}{dt}_{i,ref})^2) \mathrm{d}t </math>

subject to:

<math> \frac{d^2q_1}{dt^2} = u_1 </math>

<math> \frac{d^2q_2}{dt^2} = u_2 </math>

A transformation gives:

<math> \frac{dx_1}{dt} = x_3 </math>

<math> \frac{dx_2}{dt} = x_4 </math>

<math> \frac{dx_3}{dt} = u_1 </math>

<math> \frac{dx_4}{dt} = u_2 </math>

<math> x_{1:4}(0) = [0 \ 0 \ 0.5 \ 0] </math>

<math> x_{1:4}(2) = [0.5 \ 0.5 \ 0 \ 0.5] </math>

<math> w_{1:4} = [100 \ 100 \ 500 \ 500] </math>

<math> x_1{1,ref} = \frac{t}{2} \ (0<t<1), \ \frac{1}{2} \ (1<t<2) </math>

<math> x_2{1,ref} = 0 \ (0<t<1), \ \frac{t-1}{2} \ (1<t<2) </math>

<math> x_3{1,ref} = \frac{1}{2} \ (0<t<1), \ 0 \ (1<t<2) </math>

<math> x_4{1,ref} = 0 \ (0<t<1), \ \frac{1}{2} \ (1<t<2) </math>

<math> |u| < 10 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 2, 100, [], 'fem1s'); % Use splines with FEM constraints
%p = tomPhase('p', t, 0, 2, 100, [], 'fem1'); % Use linear finite elements
%p = tomPhase('p', t, 0, 2, 100); % Use Gauss point collocation
setPhase(p);

tomStates x1 x2 x3 x4
tomControls u1 u2

% Box constraints
cbox = {
-10 <= collocate(u1) <= 10
-10 <= collocate(u2) <= 10};

% Boundary constraints
cbnd = {initial({x1 == 0; x2 == 0
x3 == 0.5; x4 == 0})
final({x1 == 0.5; x2 == 0.5
x3 == 0; x4 == 0.5})};

% ODEs and path constraints
w1 = 100; w2 = 100;
w3 = 500; w4 = 500;

err1 = w1*(x1-t/2.*(t<1)-1/2*(t>=1)).^2;
err2 = w2*(x2-(t-1)/2.*(t>=1)).^2;
err3 = w3*(x3-1/2*(t<1)).^2;
err4 = w4*(x4-1/2*(t>=1)).^2;

toterr = integrate(err1+err2+err3+err4);

ceq = collocate({
dot(x1) == x3
dot(x2) == x4
dot(x3) == u1
dot(x4) == u2});

% Objective
objective = toterr;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Path Tracking Robot';
solution = ezsolve(objective, {cbox, cbnd, ceq}, [], options);
t = subs(icollocate(t),solution);
x1 = subs(icollocate(x1),solution);
x2 = subs(icollocate(x2),solution);
x3 = subs(icollocate(x3),solution);
x4 = subs(icollocate(x4),solution);
u1 = subs(icollocate(u1),solution);
u2 = subs(icollocate(u2),solution);
</source>

<pre>
Problem type appears to be: qp
Time for symbolic processing: 0.16121 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: 1: Path Tracking Robot f_k 1.031157512955303400
sum(|constr|) 0.000000055682846425
f(x_k) + sum(|constr|) 1.031157568638149800
f(x_0) 0.000000000000000000

Solver: CPLEX. EXIT=0. INFORM=1.
CPLEX Barrier QP solver
Optimal solution found

FuncEv 10 GradEv 10 ConstrEv 10 Iter 10
CPU time: 0.249602 sec. Elapsed time: 0.061000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1);
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Path Tracking Robot state variables');

subplot(2,1,2);
plot(t,u1,'*-',t,u2,'*-');
legend('u1','u2');
title('Path Tracking Robot control variables');
</source>

[[File:pathTrackingRobot_01.png]]

PROPT Park-Ramirez bioreactor

2011-11-09T12:11:06Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Dynamic optimization of chemical and biochemical processes using restricted second-order information 2001, Eva Balsa-Canto, Julio R. Banga, Antonio A. Alonso Vassilios S. Vassiliadis

Case Study I: Park-Ramirez bioreactor

==Problem description==

This case study deals with the optimal production of secreted protein in a fed-batch reactor. It was originally formulated by Park and Ramirez (Park & Ramirez, 1988) and it has also been considered by other researchers (Vassiliadis, 1993; Yang & Wu, 1994; Banga, Irizarry & Seider, 1995; Luus, 1995; Tholudur & Ramirez, 1997). The objective is to maximize the secreted heterologous protein by a yeast strain in a fed-batch culture. The dynamic model accounts for host-cell growth, gene expression, and the secretion of expressed polypeptides. The mathematical statement is as follows:

Find u(t) over t in [t0,t_f] to maximize:

<math> J = x_1(t_f)*x_5(t_f) </math>

subject to:

<math> \frac{dx_1}{dt} = g_1*(x_2-x_1)-\frac{u}{x_5}*x_1 ,</math>

<math> \frac{dx_2}{dt} = g_2*x_3-\frac{u}{x_5}*x_2 ,</math>

<math> \frac{dx_3}{dt} = g_3*x_3-\frac{u}{x_5}*x_3 ,</math>

<math> \frac{dx_4}{dt} = -7.3*g_3*x_3+\frac{u}{x_5}*(20-x_4) ,</math>

<math> \frac{dx_5}{dt} = u ,</math>

with:

<math> g_1 = 4.75*\frac{g_3}{(0.12+g_3)} ,</math>

<math> g_2 = \frac{x_4}{(0.1+x4)}*exp(-5*x_4) ,</math>

<math> g_3 = 21.87*\frac{x_4}{(x_4+0.4)*(x_4+62.5)} ,</math>

where x1 and x2 are, respectively, the concentration of the secreted protein and the total protein (l-1), x3 is the culture cell density (g l-1), x4 is the substrate concentration (g l-1), x5 is the holdup volume (l), u is the nutrient (glucose) feed rate (l h-1), and J is the mass of protein produced (in arbitrary units). The initial conditions are:

<math> x(t_0) = [0 \ 0 \ 1 \ 5 \ 1]' </math>

For final time t_f = 15 h, and the following constraints on the control variable:

<math> 0 <= u <= 2 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
</source>

==Solve the problem, using a successively larger number collocation points==

<source lang="matlab">
for n=[20 40 80 120]
</source>

<source lang="matlab">
p = tomPhase('p', t, 0, 15, n);
setPhase(p);

tomStates z1 z2 z3 z4 z5
tomControls u

if n>20
% Interpolate an initial guess for the n collocation points
x0 = {icollocate({z1 == z1opt; z2 == z2opt
z3 == z3opt; z4 == z4opt; z5 == z5opt})
collocate(u == uopt)};
else
x0 = {};
end

% Box constraints
cbox = {icollocate({z1 <= 3
z2 <= 3; 0 <= z3 <= 4
0 <= z4 <= 10; 0.5 <= z5 <= 25})
0 <= collocate(u) <= 2.5};

% Boundary constraints
cbnd = initial({z1 == 0; z2 == 0; z3 == 1
z4 == 5; z5 == 1});

% Various constants and expressions
g3 = 21.87*z4./(z4+.4)./(z4+62.5);
g1 = 4.75*g3./(0.12+g3);
g2 = z4./(0.1+z4).*exp(-5*z4);

% ODEs and path constraints
ceq = collocate({
dot(z1) == g1.*(z2-z1)-u./z5.*z1
dot(z2) == g2.*z3-u./z5.*z2
dot(z3) == g3.*z3-u./z5.*z3
dot(z4) == -7.3*g3.*z3+u./z5.*(20-z4)
dot(z5) == u});

% Secreted protein must be less than total protein
% proteinlimit = {z1 <= z2};

% Objective
if n == 120
objective = -final(z1)*final(z5)+var(diff(collocate(u)));
else
objective = -final(z1)*final(z5);
end
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Park Bio Reactor';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

% Optimal z and u, to use as starting guess in the
% next iteration
z1opt = subs(z1, solution);
z2opt = subs(z2, solution);
z3opt = subs(z3, solution);
z4opt = subs(z4, solution);
z5opt = subs(z5, solution);
uopt = subs(u, solution);
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.48856 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Park Bio Reactor f_k -31.891604942090204000
sum(|constr|) 0.000000000222835644
f(x_k) + sum(|constr|) -31.891604941867367000
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 102 ConJacEv 102 Iter 61 MinorIter 911
CPU time: 0.171601 sec. Elapsed time: 0.172000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.48207 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Park Bio Reactor f_k -29.523702508277267000
sum(|constr|) 0.000000000531678452
f(x_k) + sum(|constr|) -29.523702507745590000
f(x_0) -31.891604942089948000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 75 ConJacEv 75 Iter 58 MinorIter 588
CPU time: 0.436803 sec. Elapsed time: 0.434000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.4837 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Park Bio Reactor f_k -31.866762161406442000
sum(|constr|) 0.000000000726937524
f(x_k) + sum(|constr|) -31.866762160679503000
f(x_0) -29.523702508277417000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 137 ConJacEv 137 Iter 120 MinorIter 784
CPU time: 4.586429 sec. Elapsed time: 4.584000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.51719 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Park Bio Reactor f_k -32.616144466384739000
sum(|constr|) 0.000000000312031161
f(x_k) + sum(|constr|) -32.616144466072711000
f(x_0) -31.477810607457531000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 134 ConJacEv 134 Iter 116 MinorIter 930
CPU time: 46.815900 sec. Elapsed time: 16.326000 sec.

</pre>

<source lang="matlab">
end

t = subs(collocate(t),solution);
z1 = collocate(z1opt);
z2 = collocate(z2opt);
z3 = collocate(z3opt);
z4 = collocate(z4opt);
z5 = collocate(z5opt);
u = collocate(uopt);
</source>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,z1,'*-',t,z2,'*-',t,z3,'*-',t,z4,'*-',t,z5,'*-');
legend('z1','z2','z3','z4','z5');
title('Park Bio Reactor state variables');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Park Bio Reactor control');
</source>

[[File:parkBioReactor_01.png]]

PROPT Parameter Estimation Problem

2011-11-09T12:10:52Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Example 5: DYNOPT User's Guide version 4.1.0

M. Cizniar, M. Fikar, M. A. Latifi, MATLAB Dynamic Optimisation Code DYNOPT. User's Guide, Technical Report, KIRP FCHPT STU Bratislava, Slovak Republic, 2006.

==Problem description==

Find p1 and p2 over t in [0; 6 ] to minimize

<math> J = \sum_{i=1,2,3,5} ( x_1(t_i) - x_1^{m}(t_i) )^2 </math>

subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = 1-2*x_2-x_1 </math>

where

<math> x_0 = [p_1 \ p_2] </math>

<math> t_i = [1 \ 2 \ 3 \ 5] </math>

<math> x_1^{m}(t_i) = [0.264 \ 0.594 \ 0.801 \ 0.959] </math>

<math> -1.5 <= p_{1:2} <= 1.5 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t p1 p2
x1meas = [0.264;0.594;0.801;0.959];
tmeas = [1;2;3;5];

% Box constraints
cbox = {-1.5 <= p1 <= 1.5
-1.5 <= p2 <= 1.5};
</source>

==Solve the problem, using a successively larger number collocation points==

<source lang="matlab">
for n=[10 40]
</source>

<source lang="matlab">
p = tomPhase('p', t, 0, 6, n);
setPhase(p);
tomStates x1 x2

% Initial guess
if n == 10
x0 = {p1 == 0; p2 == 0};
else
x0 = {p1 == p1opt; p2 == p2opt
icollocate({x1 == x1opt; x2 == x2opt})};
end

% Boundary constraints
cbnd = initial({x1 == p1; x2 == p2});

% ODEs and path constraints
x1err = sum((atPoints(tmeas,x1) - x1meas).^2);
ceq = collocate({dot(x1) == x2; dot(x2) == 1-2*x2-x1});

% Objective
objective = x1err;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Parameter Estimation';
options.solver = 'snopt';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

% Optimal x, p for starting point
x1opt = subs(x1, solution);
x2opt = subs(x2, solution);
p1opt = subs(p1, solution);
p2opt = subs(p2, solution);
</source>

<pre>
Problem type appears to be: qp
Time for symbolic processing: 0.049765 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: 1: Parameter Estimation f_k 0.000000352979271145
sum(|constr|) 0.000000000000054134
f(x_k) + sum(|constr|) 0.000000352979325279
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 Iter 5 MinorIter 17
Elapsed time: 0.002000 sec.

</pre>

<pre>
Problem type appears to be: qp
Time for symbolic processing: 0.051384 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: 1: Parameter Estimation f_k 0.000000355693079657
sum(|constr|) 0.000000000000071929
f(x_k) + sum(|constr|) 0.000000355693151586
f(x_0) -1.983813647020728800

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 MinorIter 40
Elapsed time: 0.003000 sec.

</pre>

<source lang="matlab">
end

t = subs(collocate(t),solution);
x1 = collocate(x1opt);
</source>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,x1,'*-',tmeas,x1meas,'ro');
legend('x1','Meas');
title('Parameter Estimation state variable');
</source>

[[File:parameterEstimation_01.png]]

PROPT Parallel Reactions in Tubular Reactor

2011-11-09T12:10:45Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Problem 4: DYNOPT User's Guide version 4.1.0

Batch reactor with reactions: A -> B and A -> C.

M. Cizniar, M. Fikar, M. A. Latifi, MATLAB Dynamic Optimisation Code DYNOPT. User's Guide, Technical Report, KIRP FCHPT STU Bratislava, Slovak Republic, 2006.

==Problem description==

Find T over t in [0; 1 ] to maximize

<math> J = x_2(t_f) </math>

subject to:

<math> \frac{dx_1}{dt} = -(u+0.5*u^2)*x_1 </math>

<math> \frac{dx_2}{dt} = u*x_1 </math>

where

<math> x(0) = [1 \ 0] </math>

<math> 0 <= u <= 5 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
p = tomPhase('p', t, 0, 1, 30);
setPhase(p);

tomStates x1 x2
tomControls u

% Initial guess
x0 = {icollocate({x1 == 1; x2 == 0})
collocate(u == 5*t)};

% Box constraints
cbox = {0 <= collocate(u) <= 5};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 0});

% ODEs and path constraints
ceq = collocate({
dot(x1) == -(u+0.5*u.^2).*x1
dot(x2) == u.*x1});

% Objective
objective = -final(x2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Parallel Reactions';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
t = subs(collocate(t),solution);
x1 = subs(collocate(x1),solution);
x2 = subs(collocate(x2),solution);
u = subs(collocate(u),solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.087774 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Parallel Reactions f_k -0.573540787113988930
sum(|constr|) 0.000000228705850230
f(x_k) + sum(|constr|) -0.573540558408138670
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 35 ConJacEv 35 Iter 34 MinorIter 118
CPU time: 0.046800 sec. Elapsed time: 0.042000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-');
legend('x1','x2');
title('Parallel Reactions state variables');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Parallel Reactions control');
</source>

[[File:parallelReactions_01.png]]

PROPT Orbit Raising Minimum Time

2011-11-09T12:10:42Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

==Problem description==

Minimize:

<math> J = t_f </math>

subject to the dynamic constraints

<math> \frac{d_r}{dt} = u </math>

<math> \frac{du}{dt} = \frac{v^2}{r}-\frac{mmu}{r^2}+T*\frac{w1}{m} </math>

<math> \frac{dv}{dt} = -u*\frac{v}{r}+T*\frac{w2}{m} </math>

<math> \frac{dm}{dt} = -\frac{T}{g0*ISP} </math>

the boundary conditions

<math> r(t_0) = 1 </math>

<math> u(t_0) = 1 </math>

<math> u(t_f) = 0 </math>

<math> v(t_0) = (\frac{mmu}{r(t_0)})^{0.5} </math>

<math> v(t_f) = (\frac{mmu}{r(t_f)})^{0.5} </math>

<math> m(t_0) = 1 </math>

where w1 = sin(phi) and w2 = cos(phi)

At t_f, r and m are free.

<source lang="matlab">
% Copyright (c) 2007-2009 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
mmu = 1;
t_f = 3.32; m_0 = 1; r_0 = 1; u_0 = 0;
u_f = 0; v_0 = sqrt(mmu/r_0); rmin = 0.9; rmax = 5;
umin = -5; umax = 5; vmin = -5; vmax = 5;
mmax = m_0; mmin = 0.1; tf_min = 0.5; tf_max = 10;
r_f = 1.5;

T = 0.1405;
Ve = 1.8758;

toms t t_f
p1 = tomPhase('p1', t, 0, t_f, 50);
setPhase(p1);

tomStates r u v m

% The problem becomes less nonlinear if w1 and w2 are control variables
% (with the constraints w1^2+w2^2==1) than if phi is the control varialbe
% (with w1 and w2 being nonlinear functions of phi).
tomControls w1 w2
phi = atan2(w1,w2);

% Initial guess
x0 = {t_f == 3.32
icollocate({
r == r_0+(r_f-r_0)*t/t_f
u == 0.1
v == v_0
m == m_0-(T/Ve)*t})
collocate({
w1 == -0.7*sign(t-t_f/2)
w2 == 0.4
})};

% Boundary constraints
cbnd = {initial({
r == r_0
u == u_0
v == v_0
m == m_0
})
final({
r == r_f
u == u_f
v == sqrt(mmu/r)})};

% Box constraints
cbox = {0.5 <= t_f <= 10
rmin <= icollocate(r) <= rmax
umin <= icollocate(u) <= umax
vmin <= icollocate(v) <= vmax
};

% ODEs and path constraints
ceq = collocate({
dot(r) == u
dot(u) == v^2/r-mmu/r^2+T*w1/m
dot(v) == -u*v/r+T*w2/m
dot(m) == -T/Ve
w1^2+w2^2 == 1
});

% Objective
objective = t_f;
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Orbit Raising Problem Min Time';
options.scale = 'manual'; % Auto-scaling is not really needed as all variables are already reasonably scaled.
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.26997 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Orbit Raising Problem Min Time f_k 3.248079535631201800
sum(|constr|) 0.000032253387193478
f(x_k) + sum(|constr|) 3.248111789018395300
f(x_0) 3.319999999999999800

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 68 ConJacEv 68 Iter 35 MinorIter 283
CPU time: 1.606810 sec. Elapsed time: 0.414000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
ezplot([r u v m]);
legend('r','u','v','m');
title('Orbit Raising Problem Min Time state variables');

subplot(2,1,2)
ezplot([w1 w2 phi])
legend('w_1', 'w_2', '\phi');
title('Orbit Raising Problem Min Time control');
</source>

[[File:orbitRaisingMinTime_01.png]]

PROPT Orbit Raising Maximum Radius

2011-11-09T12:10:40Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

==Problem description==

Maximize:

<math> J = r(t_f) </math>

subject to the dynamic constraints

<math> \frac{d_r}{dt} = u </math>

<math> \frac{du}{dt} = \frac{v^2}{r}-\frac{mmu}{r^2}+T*\frac{w1}{m} </math>

<math> \frac{dv}{dt} = -u*\frac{v}{r}+T*\frac{w2}{m} </math>

<math> \frac{dm}{dt} = -\frac{T}{g0*ISP} </math>

the boundary conditions

<math> r(t_0) = 1 </math>

<math> u(t_0) = 1 </math>

<math> u(t_f) = 0 </math>

<math> v(t_0) = (\frac{mmu}{r(t_0)})^{0.5} </math>

<math> v(t_f) = (\frac{mmu}{r(t_f)})^{0.5} </math>

<math> m(t_0) = 1 </math>

and the path constraint

<math> w_1^2+w_2^2 = 1 </math>

The control pitch angel is not being used in this formulation. Instead two control variables (w1,w2) are used to for the thrust direction. A path constraint ensures that (w1,w2) is a unit vector.

<source lang="matlab">
% Copyright (c) 2007-2009 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
t_F = 3.32;
mmu = 1;
m_0 = 1; r_0 = 1; u_0 = 0;
u_f = 0; v_0 = sqrt(mmu/r_0); rmin = 0.9;
rmax = 5; umin = -5; umax = 5;
vmin = -5; vmax = 5; mmax = m_0;
mmin = 0.1;

T = 0.1405;
Ve = 1.8758;

toms t
p1 = tomPhase('p1', t, 0, t_F, 40);
setPhase(p1);

tomStates r u v m
tomControls w1 w2

% Initial guess
x0 = {icollocate({
r == r_0
u == u_0 + (u_f-u_0)*t/t_F
v == v_0
m == m_0})
collocate({w1 == 0; w2 == 1})};

% Boundary constraints
cbnd = {initial({
r == r_0
u == u_0
v == v_0
m == m_0
})
final({u == u_f
v - sqrt(mmu/r) == 0})};

% Box constraints
cbox = {
rmin <= icollocate(r) <= rmax
umin <= icollocate(u) <= umax
vmin <= icollocate(v) <= vmax
mmin <= icollocate(m) <= mmax
-1 <= collocate(w1) <= 1
-1 <= collocate(w2) <= 1};

% ODEs and path constraints
ceq = collocate({
dot(r) == u
dot(u) == v^2/r-mmu/r^2+T*w1/m
dot(v) == -u*v/r+T*w2/m
dot(m) == -T/Ve
w1^2+w2^2 == 1});

% Objective
objective = -final(r);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Orbit Raising Problem Max Radius';
options.solver = 'snopt';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.22212 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Orbit Raising Problem Max Radius f_k -1.518744202740535800
sum(|constr|) 0.000017265830060524
f(x_k) + sum(|constr|) -1.518726936910475200
f(x_0) -0.999999999999991120

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 59 ConJacEv 59 Iter 36 MinorIter 346
CPU time: 0.202801 sec. Elapsed time: 0.226000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
ezplot([r u v m]);
legend('r','u','v','m');
title('Orbit Raising Problem Max Radius state variables');

subplot(2,1,2)
ezplot([w1 w2 atan2(w1,w2)])
legend('w_1', 'w_2', '\phi');
title('Orbit Raising Problem Max Radius control');
</source>

[[File:orbitRaisingMaxRadius_01.png]]

PROPT One Dimensional Rocket Ascent

2011-11-09T12:10:31Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

User's Guide for DIRCOL

Problem 2.3 One-dimensional ascent of a rocket

==Problem Formulation==

Find tCut over t in [0; t_F ] to minimize

<math> J = tCut </math>

subject to:

<math> \frac{dx_1}{dt} = x2 </math>

<math> \frac{dx_2}{dt} = a-g \ (0<t<tCut)</math>

<math> \frac{dx_2}{dt} = -g \ (tCut<t<t_F)</math>

<math> [x_1 \ x_2] = 0 </math>

<math> g = 1 </math>

<math> a = 2 </math>

<math> x_1(t_f) = 100 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
toms tCut tp2
p1 = tomPhase('p1', t, 0, tCut, 20);
p2 = tomPhase('p2', t, tCut, tp2, 20);

t_f = tCut+tp2;

x1p1 = tomState(p1,'x1p1');
x2p1 = tomState(p1,'x2p1');
x1p2 = tomState(p2,'x1p2');
x2p2 = tomState(p2,'x2p2');

% Initial guess
x0 = {tCut==10
t_f==15
icollocate(p1,{x1p1 == 50*tCut/10;x2p1 == 0;})
icollocate(p2,{x1p2 == 50+50*t/100;x2p2 == 0;})};

% Box constraints
cbox = {
1 <= tCut <= t_f-0.00001
t_f <= 100
0 <= icollocate(p1,x1p1)
0 <= icollocate(p1,x2p1)
0 <= icollocate(p2,x1p2)
0 <= icollocate(p2,x2p2)};

% Boundary constraints
cbnd = {initial(p1,{x1p1 == 0;x2p1 == 0;})
final(p2,x1p2 == 100)};

% ODEs and path constraints
a = 2; g = 1;
ceq = {collocate(p1,{
dot(p1,x1p1) == x2p1
dot(p1,x2p1) == a-g})
collocate(p2,{
dot(p2,x1p2) == x2p2
dot(p2,x2p2) == -g})};

% Objective
objective = tCut;

% Link phase
link = {final(p1,x1p1) == initial(p2,x1p2)
final(p1,x2p1) == initial(p2,x2p2)};
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'One Dim Rocket';
constr = {cbox, cbnd, ceq, link};
solution = ezsolve(objective, constr, x0, options);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.11449 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: One Dim Rocket f_k 9.999998166162907200
sum(|constr|) 0.000733735151552596
f(x_k) + sum(|constr|) 10.000731901314460000
f(x_0) 10.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 10 ConJacEv 10 Iter 8 MinorIter 91
CPU time: 0.015600 sec. Elapsed time: 0.015000 sec.

</pre>

==Plot result==

<source lang="matlab">
t = [subs(collocate(p1,t),solution);subs(collocate(p2,t),solution)];
x1 = subs(collocate(p1,x1p1),solution);
x1 = [x1;subs(collocate(p2,x1p2),solution)];
x2 = subs(collocate(p1,x2p1),solution);
x2 = [x2;subs(collocate(p2,x2p2),solution)];

figure(1)
plot(t,x1,'*-',t,x2,'*-');
legend('x1','x2');
title('One Dim Rocket state variables');
</source>

[[File:onedimRocket_01.png]]

PROPT Oil Shale Pyrolysis

2011-11-09T12:10:29Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Dynamic Optimization of Batch Reactors Using Adaptive Stochastic Algorithms 1997, Eugenio F. Carrasco, Julio R. Banga

Case Study II: Oil Shale Pyrolysis

==Problem description==

A very challenging optimal control problem is the oil shale pyrolysis case study, as considered by Luus (1994). The system consists of a series of five chemical reactions:

A1 -> A2

A2 -> A3

A1+A2 -> A2+A2

A1+A2 -> A3+A2

A1+A2 -> A4+A2

This system is described by the differential equations

<math> \frac{dx_1}{dt} = -k_1*x_1-(k_3+k_4+k_5)*x_1*x_2 </math>

<math> \frac{dx_2}{dt} = k_1*x_1-k_2*x_2+k_3*x_1*x_2 </math>

<math> \frac{dx_3}{dt} = k_2*x_2+k_4*x_1*x_2 </math>

<math> \frac{dx_4}{dt} = k_5*x_1*x_2 </math>

where the state variables are the concentrations of species, Ai, i = 1, ..., 4. The initial condition is

<math> x(t_0) = [1 \ 0 \ 0 \ 0]' </math>

The rate expressions are given by:

<math> k_i = k_{i0}*exp(-\frac{Ei}{R*T}), i=1,2,3,4,5 </math>

where the values of ki0 and Ei are given by Luus (1994). The optimal control problem is to find the residence time t_f and the temperature profile T(t) in the time interval 0 <= t <= t_f so that the production of pyrolytic bitumen, given by x2, is maximized. Therefore, the performance index is

<math> J = x_2(t_f) </math>

The constraints on the control variable (temperature) are:

<math> 698.15 <= T <= 748.15 </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
toms t_f

ai = [8.86; 24.25; 23.67; 18.75; 20.70];
bi = [20300; 37400; 33800; 28200; 31000]/1.9872;

for n=[4 10 20 30 35]
</source>

<source lang="matlab">
p = tomPhase('p', t, 0, t_f, n);
setPhase(p);

tomStates x1 x2 x3 x4
tomControls T

% Initial guess
if n == 4
x0 = {t_f == 9.3
collocate(T == 725)};
else
x0 = {t_f == tfopt
icollocate({
x1 == x1opt; x2 == x2opt
x3 == x3opt; x4 == x4opt
})
collocate(T == Topt)};
end

% Box constraints
cbox = {9.1 <= t_f <= 12
icollocate({0 <= x1 <= 1; 0 <= x2 <= 1
0 <= x3 <= 1; 0 <= x4 <= 1})
698.15 <= collocate(T) <= 748.15};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 0; x3 == 0; x4 == 0});

% ODEs and path constraints

ki1 = exp(ai(1)-bi(1)./T);
ki2 = exp(ai(2)-bi(2)./T);
ki3 = exp(ai(3)-bi(3)./T);
ki4 = exp(ai(4)-bi(4)./T);
ki5 = exp(ai(5)-bi(5)./T);

ceq = collocate({
dot(x1) == -ki1.*x1-(ki3+ki4+ki5).*x1.*x2
dot(x2) == ki1.*x1-ki2.*x2+ki3.*x1.*x2
dot(x3) == ki2.*x2+ki4.*x1.*x2
dot(x4) == ki5.*x1.*x2});

% Objective
objective = -final(x2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Oil Pyrolysis';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

x1opt = subs(x1, solution);
x2opt = subs(x2, solution);
x3opt = subs(x3, solution);
x4opt = subs(x4, solution);
Topt = subs(T, solution);
tfopt = subs(final(t), solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.42383 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Oil Pyrolysis f_k -0.357327805323273240
sum(|constr|) 0.000000000957541036
f(x_k) + sum(|constr|) -0.357327804365732190
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 93 ConJacEv 93 Iter 50 MinorIter 197
CPU time: 0.062400 sec. Elapsed time: 0.064000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.42239 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Oil Pyrolysis f_k -0.354368283541904860
sum(|constr|) 0.000000002719365042
f(x_k) + sum(|constr|) -0.354368280822539790
f(x_0) -0.357327805323273070

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 208 ConJacEv 208 Iter 112 MinorIter 305
CPU time: 0.156001 sec. Elapsed time: 0.153000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.43305 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Oil Pyrolysis f_k -0.351747594492437030
sum(|constr|) 0.000000280199856233
f(x_k) + sum(|constr|) -0.351747314292580830
f(x_0) -0.354368283541905970

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 124 ConJacEv 124 Iter 71 MinorIter 281
CPU time: 0.140401 sec. Elapsed time: 0.142000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.42125 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Oil Pyrolysis f_k -0.352833701465704920
sum(|constr|) 0.000000018165177960
f(x_k) + sum(|constr|) -0.352833683300526950
f(x_0) -0.351747594492436640

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 409 ConJacEv 409 Iter 194 MinorIter 697
CPU time: 0.639604 sec. Elapsed time: 0.645000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.42103 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Oil Pyrolysis f_k -0.352618526247765740
sum(|constr|) 0.000016167955157926
f(x_k) + sum(|constr|) -0.352602358292607830
f(x_0) -0.352833701465704590

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 64 ConJacEv 64 Iter 46 MinorIter 364
CPU time: 0.187201 sec. Elapsed time: 0.180000 sec.

</pre>

<source lang="matlab">
end

t = subs(collocate(t),solution);
x1 = subs(collocate(x1opt),solution);
x2 = subs(collocate(x2opt),solution);
x3 = subs(collocate(x3opt),solution);
x4 = subs(collocate(x4opt),solution);
T = subs(collocate(Topt),solution);
</source>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Oil Pyrolysis state variables');

subplot(2,1,2)
plot(t,T,'+-');
legend('T');
title('Oil Pyrolysis control');
</source>

[[File:oilPyrolysis_01.png]]

PROPT Obstacle Avoidance

2011-11-09T12:10:27Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

OPTRAGEN 1.0 A MATLAB Toolbox for Optimal Trajectory Generation, Raktim Bhattacharya, Texas A&M University (Note: There is typographical error in the OPTRAGEN documentation. The objective involves second derivatives of x and y.)

A robot with obstacles in 2D space. Travel from point A to B using minimum energy.

==Problem Formulation==

Find theta(t) and V over t in [0; 1 ] to minimize

<math> \int_0^{1} ( (\frac{d^2x}{dt^2})^2 + (\frac{d^2y}{dt^2})^2 ) \mathrm{d}t </math>

subject to:

<math> \frac{dx}{dt} = V*cos(theta) </math>

<math> \frac{dy}{dt} = V*sin(theta) </math>

<math> [x_0 \ y_0] = [0 \ 0] </math>

<math> [x_1 \ y_1] = [1.2 \ 1.6] </math>

<math> (x-0.4)^2 + (y-0.5)^2 >= 0.1 </math>

<math> (x-0.8)^2 + (y-1.5)^2 >= 0.1 </math>

Where V is a constant scalar speed.

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.

clear pi % 3.14159... unless someone has set it to something else!

toms t V

t_f = 1; % final time

% Starting guess
speed = 5;
xopt = 1.2*t;
yopt = 1.6*t;
vxopt = 1.2;
vyopt = 1.6;
thetaopt = pi/4;
</source>

==Solve the problem, using a successively larger number collocation points==

<source lang="matlab">
for n=[4 15 30]
</source>

<source lang="matlab">
% Create a new phase and states, using n collocation points
p = tomPhase('p', t, 0, t_f, n);
setPhase(p);
tomStates x y vx vy
tomControls theta

% Interpolate an initial guess for the n collocation points
x0 = {V == speed
icollocate({x == xopt; y == yopt; vx == vxopt; vy == vyopt})
collocate(theta == thetaopt)};

% Box constraints
cbox = {0 <= V <= 100 };

% Boundary constraints
cbnd = {initial({x == 0; y == 0})
final({x == 1.2; y == 1.6})};

% ODEs and path constraints
ode = collocate({
dot(x) == vx == V*cos(theta)
dot(y) == vy == V*sin(theta)
});

% A 30th order polynomial is more than sufficient to give good
% accuracy. However, that means that mcollocate would only check
% about 60 points. In order to make sure we don't hit an obstacle,
% we check 300 evenly spaced points instead, using atPoints.
obstacles = atPoints(linspace(0,t_f,300), {
(x-0.4)^2 + (y-0.5)^2 >= 0.1
(x-0.8)^2 + (y-1.5)^2 >= 0.1});

% Objective: minimum energy.
objective = integrate(dot(vx)^2+dot(vy)^2);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Obstacle avoidance';
constr = {cbox, cbnd, ode, obstacles};
solution = ezsolve(objective, constr, x0, options);

% Optimal x, y, and speed, to use as starting guess in the next iteration
xopt = subs(x, solution);
yopt = subs(y, solution);
vxopt = subs(vx, solution);
vyopt = subs(vy, solution);
thetaopt = subs(theta, solution);
speed = subs(V,solution);
</source>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.16968 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Obstacle avoidance f_k 29.812856165009947000
sum(|constr|) 0.000000001309307815
f(x_k) + sum(|constr|) 29.812856166319254000
f(x_0) 0.000000000000062528

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 22 ConJacEv 22 Iter 20 MinorIter 2732
CPU time: 0.218401 sec. Elapsed time: 0.218000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.17036 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Obstacle avoidance f_k 22.128728366252528000
sum(|constr|) 0.000000000006755102
f(x_k) + sum(|constr|) 22.128728366259281000
f(x_0) 29.812856165012775000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 151 ConJacEv 151 Iter 136 MinorIter 486
CPU time: 0.920406 sec. Elapsed time: 0.915000 sec.

</pre>

<pre>
Problem type appears to be: qpcon
Time for symbolic processing: 0.16882 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Obstacle avoidance f_k 22.091923284169869000
sum(|constr|) 0.000000000011094717
f(x_k) + sum(|constr|) 22.091923284180965000
f(x_0) 22.128728366262933000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 309 ConJacEv 309 Iter 281 MinorIter 720
CPU time: 4.071626 sec. Elapsed time: 3.770000 sec.

</pre>

<source lang="matlab">
end
</source>

==Plot result==

<source lang="matlab">
figure(1)
th = linspace(0,2*pi,500);
x1 = sqrt(0.1)*cos(th)+0.4;
y1 = sqrt(0.1)*sin(th)+0.5;
x2 = sqrt(0.1)*cos(th)+0.8;
y2 = sqrt(0.1)*sin(th)+1.5;
ezplot(x,y);
hold on
plot(x1,y1,'r',x2,y2,'r');
hold off
xlabel('x');
ylabel('y');
title(sprintf('Obstacle avoidance state variables, Speed = %2.4g',speed));
axis image
</source>

[[File:obstacleAvoidance_01.png]]

PROPT Nonlinear CSTR

2011-11-09T12:10:25Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Dynamic optimization of chemical and biochemical processes using restricted second-order information 2001, Eva Balsa-Canto, Julio R. Banga, Antonio A. Alonso Vassilios S. Vassiliadis

Case Study III: Nonlinear CSTR

==Problem description==

The problem was first introduced by Jensen (1964) and consists of determining the four optimal controls of a chemical reactor in order to obtain maximum economic benefit. The system dynamics describe four simultaneous chemical reactions taking place in an isothermal continuous stirred tank reactor. The controls are the flow rates of three feed streams and an electrical energy input used to promote a photochemical reaction. Luus (1990) and Bojkov, Hansel, and Luus (1993) considered two sub-cases using three and four control variables respectively.

The problem is formulated as follows:Find u1(t), u2(t), u3(t) and u4(t) over t in [t0,t_f] to maximize:

<math> J = x_8(t_f) </math>

Subject to:

<math> \frac{dx_1}{dt} = u_4-q*x_1-17.6*x_1*x_2-23*x_1*x_6*u_3 </math>

<math> \frac{dx_2}{dt} = u_1-q*x_2-17.6*x_1*x_2-146*x_2*x_3 </math>

<math> \frac{dx_3}{dt} = u_2-q*x_3-73*x_2*x_3 </math>

<math> \frac{dx_4}{dt} = -q*x_4+35.2*x_1*x_2-51.3*x_4*x_5 </math>

<math> \frac{dx_5}{dt} = -q*x_5+219*x_2*x_3-51.3*x_4*x_5 </math>

<math> \frac{dx_6}{dt} = -q*x_6+102.6*x_4.*x_5-23*x_1*x_6*u_3 </math>

<math> \frac{dx_7}{dt} = -q*x_7+46*x_1*x_6*u_3 </math>

<math> \frac{dx_8}{dt} = 5.8*(q*x_1-u_4)-3.7*u_1-4.1*u_2+ </math>

<math> q*(23*x_4+11*x_5+28*x_6+35*x_7)-5*u_3^2-0.099 </math>

where:

<math> q = u_1+u_2+u_4; </math>

with the initial conditions:

<math> x(t_0) = [0.1883 \ 0.2507 \ 0.0467 \ 0.0899 \ 0.1804 \ 0.1394 \ 0.1046 \ 0.000]' </math>

And the following bounds on the control variables:

<math> 0 <= u1 <= 20 </math>

<math> 0 <= u2 <= 6 </math>

<math> 0 <= u3 <= 4 </math>

<math> 0 <= u4 <= 20 </math>

The final time is considered fixed as t_f = 0.2.

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t
</source>

==Solve the problem, using a successively larger number collocation points==

<source lang="matlab">
for n=[5 20 60]
</source>

<source lang="matlab">
p = tomPhase('p', t, 0, 0.2, n);
setPhase(p);

tomStates x1 x2 x3 x4 x5 x6 x7 x8
tomControls u1 u2 u3 u4

% Interpolate an initial guess for the n collocation points
if n == 5
x0 = {};
else
x0 = {icollocate({x1 == x1opt; x2 == x2opt
x3 == x3opt; x4 == x4opt; x5 == x5opt
x6 == x6opt; x7 == x7opt; x8 == x8opt})
collocate({u1 == u1opt; u2 == u2opt
u3 == u3opt; u4 == u4opt})};
end

% Box constraints
cbox = {icollocate({
0 <= x1; 0 <= x2; 0 <= x3
0 <= x4; 0 <= x5; 0 <= x6
0 <= x7; 0 <= x8})
collocate({
0 <= u1 <= 20; 0 <= u2 <= 6
0 <= u3 <= 4; 0 <= u4 <= 20})};

% Boundary constraints
cbnd = initial({x1 == 0.1883; x2 == 0.2507
x3 == 0.0467; x4 == 0.0899; x5 == 0.1804
x6 == 0.1394; x7 == 0.1064; x8 == 0});

% ODEs and path constraints
% 4.1*u2+(u1+u2.*u4) in another paper, -0.09 instead of -0.099
q = u1+u2+u4;
ceq = collocate({
dot(x1) == (u4-q.*x1-17.6*x1.*x2-23*x1.*x6.*u3)
dot(x2) == (u1-q.*x2-17.6*x1.*x2-146*x2.*x3)
dot(x3) == (u2-q.*x3-73*x2.*x3)
dot(x4) == (-q.*x4+35.2*x1.*x2-51.3*x4.*x5)
dot(x5) == (-q.*x5+219*x2.*x3-51.3*x4.*x5)
dot(x6) == (-q.*x6+102.6*x4.*x5-23*x1.*x6.*u3)
dot(x7) == (-q.*x7+46*x1.*x6.*u3)
dot(x8) == (5.8*(q.*x1-u4)-3.7*u1-4.1*u2+q.*...
(23*x4+11*x5+28*x6+35*x7)-5*u3.^2-0.099)});

% Objective
objective = -final(x8);
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Nonlinear CSTR';
solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

% Optimal x and u as starting point
x1opt = subs(x1, solution);
x2opt = subs(x2, solution);
x3opt = subs(x3, solution);
x4opt = subs(x4, solution);
x5opt = subs(x5, solution);
x6opt = subs(x6, solution);
x7opt = subs(x7, solution);
x8opt = subs(x8, solution);
u1opt = subs(u1, solution);
u2opt = subs(u2, solution);
u3opt = subs(u3, solution);
u4opt = subs(u4, solution);
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.53917 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Nonlinear CSTR f_k -21.841502289865439000
sum(|constr|) 0.000000000210252562
f(x_k) + sum(|constr|) -21.841502289655185000
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 53 ConJacEv 53 Iter 41 MinorIter 340
CPU time: 0.093601 sec. Elapsed time: 0.085000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.5364 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Nonlinear CSTR f_k -21.896802275281715000
sum(|constr|) 0.000000001588897791
f(x_k) + sum(|constr|) -21.896802273692817000
f(x_0) -21.841502289865474000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 96 ConJacEv 96 Iter 91 MinorIter 378
CPU time: 0.436803 sec. Elapsed time: 0.436000 sec.

</pre>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.53229 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Nonlinear CSTR f_k -21.887245712611271000
sum(|constr|) 0.000000000382798748
f(x_k) + sum(|constr|) -21.887245712228474000
f(x_0) -21.896802275281715000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 288 ConJacEv 288 Iter 272 MinorIter 1059
CPU time: 15.802901 sec. Elapsed time: 15.853000 sec.

</pre>

<source lang="matlab">
end

t = subs(collocate(t),solution);
x1 = collocate(x1opt);
x2 = collocate(x2opt);
x3 = collocate(x3opt);
x4 = collocate(x4opt);
x5 = collocate(x5opt);
x6 = collocate(x6opt);
x7 = collocate(x7opt);
x8 = collocate(x8opt);
u1 = collocate(u1opt);
u2 = collocate(u2opt);
u3 = collocate(u3opt);
u4 = collocate(u4opt);
</source>

==Plot result==

<source lang="matlab">
figure(1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-' ...
,t,x5,'*-',t,x6,'*-',t,x7,'*-',t,x8/10,'*-');
legend('x1','x2','x3','x4','x5','x6','x7','x8/10');
title('Nonlinear CSTR state variables');

figure(2)
plot(t,u1,'+-',t,u2,'+-',t,u3,'+-',t,u4,'+-');
legend('u1','u2','u3','u4');
title('Nonlinear CSTR control');
</source>

[[File:nonlinearCSTR_01.png]]

[[File:nonlinearCSTR_02.png]]

PROPT Nondifferentiable system

2011-11-09T12:10:21Z

Mbot:

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Global Optimization of Chemical Processes using Stochastic Algorithms 1996, Julio R. Banga, Warren D. Seider

Case Study IV: Optimal control of a nondifferentiable system

==Problem description==

This problem has been studied by Thomopoulos and Papadakis who report convergence difficulties using several optimization algorithms and by Luus using Iterative Dynamic Programming. The optimal control problem is:

Find u(t) to minimize:

<math> J = x_3(t_f) </math>

Subject to:

<math> \frac{dx_1}{dt} = x_2 </math>

<math> \frac{dx_2}{dt} = -x_1-x_2+u+d </math>

<math> \frac{dx_3}{dt} = 5*x_1^2+2.5*x_2^2+0.5*u^2 </math>

with

<math> d = 100*[U(t-0.5)-U(t-0.6)] </math>

where U(t-alpha) is the unit function such that U = 0 for t - alpha < 0 and U = 1 for t - alpha > 0. Hence d is a rectangular pulse of magnitude 100 from t=0.5 until t=0.6. These authors consider t_f = 2.0s and the initial conditions:

<math> x(t_0) = [0 \ 0 \ 0]' </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
</source>

==Problem setup==

<source lang="matlab">
toms t1
p1 = tomPhase('p1', t1, 0, 0.5, 30);
toms t2
p2 = tomPhase('p2', t2, 0.5, 0.1, 20);
toms t3
p3 = tomPhase('p3', t3, 0.6, 2-0.6, 50);

x1p1 = tomState(p1,'x1p1');
x2p1 = tomState(p1,'x2p1');
x3p1 = tomState(p1,'x3p1');
up1 = tomControl(p1,'up1');

x1p2 = tomState(p2,'x1p2');
x2p2 = tomState(p2,'x2p2');
x3p2 = tomState(p2,'x3p2');
up2 = tomControl(p2,'up2');

x1p3 = tomState(p3,'x1p3');
x2p3 = tomState(p3,'x2p3');
x3p3 = tomState(p3,'x3p3');
up3 = tomControl(p3,'up3');

% Initial guess
x0 = {icollocate(p1,{x1p1 == 0;x2p1 == 0;x3p1 == 0})
icollocate(p2,{x1p2 == 0;x2p2 == 0;x3p2 == 0})
icollocate(p3,{x1p3 == 0;x2p3 == 0;x3p3 == 0})
collocate(p1,up1==-4-8*t1/0.5)
collocate(p2,up2==-12)
collocate(p3,up3==-12+14*t3/2)};

% Box constraints
cbox = {
icollocate(p1,{-100 <= x1p1 <= 100
-100 <= x2p1 <= 100
-100 <= x3p1 <= 100})
icollocate(p2,{-100 <= x1p2 <= 100
-100 <= x2p2 <= 100
-100 <= x3p2 <= 100})
icollocate(p3,{-100 <= x1p3 <= 100
-100 <= x2p3 <= 100
-100 <= x3p3 <= 100})
collocate(p1,-15 <= up1 <= 2)
collocate(p2,-15 <= up2 <= 2)
collocate(p3,-15 <= up3 <= 2)};

% Boundary constraints
cbnd = {initial(p1,{x1p1 == 0;x2p1 == 0;x3p1 == 0})
final(p3,x3p3) <= 60};

% ODEs and path constraints
ceq = {collocate(p1,{
dot(p1,x1p1) == x2p1
dot(p1,x2p1) == -x1p1-x2p1+up1
dot(p1,x3p1) == 5*x1p1.^2+2.5*x2p1.^2+0.5*up1.^2})
collocate(p2,{
dot(p2,x1p2) == x2p2
dot(p2,x2p2) == -x1p2-x2p2+up2+100
dot(p2,x3p2) == 5*x1p2.^2+2.5*x2p2.^2+0.5*up2.^2})
collocate(p3,{
dot(p3,x1p3) == x2p3
dot(p3,x2p3) == -x1p3-x2p3+up3
dot(p3,x3p3) == 5*x1p3.^2+2.5*x2p3.^2+0.5*up3.^2})};

% Objective
objective = final(p3,x3p3);

% Link phase
link = {final(p1,x1p1) == initial(p2,x1p2)
final(p1,x2p1) == initial(p2,x2p2)
final(p1,x3p1) == initial(p2,x3p2)
final(p2,x1p2) == initial(p3,x1p3)
final(p2,x2p2) == initial(p3,x2p3)
final(p2,x3p2) == initial(p3,x3p3)};
</source>

==Solve the problem==

<source lang="matlab">
options = struct;
options.name = 'Nondiff System';
constr = {cbox, cbnd, ceq, link};
solution = ezsolve(objective, constr, x0, options);

t = subs(collocate(p1,t1),solution);
t = [t;subs(collocate(p2,t2),solution)];
t = [t;subs(collocate(p3,t3),solution)];
x1 = subs(collocate(p1,x1p1),solution);
x1 = [x1;subs(collocate(p2,x1p2),solution)];
x1 = [x1;subs(collocate(p3,x1p3),solution)];
x2 = subs(collocate(p1,x2p1),solution);
x2 = [x2;subs(collocate(p2,x2p2),solution)];
x2 = [x2;subs(collocate(p3,x2p3),solution)];
x3 = subs(collocate(p1,x3p1),solution);
x3 = [x3;subs(collocate(p2,x3p2),solution)];
x3 = [x3;subs(collocate(p3,x3p3),solution)];
u = subs(collocate(p1,up1),solution);
u = [u;subs(collocate(p2,up2),solution)];
u = [u;subs(collocate(p3,up3),solution)];
</source>

<pre>
Problem type appears to be: lpcon
Time for symbolic processing: 0.34745 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Nondiff System f_k 58.065028469582899000
sum(|constr|) 0.000030760924864613
f(x_k) + sum(|constr|) 58.065059230507764000
f(x_0) 0.000000000000000000

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 1 ConstrEv 45 ConJacEv 45 Iter 38 MinorIter 1056
CPU time: 0.530403 sec. Elapsed time: 0.534000 sec.

</pre>

==Plot result==

<source lang="matlab">
subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-');
legend('x1','x2','x3');
title('Nondiff System state variable');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Nondiff System control');
</source>

[[File:nondiffSystem_01.png]]