Solving Optimal Control Problems With MATLAB
Solving Optimal Control Problems With MATLAB
Solving Optimal Control Problems With MATLAB
1 Introduction
The theory of optimal control has been well developed for over forty years.
With the advances of computer technique, optimal control is now widely
used in multi-disciplinary applications such as biological systems, communi-
cation networks and socio-economic systems etc. As a result, more and more
people will benefit greatly by learning to solve the optimal control prob-
lems numerically. Realizing such growing needs, books on optimal control
put more weight on numerical methods. In retrospect, [1] was the first and
the “classic” book for studying the theory as well as many interesting cases
(time-optimal, fuel-optimal and linear quadratic regulator(LQR) problems).
Necessary conditions for various systems were derived and explicit solutions
were given when possible. Later, [2] proved to be a concise yet excellent book
with more engineering examples. One of the distinguish features of this book
is that it introduced several iterative algorithms for solving problems numer-
ically. More recently, [3] uses MATLAB to solve problems which is easier
and more precise. However, the numerical methods covered in these books
are insufficient for the wide range of problems emerging from various fields.
Especially, for those problems with free final time and nonlinear dynamics.
This tutorial shows common routines in MATLAB to solve both fixed and
free final time problems. Specifically, the following subjects are discussed
with examples:
1
3. How to reformulate the original problem as a Boundary Value Problem
(BVP) and solve it with bvp4c?
4. How to get ‘good enough’ solutions with bvp4c and under what condi-
tions will bvp4c fail to find a solution?
It should be noted that all the routines (except the steepest descent
method) discussed in this tutorial belong to the “indirect methods” category.
This means that constraints on the controls and states are not considered.1
In other words, the control can be solved in terms of states and costates and
the problem is equivalently to a BVP. The reference of all the examples used
in this tutorial are stated such that the results can be compared and verified.
2
with boundary conditions:
T
x(0) = 0 x(2) = 5 2
(b) Consider the performance measure:
Z 2
1 1 1 2
J(u) = (x1 (2) − 5)2 + (x2 (2) − 2)2 + u (t) dt
2 2 0 2
% Hamiltonian
syms p1 p2 H;
H = g + p1*Dx1 + p2*Dx2;
% Costate equations
Dp1 = -diff(H,x1);
Dp2 = -diff(H,x2);
The MATLAB commands we used here are diff and solve. diff differen-
tiates a symbolic expression and solve gives symbolic solution to algebraic
3
equations. For more details about Symbolic Math Toolbox, please refer to
[5]. Applying the necessary conditions for optimality we get two equations:2
∂H
ṗ∗i = − (3)
∂x∗i
∂H
= 0 (4)
∂u∗
The first equation gives costate equations. From the second equation, we
solve for control u in terms of states and costates. The second step is to
substitute u from (4) back to the state and costate equations to get a set of
2n first-order ordinary differential equations (ODE’s). A solution (with 2n
arbitrary coefficients) can be obtained by using the dsolve command without
any boundary conditions. The symbolic solution looks different from the one
in [2]. By simplifying the expression and rearrange the arbitrary coefficients,
it is not difficult to see that the two are the same.
% Substitute u to state equations
Dx2 = subs(Dx2, u, sol_u);
sol_h = dsolve(eq1,eq2,eq3,eq4);
As stated in [2], the differences of the three cases in this problem are
merely the boundary conditions. For (a), the arbitrary coefficients can be
determined by supplying the 2n boundary conditions to dsovle:
% case a: (a) x1(0)=x2(0)=0; x1(2) = 5; x2(2) = 2;
conA1 = ’x1(0) = 0’;
conA2 = ’x2(0) = 0’;
conA3 = ’x1(2) = 5’;
conA4 = ’x2(2) = 2’;
sol_a = dsolve(eq1,eq2,eq3,eq4,conA1,conA2,conA3,conA4);
Again the solutions given by MATLAB and [2] look different from each
other. Yet Figure 1 shows that the two are in fact equivalent. For all the
figures in this problem, * represent state trajectory from [2] while symbolic
solution from MATLAB is plotted with a continuous line.
2
We use * to indicate the optimal state trajectory or control.
4
Figure 1: Trajectories for Example 1 (a)
sol_b = solve(eq1b,eq2b,eq3b,eq4b);
5
Figure 2: Trajectories for Example 1 (b)
Case (c) is almost the same as case (b) with slightly different bound-
ary conditions. From [2], the boundary conditions are: x∗1 (2) + 5x∗2 (2) =
15, p∗2 (2) = 5p∗1 (2). Figure 3 shows the results of both solutions from MAT-
LAB and [2].
% case c: x1(0)=x2(0)=0;x1(2)+5*x2(2)=15;p2(2)= 5*p1(2);
eq1c = char(subs(sol_h.x1,’t’,0));
eq2c = char(subs(sol_h.x2,’t’,0));
eq3c = strcat(char(subs(sol_h.p2,’t’,2)),...
’-(’,char(subs(sol_h.p1,’t’,2)),’)*5’);
eq4c = strcat(char(subs(sol_h.x1,’t’,2)),...
’+(’,char(subs(sol_h.x2,’t’,2)),’)*5-15’);
sol_c = solve(eq1c,eq2c,eq3c,eq4c);
% Substitute the coefficients
C1 = double(sol_c.C1);
6
Figure 3: Trajectories for Example 1 (c)
C2 = double(sol_c.C2);
C3 = double(sol_c.C3);
C4 = double(sol_c.C4);
sol_c2 = struct(’x1’,{subs(sol_h.x1)},’x2’,{subs(sol_h.x2)}, ...
’p1’,{subs(sol_h.p1)},’p2’,{subs(sol_h.p2)});
7
C(t) is the deviation from the steady-state concentration. u(t) is the normal-
ized control variable, representing the effect of coolant flow on the chemical
reaction.
The performance measure to be minimized is:
Z 0.78
J= [x21 (t) + x22 (t) + Ru2 (t)] dt,
0
First, we will implement the steepest descent method based on the scheme
outlined in [2]. The algorithm consists of 4 steps:
1. Subdivide the interval [t0 , tf ] into N equal subintervals and assume
a piecewise-constant control u(0) (t) = u(0) (tk ), t ∈ [tk , tk+1 ] k =
0, 1, · · · , N − 1
2. Applying the assumed control u(i) to integrate the state equations from
t0 to tf with initial conditions x(t0 ) = x0 and store the state trajectory
x(i) .
3. Applying u(i) and x(i) to integrate costate equations backward, i.e.,
from [tf , t0 ]. The “initial value” p(i) (tf ) can be obtained by:
∂h (i)
p(i) (tf ) =
x (tf ) .
∂x
Evaluate ∂H (i) (t)/∂u, t ∈ [t0 , tf ] and store this vector.
4. If
∂H (i)
∂u
≤ γ (6)
8
The main loop in MATLAB is as follows:
for i = 1:max_iteration
% 1) start with assumed control u and move forward
[Tx,X] = ode45(@(t,x) stateEq(t,x,u,Tu), [t0 tf], ...
initx, options);
9
dx(1) = -2*(x(1) + 0.25) + (x(2) + 0.5)*exp(25*x(1)/(x(1)+2)) ...
- (x(1) + 0.25).*u;
dx(2) = 0.5 - x(2) -(x(2) + 0.5)*exp(25*x(1)/(x(1)+2));
% Costate equations
function dp = costateEq(t,p,u,Tu,x1,x2,xt)
dp = zeros(2,1);
x1 = interp1(xt,x1,t); % Interploate the state varialbes
x2 = interp1(xt,x2,t);
u = interp1(Tu,u,t); % Interploate the control
dp(1) = p(1).*(u + exp((25*x1)/(x1 + 2)).*((25*x1)/(x1 + 2)^2 - ...
25/(x1 + 2))*(x2 + 1/2) + 2) - ...
2*x1 - p(2).*exp((25*x1)/(x1 + 2))*((25*x1)/(x1 + 2)^2 - ...
25/(x1 + 2))*(x2 + 1/2);
dp(2) = p(2).*(exp((25*x1)/(x1 + 2)) + 1) - ...
p(1).*exp((25*x1)/(x1 + 2)) - 2*x2;
The step size τ is set as a constant in this example by some post hoc
investigation. A better strategy is to select τ with a line search method which
will maximize the reduction of performance measure with given ∂H (i) /∂u in
each iteration. Figure 4 shows the optimal state trajectory and control over
the time. The value of performance measure as a function of iteration number
is shown in Figure 5. In [2], two more numerical algorithms were introduced:
variation of extremals and quasilinearization. These two methods basically
reformulate and solve the original problem as a Boundary Value Problem
(BVP). In MATLAB, a BVP is typically solved with bvp4c. [4] is an excellent
reference on using bvp4c. For fix-final-time problems, u can always be solved
10
Figure 4: Example 2 Steepest descent method
with respect to x and p by applying the necessary conditions (3) and (4).
And we will have 2n ODE’s and 2n boundary conditions.
11
Figure 6: Solution from bfv4c for Problem 2
12
function dydt = BVP_ode(t,y)
global R;
t1 = y(1)+.25;
t2 = y(2)+.5;
t3 = exp(25*y(1)/(y(2)+2));
t4 = 50/(y(1)+2)^2;
u = y(3)*t1/(2*R);
dydt = [-2*t1+t2*t3-t2*u
0.5-y(2)-t2*t3
-2*y(1)+2*y(3)-y(3)*t2*t4*t3+y(3)*u+y(4)*t2*t4*t3
-2*y(2)-y(3)*t3+y(4)*(1+t3)];
% -----------------------------------------------
% The boundary conditions:
% x1(0) = 0.05, x2(0) = 0, tf = 0.78, p1(tf) = 0, p2(tf) = 0;
function res = BVP_bc(ya,yb)
res = [ ya(1) - 0.05
ya(2) - 0
yb(3) - 0
yb(4) - 0 ];
In this example, bvp4c works perfectly. It is faster and gives better re-
sults, i.e. a smaller performance measure J comparing to the steepest descent
method (see Figure 6). In the following section, we will solely use bvp4c when
numerical solutions are needed.
13
find the optimal control given the boundary conditions as:
T
x(0) = 1 2 , x1 (tf ) = 3, x2 (tf ) is free
To use the Symbolic Math Toolbox, the routine is very similar to Problem
1. We first supply the ODE’s and boundary conditions on states and costates
to dsolve. The only difference is that the final time tf itself is now a variable.
As a result, the solution is a function of tf . Next, we introduce four more
variables, namely x1 (tf ), x2 (tf ), p1 (tf ), p2 (tf ) into the solution obtained
above. With one additional boundary condition from
∂h ∗
H (x∗ (tf ), u∗ (tf ), p∗ (tf ), tf ) + (x (tf ), tf ) = 0
∂t
For this problem, h ≡ 0 and we have p1 (tf )x2 (tf ) − 0.5p22 (tf ) = 0. Now we
have 5 algebraic equations with 5 unknowns. And solve comes in handy
to solve this problem. Figure 7 shows the results from MATLAB and the
analytical solution in [3]. Although Symbolic Math Toolbox works fine in this
example, it should be pointed out that in most problems, it is impossible to
get explicit solutions. For example, there is no explicit solutions for Example
1 even though the state equations are similar to those of Example 3.
sol = dsolve(’Dx1 = x2, Dx2 = -p2, Dp1 = 0, Dp2 = -p1’,...
’x1(0) = 1, x2(0) = 2, x1(tf) = 3, p2(tf) = 0’);
eq1 = subs(sol.x1) - ’x1tf’;
eq2 = subs(sol.x2) - ’x2tf’;
eq3 = subs(sol.p1) - ’p1tf’;
eq4 = subs(sol.p2) - ’p2tf’;
eq5 = sym(’p1tf*x2tf - 0.5*p2tf^2’);
%%
sol_2 = solve(eq1, eq2, eq3, eq4, eq5);
tf = sol_2.tf;
x1tf = sol_2.x1tf;
x2tf = sol_2.x2tf;
x1 = subs(sol.x1);
x2 = subs(sol.x2);
p1 = subs(sol.p1);
p2 = subs(sol.p2);
14
Figure 7: Example 3 Symbolic method
a problem: the time interval is not known. One common treatment [4] [7]
for such a situation is to change the independent variable t to τ = t/T , the
augmented state and costate equations will then become x̃˙ = T f (x, q, τ ).3
Now the problem is posed on fixed interval [ 0, 1 ]. This can be implemented
in bvp4c by treating T as an auxiliary variable. The following code snippet
shows the details.
solinit = bvpinit(linspace(0,1),[2;3;1;1;2]);
% -------------------------------------------------------------------------
% ODE’s of augmented states
function dydt = ode(t,y)
dydt = y(5)*[ y(2);-y(4);0;-y(3);0 ];
% -------------------------------------------------------------------------
% boundary conditions: x1(0)=1;x2(0)=2, x1(tf)=3, p2(tf)=0;
3
f denotes the ODE’s for state and costates.
15
% p1(tf)*x2(tf)-0.5*p2(2)^2
function res = bc(ya,yb)
res = [ ya(1) - 1; ya(2) - 2; yb(1) - 3; yb(4);
yb(3)*yb(2)-0.5*yb(4)^2];
% -----------------------------------------------------------
% ODE’s of augmented states
function dydt = ode(t,y,T)
dydt = T*[ y(2);-y(4);0;-y(3) ];
% -----------------------------------------------------------
% boundary conditions: x1(0)=1;x2(0)=2, x1(tf)=3, p2(tf)=0;
% p1(tf)*x2(tf)-0.5*p2(2)^2
function res = bc(ya,yb,T)
res = [ ya(1) - 1; ya(2) - 2; yb(1) - 3; yb(4);
yb(3)*yb(2)-0.5*yb(4)^2];
16
Figure 8: Example 3 Numerical method
17
The boundary conditions are:
sol = bvp4c(@ode,@bc,solinit);
...
% the solution for one value of mu is used as guess for the next.
for i=2:3
if i==2
mu = 0.3;
else
mu = .1;
end;
% After creating function handles, the new value of mu
% will be used in nested functions.
sol = bvp4c(@ode,@bc,sol);
...
end
% -----------------------------------------------------------
function dydt = ode(t,y,T)
global mu;
term = sqrt(mu*y(3)^2+y(4)^2);
dydt = T*[ y(2) + mu*y(3)/term
y(4)/term
0
-y(3)];
% ------------------------------------------------------------
% boundary conditions, with 4 states and 1 parameters, 5 conditions are
% needed: x1(0) =2, x2(0) = 2; x1(1) = 0; x2(1) = 0; x3(1)^2+x4(1)^2 = 1;
function res = bc(ya,yb,T)
res = [ya(1)-2
ya(2)-2
18
Figure 9: Example 4 State trajectory
yb(1)
yb(2)
yb(3)^2+yb(4)^2-1];
19
Figure 10: Example 4 Control
4 Discussion
For many optimal control problems, bvp4c is the best option we have. Once
the problem is reformatted into a BVP, bvp4c can usually solve it efficiently.
However, we must pay attention to the limitations of bvp4c, more specifically
• Good initial guess is important to get an accurate solution, if a solution
exists.
20
and improves with time. More information on GPOPS can be found in
http://www.gpops.org/.
References
[1] M. Athans and P. Falb (2007). Optimal control: an introduction to the
theory and its applications. Dover Publications, Inc.
[6] S. N. Avvakumov and Yu.N. Kiselev Boundary value problem for or-
dinary differential equations with applications to optimal control.
21