Global Optimization of Chemical Processes using Stochastic Algorithms 1996, Julio R. Banga, Warren D. Seider
Case Study IV: Optimal control of a nondifferentiable system
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:
J = x3(tf) |
Subject to:
| = x2 |
| = −x1−x2+u+d |
| = 5*x12+2.5*x22+0.5*u2 |
with
d = 100*[U(t−0.5)−U(t−0.6)] |
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:
x(t0) = [0 0 0]′ |
Reference: [4]
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)};
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)];
Problem type appears to be: lpcon Starting numeric solver ===== * * * =================================================================== * * * TOMLAB - Tomlab Optimization Inc. Development license 999001. Valid to 2011-02-05 ===================================================================================== Problem: --- 1: Nondiff System f_k 58.065028469582764000 sum(|constr|) 0.000030760875102477 f(x_k) + sum(|constr|) 58.065059230457869000 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: 1.437500 sec. Elapsed time: 1.516000 sec.
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');