PROPT Bryson-Denham Two Phase Problem: Difference between revisions
From TomWiki
Jump to navigationJump to search
No edit summary |
|||
Line 202: | Line 202: | ||
[[File:brysonDenhamTwoPhase_01.png]] | [[File:brysonDenhamTwoPhase_01.png]] | ||
[[Category:PROPT Examples]] |
Latest revision as of 04:55, 14 February 2012
This page is part of the PROPT Manual. See PROPT Manual. |
An example of how the input could look for PROPT.
Problem description
In this example we also take advantage of the advance knowledge that the solution reaches x1=x1max with x2=0, to introduce an event that divides the time interval into two phases. This increases the accuracy of the solution.
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.
Problem setup
for n=[5 21]
toms t1 tcut
p1 = tomPhase('p1', t1, 0, tcut, n);
toms t2 tmax
p2 = tomPhase('p2', t2, tcut, tmax-tcut, n);
setPhase(p1);
tomStates x1p1 x2p1
tomControls up1
setPhase(p2);
tomStates x1p2 x2p2
tomControls up2
% Constant
x1max = 1/9;
setPhase(p1);
% Initial guess
if n==5
x01 = {tcut == 0.25
tmax == 0.5
icollocate({
x1p1 == 0
x2p1 == 1-2*t1/tcut
})
collocate(up1==0)};
else
x01 = {tcut == tcut_opt
tmax == tmax_opt
icollocate({
x1p1 == x1p1_opt
x2p1 == x2p1_opt
})
collocate(up1==up1_opt)};
end
% Box constraints
cbox1 = {0.001 <= tcut <= tmax-0.01
tmax <= 50
collocate({0 <= x1p1 <= x1max
-10 <= x2p1 <= 10})};
% Set up initial conditions in phase p1
% Initial constraints
cbnd1 = initial({x1p1 == 0; x2p1 == 1});
% ODEs and path constraints
ceq1 = collocate({
dot(x1p1) == x2p1
dot(x2p1) == up1});
% We take advantage of the fact that we've determined that a good place to
% shift between phases is when x1 reaches x1max, and that x2 must equal 0
% there (Later, we want the solver to be able to figure this out for
% itself).
% Final constraints
cbnd1 = {cbnd1
final({x1p1 == x1max; x2p1 == 0})};
% Using integral gives the integral over the phase of an expression -
% in this case 0.5 times the square of u.
% Objective
objective1 = integrate(0.5*up1.^2);
setPhase(p2);
% Initial guess
if n==5
x02 = {icollocate({
x1p2 == 0
x2p2 == 1-2*t2/tmax
})
collocate(up2==0)};
else
x02 = {icollocate({
x1p2 == x1p2_opt
x2p2 == x2p2_opt
})
collocate(up2==up2_opt)};
end
% Box constraints
cbox2 = collocate({0 <= x1p2 <= x1max
-10 <= x2p2 <= 10});
% ODEs and path constraints
ceq2 = collocate({
dot(x1p2) == x2p2
dot(x2p2) == up2});
% x2_i of p2 is already linked to x2_f of p1, but linking it to a constant
% helps convergence.
% Final conditions in phase p2.
cbnd2 = {initial(x2p2 == 0)
final(x1p2 == 0)
final(x2p2 == -1)};
objective2 = integrate(0.5*up2.^2);
% Link the phases
link = {final(p1,x1p1) == initial(p2,x1p2)
final(p1,x2p1) == initial(p2,x2p2)};
Solve the problem
options = struct;
options.name = 'Bryson Denham Two Phase';
objective = objective1+objective2;
constr = {cbox1, cbnd1, ceq1, cbox2, cbnd2, ceq2, link};
solution = ezsolve(objective, constr, {x01, x02}, options);
x1p1_opt = subs(x1p1, solution);
x2p1_opt = subs(x2p1, solution);
up1_opt = subs(up1, solution);
x1p2_opt = subs(x1p2, solution);
x2p2_opt = subs(x2p2, solution);
up2_opt = subs(up2, solution);
tcut_opt = subs(tcut, solution);
tmax_opt = subs(tmax, solution);
Problem type appears to be: con Time for symbolic processing: 0.20787 seconds Starting numeric solver ===== * * * =================================================================== * * * TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31 ===================================================================================== Problem: --- 1: Bryson Denham Two Phase f_k 3.999999993412657800 sum(|constr|) 0.000000023172966002 f(x_k) + sum(|constr|) 4.000000016585623500 f(x_0) 0.000000000000000000 Solver: snopt. EXIT=0. INFORM=1. SNOPT 7.2-5 NLP code Optimality conditions satisfied FuncEv 39 GradEv 37 ConstrEv 37 ConJacEv 37 Iter 34 MinorIter 70 CPU time: 0.046800 sec. Elapsed time: 0.037000 sec.
Problem type appears to be: con Time for symbolic processing: 0.21006 seconds Starting numeric solver ===== * * * =================================================================== * * * TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31 ===================================================================================== Problem: --- 1: Bryson Denham Two Phase f_k 3.999999993239003300 sum(|constr|) 0.000000001580940421 f(x_k) + sum(|constr|) 3.999999994819943600 f(x_0) 3.999999738838719000 Solver: snopt. EXIT=0. INFORM=1. SNOPT 7.2-5 NLP code Optimality conditions satisfied FuncEv 4 GradEv 2 ConstrEv 2 ConJacEv 2 Iter 1 MinorIter 76 CPU time: 0.015600 sec. Elapsed time: 0.017000 sec.
end
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)];
Plot the result
figure(1)
plot(t,x1,'*-',t,x2,'*-');
title('Bryson Denham Two Phase state variables');