Analytical Direct Shooting and Problem Simplification
This section revisits a minimum-time optimal control problem to emphasize an important modeling principle: not every problem requires heavy numerical machinery. Before writing code, the governing equations and optimality conditions should be simplified as far as possible analytically.
The example considered here is
min u ( ⋅ ) , t f J = t f \min_{u(\cdot),\,t_f} J=t_f u ( ⋅ ) , t f min J = t f subject to
x ˙ = v sin u , y ˙ = − v cos u , v ˙ = g cos u , \begin{aligned}
\dot{x} &= v\sin u,\\
\dot{y} &= -v\cos u,\\
\dot{v} &= g\cos u,
\end{aligned} x ˙ y ˙ v ˙ = v sin u , = − v cos u , = g cos u , with boundary conditions
x ( 0 ) = 0 , y ( 0 ) = 0 , v ( 0 ) = 0 , x ( t f ) = 2 , y ( t f ) = − 2 , \begin{aligned}
x(0)&=0, & y(0)&=0, & v(0)&=0,\\
x(t_f)&=2, & y(t_f)&=-2,
\end{aligned} x ( 0 ) x ( t f ) = 0 , = 2 , y ( 0 ) y ( t f ) = 0 , = − 2 , v ( 0 ) = 0 , and free terminal velocity.
Hamiltonian and Costate Equations ¶ The Hamiltonian is
H = 1 + λ x v sin u − λ y v cos u + λ v g cos u . \mathcal{H}
=
1
+
\lambda_x v\sin u
-
\lambda_y v\cos u
+
\lambda_v g\cos u. H = 1 + λ x v sin u − λ y v cos u + λ v g cos u . The costate equations are
λ ˙ x = − ∂ H ∂ x = 0 , λ ˙ y = − ∂ H ∂ y = 0 , λ ˙ v = − ∂ H ∂ v = − λ x sin u + λ y cos u . \begin{aligned}
\dot{\lambda}_x
&=
-\frac{\partial \mathcal{H}}{\partial x}
=
0,\\
\dot{\lambda}_y
&=
-\frac{\partial \mathcal{H}}{\partial y}
=
0,\\
\dot{\lambda}_v
&=
-\frac{\partial \mathcal{H}}{\partial v}
=
-\lambda_x\sin u+\lambda_y\cos u.
\end{aligned} λ ˙ x λ ˙ y λ ˙ v = − ∂ x ∂ H = 0 , = − ∂ y ∂ H = 0 , = − ∂ v ∂ H = − λ x sin u + λ y cos u . Thus,
λ x = constant , λ y = constant . \lambda_x=\text{constant},
\qquad
\lambda_y=\text{constant}. λ x = constant , λ y = constant . Stationarity Condition ¶ For an interior optimal control,
∂ H ∂ u = 0. \frac{\partial \mathcal{H}}{\partial u}=0. ∂ u ∂ H = 0. Hence,
λ x v cos u + λ y v sin u − λ v g sin u = 0. \lambda_x v\cos u
+
\lambda_y v\sin u
-
\lambda_v g\sin u
=
0. λ x v cos u + λ y v sin u − λ v g sin u = 0. Rearranging,
λ x v cos u = ( λ v g − λ y v ) sin u . \lambda_x v\cos u
=
\left(
\lambda_v g-\lambda_y v
\right)
\sin u. λ x v cos u = ( λ v g − λ y v ) sin u . Therefore,
tan u = λ x v λ v g − λ y v . \boxed{
\tan u
=
\frac{\lambda_x v}
{\lambda_v g-\lambda_y v}.
} tan u = λ v g − λ y v λ x v . By differentiating this relation and using the state and costate equations, one obtains
u ¨ = 0. \boxed{
\ddot{u}=0.
} u ¨ = 0. Therefore, the optimal control is affine in time:
u ( t ) = c 0 + c 1 t . \boxed{
u(t)=c_0+c_1 t.
} u ( t ) = c 0 + c 1 t . This structural result is the key simplification.
Direct-Shooting Parameterization ¶ A conventional direct-shooting implementation would parameterize u ( t ) u(t) u ( t )
Here, however, the affine control allows the dynamics to be integrated analytically.
The unknown parameter vector is
p = [ c 0 c 1 t f ] T . \boxed{
\boldsymbol{p}
=
\begin{bmatrix}
c_0 & c_1 & t_f
\end{bmatrix}^{\mathsf T}.
} p = [ c 0 c 1 t f ] T . Analytical Integration of the Velocity ¶ Using
u ( t ) = c 0 + c 1 t , u(t)=c_0+c_1t, u ( t ) = c 0 + c 1 t , the velocity equation becomes
v ˙ = g cos ( c 0 + c 1 t ) . \dot{v}
=
g\cos(c_0+c_1t). v ˙ = g cos ( c 0 + c 1 t ) . Integrating from 0 to t t t
v ( t ) = v ( 0 ) + ∫ 0 t g cos ( c 0 + c 1 τ ) d τ = g c 1 [ sin ( c 0 + c 1 t ) − sin c 0 ] . \begin{aligned}
v(t)
&=
v(0)
+
\int_0^t
g\cos(c_0+c_1\tau)\,\mathrm{d}\tau\\
&=
\frac{g}{c_1}
\left[
\sin(c_0+c_1t)-\sin c_0
\right].
\end{aligned} v ( t ) = v ( 0 ) + ∫ 0 t g cos ( c 0 + c 1 τ ) d τ = c 1 g [ sin ( c 0 + c 1 t ) − sin c 0 ] . Since v ( 0 ) = 0 v(0)=0 v ( 0 ) = 0
v ( t ) = g c 1 [ sin ( c 0 + c 1 t ) − sin c 0 ] . \boxed{
v(t)
=
\frac{g}{c_1}
\left[
\sin(c_0+c_1t)-\sin c_0
\right].
} v ( t ) = c 1 g [ sin ( c 0 + c 1 t ) − sin c 0 ] . The case c 1 = 0 c_1=0 c 1 = 0 c 1 c_1 c 1
Analytical Integration of the Horizontal Position ¶ The horizontal dynamics are
x ˙ = v sin u . \dot{x}
=
v\sin u. x ˙ = v sin u . Substituting v ( t ) v(t) v ( t ) u ( t ) u(t) u ( t )
x ˙ = g c 1 [ sin ( c 0 + c 1 t ) − sin c 0 ] sin ( c 0 + c 1 t ) . \dot{x}
=
\frac{g}{c_1}
\left[
\sin(c_0+c_1t)-\sin c_0
\right]
\sin(c_0+c_1t). x ˙ = c 1 g [ sin ( c 0 + c 1 t ) − sin c 0 ] sin ( c 0 + c 1 t ) . Expanding,
x ˙ = g c 1 sin 2 ( c 0 + c 1 t ) − g sin c 0 c 1 sin ( c 0 + c 1 t ) . \dot{x}
=
\frac{g}{c_1}
\sin^2(c_0+c_1t)
-
\frac{g\sin c_0}{c_1}
\sin(c_0+c_1t). x ˙ = c 1 g sin 2 ( c 0 + c 1 t ) − c 1 g sin c 0 sin ( c 0 + c 1 t ) . Using
sin 2 α = 1 − cos 2 α 2 , \sin^2\alpha
=
\frac{1-\cos 2\alpha}{2}, sin 2 α = 2 1 − cos 2 α , we obtain
x ˙ = g 2 c 1 − g 2 c 1 cos ( 2 c 0 + 2 c 1 t ) − g sin c 0 c 1 sin ( c 0 + c 1 t ) . \dot{x}
=
\frac{g}{2c_1}
-
\frac{g}{2c_1}
\cos(2c_0+2c_1t)
-
\frac{g\sin c_0}{c_1}
\sin(c_0+c_1t). x ˙ = 2 c 1 g − 2 c 1 g cos ( 2 c 0 + 2 c 1 t ) − c 1 g sin c 0 sin ( c 0 + c 1 t ) . Integrating and enforcing x ( 0 ) = 0 x(0)=0 x ( 0 ) = 0
x ( t ) = g t 2 c 1 − g 4 c 1 2 [ sin ( 2 c 0 + 2 c 1 t ) − sin ( 2 c 0 ) ] + g sin c 0 c 1 2 [ cos ( c 0 + c 1 t ) − cos c 0 ] . \begin{aligned}
x(t)
={}&
\frac{gt}{2c_1}
-
\frac{g}{4c_1^2}
\left[
\sin(2c_0+2c_1t)-\sin(2c_0)
\right]
\nonumber\\
&+
\frac{g\sin c_0}{c_1^2}
\left[
\cos(c_0+c_1t)-\cos c_0
\right].
\end{aligned} x ( t ) = 2 c 1 g t − 4 c 1 2 g [ sin ( 2 c 0 + 2 c 1 t ) − sin ( 2 c 0 ) ] + c 1 2 g sin c 0 [ cos ( c 0 + c 1 t ) − cos c 0 ] . Therefore,
x ( t ) = g t 2 c 1 − g 4 c 1 2 [ sin ( 2 c 0 + 2 c 1 t ) − sin 2 c 0 ] + g sin c 0 c 1 2 [ cos ( c 0 + c 1 t ) − cos c 0 ] . \boxed{
x(t)
=
\frac{gt}{2c_1}
-
\frac{g}{4c_1^2}
\left[
\sin(2c_0+2c_1t)-\sin 2c_0
\right]
+
\frac{g\sin c_0}{c_1^2}
\left[
\cos(c_0+c_1t)-\cos c_0
\right].
} x ( t ) = 2 c 1 g t − 4 c 1 2 g [ sin ( 2 c 0 + 2 c 1 t ) − sin 2 c 0 ] + c 1 2 g sin c 0 [ cos ( c 0 + c 1 t ) − cos c 0 ] . Analytical Integration of the Vertical Position ¶ The vertical dynamics are
y ˙ = − v cos u . \dot{y}
=
-v\cos u. y ˙ = − v cos u . Substituting,
y ˙ = − g c 1 [ sin ( c 0 + c 1 t ) − sin c 0 ] cos ( c 0 + c 1 t ) . \dot{y}
=
-\frac{g}{c_1}
\left[
\sin(c_0+c_1t)-\sin c_0
\right]
\cos(c_0+c_1t). y ˙ = − c 1 g [ sin ( c 0 + c 1 t ) − sin c 0 ] cos ( c 0 + c 1 t ) . Thus,
y ˙ = − g c 1 sin ( c 0 + c 1 t ) cos ( c 0 + c 1 t ) + g sin c 0 c 1 cos ( c 0 + c 1 t ) . \begin{aligned}
\dot{y}
={}&
-\frac{g}{c_1}
\sin(c_0+c_1t)\cos(c_0+c_1t)
\nonumber\\
&+
\frac{g\sin c_0}{c_1}
\cos(c_0+c_1t).
\end{aligned} y ˙ = − c 1 g sin ( c 0 + c 1 t ) cos ( c 0 + c 1 t ) + c 1 g sin c 0 cos ( c 0 + c 1 t ) . Integrating and enforcing y ( 0 ) = 0 y(0)=0 y ( 0 ) = 0
y ( t ) = − g 2 c 1 2 [ sin 2 ( c 0 + c 1 t ) − sin 2 c 0 ] + g sin c 0 c 1 2 [ sin ( c 0 + c 1 t ) − sin c 0 ] . \begin{aligned}
y(t)
={}&
-\frac{g}{2c_1^2}
\left[
\sin^2(c_0+c_1t)-\sin^2 c_0
\right]
\nonumber\\
&+
\frac{g\sin c_0}{c_1^2}
\left[
\sin(c_0+c_1t)-\sin c_0
\right].
\end{aligned} y ( t ) = − 2 c 1 2 g [ sin 2 ( c 0 + c 1 t ) − sin 2 c 0 ] + c 1 2 g sin c 0 [ sin ( c 0 + c 1 t ) − sin c 0 ] . Therefore,
y ( t ) = − g 2 c 1 2 [ sin 2 ( c 0 + c 1 t ) − sin 2 c 0 ] + g sin c 0 c 1 2 [ sin ( c 0 + c 1 t ) − sin c 0 ] . \boxed{
y(t)
=
-\frac{g}{2c_1^2}
\left[
\sin^2(c_0+c_1t)-\sin^2 c_0
\right]
+
\frac{g\sin c_0}{c_1^2}
\left[
\sin(c_0+c_1t)-\sin c_0
\right].
} y ( t ) = − 2 c 1 2 g [ sin 2 ( c 0 + c 1 t ) − sin 2 c 0 ] + c 1 2 g sin c 0 [ sin ( c 0 + c 1 t ) − sin c 0 ] . Terminal Constraints ¶ The terminal conditions require
x ( t f ; c 0 , c 1 ) = 2 , y ( t f ; c 0 , c 1 ) = − 2. \begin{aligned}
x(t_f;c_0,c_1)&=2,\\
y(t_f;c_0,c_1)&=-2.
\end{aligned} x ( t f ; c 0 , c 1 ) y ( t f ; c 0 , c 1 ) = 2 , = − 2. Define
F 1 ( c 0 , c 1 , t f ) = x ( t f ; c 0 , c 1 ) − 2 , F 2 ( c 0 , c 1 , t f ) = y ( t f ; c 0 , c 1 ) + 2. \begin{aligned}
F_1(c_0,c_1,t_f)
&=
x(t_f;c_0,c_1)-2,\\
F_2(c_0,c_1,t_f)
&=
y(t_f;c_0,c_1)+2.
\end{aligned} F 1 ( c 0 , c 1 , t f ) F 2 ( c 0 , c 1 , t f ) = x ( t f ; c 0 , c 1 ) − 2 , = y ( t f ; c 0 , c 1 ) + 2. The optimization problem becomes
min c 0 , c 1 , t f t f \boxed{
\min_{c_0,c_1,t_f}
t_f
} c 0 , c 1 , t f min t f subject to
F 1 ( c 0 , c 1 , t f ) = 0 , F 2 ( c 0 , c 1 , t f ) = 0. \boxed{
F_1(c_0,c_1,t_f)=0,
\qquad
F_2(c_0,c_1,t_f)=0.
} F 1 ( c 0 , c 1 , t f ) = 0 , F 2 ( c 0 , c 1 , t f ) = 0. No numerical integration of the state equations is required.
Degrees-of-Freedom Analysis ¶ There are three unknowns:
c 0 , c 1 , t f . c_0,\quad c_1,\quad t_f. c 0 , c 1 , t f . There are two terminal equality constraints:
F 1 = 0 , F 2 = 0. F_1=0,
\qquad
F_2=0. F 1 = 0 , F 2 = 0. Hence, the number of free degrees of freedom is
3 − 2 = 1. \boxed{
3-2=1.
} 3 − 2 = 1. That remaining degree of freedom allows minimization of t f t_f t f
If t f t_f t f
The analytical reduction provides several advantages:
no numerical ODE integration;
no integration tolerances;
no accumulated state-integration error;
fewer sources of implementation error;
faster optimization;
clearer interpretation of the unknowns;
easier verification.
The main lesson is not that numerical integration is wrong. The lesson is that it should not be used when the problem can be simplified exactly.
Treatment of the Case c 1 = 0 c_1=0 c 1 = 0 ¶ If c 1 = 0 c_1=0 c 1 = 0
u ( t ) = c 0 . u(t)=c_0. u ( t ) = c 0 . The velocity becomes
v ( t ) = g t cos c 0 . v(t)=gt\cos c_0. v ( t ) = g t cos c 0 . Then
x ( t ) = 1 2 g t 2 sin c 0 cos c 0 , y ( t ) = − 1 2 g t 2 cos 2 c 0 . \begin{aligned}
x(t)
&=
\frac{1}{2}gt^2\sin c_0\cos c_0,\\
y(t)
&=
-\frac{1}{2}gt^2\cos^2 c_0.
\end{aligned} x ( t ) y ( t ) = 2 1 g t 2 sin c 0 cos c 0 , = − 2 1 g t 2 cos 2 c 0 . This special case can be checked separately against the terminal constraints.
function solve_session31_part1
g = 9.81;
z0 = [0.5; 0.5; 1.0]; % [c0, c1, tf]
obj = @(z) z(3);
nonlcon = @(z) terminal_constraints(z,g);
options = optimoptions('fmincon', ...
'Display','iter', ...
'Algorithm','sqp', ...
'ConstraintTolerance',1e-10, ...
'OptimalityTolerance',1e-10);
lb = [-pi; -10; 1e-3];
ub = [ pi; 10; 20];
z = fmincon(obj,z0,[],[],[],[],lb,ub,nonlcon,options);
c0 = z(1);
c1 = z(2);
tf = z(3);
fprintf('c0 = %.12f\n',c0);
fprintf('c1 = %.12f\n',c1);
fprintf('tf = %.12f\n',tf);
end
function [c,ceq] = terminal_constraints(z,g)
c0 = z(1);
c1 = z(2);
tf = z(3);
if abs(c1) < 1e-8
x = 0.5*g*tf^2*sin(c0)*cos(c0);
y = -0.5*g*tf^2*cos(c0)^2;
else
x = g*tf/(2*c1) ...
- g/(4*c1^2)*(sin(2*c0+2*c1*tf)-sin(2*c0)) ...
+ g*sin(c0)/c1^2*(cos(c0+c1*tf)-cos(c0));
y = -g/(2*c1^2)*(sin(c0+c1*tf)^2-sin(c0)^2) ...
+ g*sin(c0)/c1^2*(sin(c0+c1*tf)-sin(c0));
end
c = [];
ceq = [x-2; y+2];
endVerification ¶ After obtaining ( c 0 , c 1 , t f ) (c_0,c_1,t_f) ( c 0 , c 1 , t f )
verify the terminal constraints;
reconstruct u ( t ) u(t) u ( t )
reconstruct v ( t ) v(t) v ( t ) x ( t ) x(t) x ( t ) y ( t ) y(t) y ( t )
optionally integrate the original dynamics numerically;
compare the numerical and analytical trajectories.
General Lessons for Direct Shooting ¶ Derive the structure of the optimal control first.
Use analytical integration whenever possible.
Count unknowns, equations, and degrees of freedom.
Do not use a numerical solver as a substitute for analysis.
Formulate the smallest nonlinear program that represents the problem.
Connection. When the reduced Hamiltonian is linear in a bounded control, switching functions distinguish bang arcs from possible singular intervals.