Hamilton函数方法是变分法应用在控制系统上的标准化方法，即使不懂变分法，简单套用表格中的公式也可以列写出方程，这个方法是最优控制理论用的最多的方法。
本篇博客是本系列最长也是最重要的一篇，持续更新，欢迎同学和朋友们提出修改建议。
Hamiltonian 目录
1. 规范化的最优控制问题
按照第一章最优控制理论 一、变分法和泛函极值问题，我们已经讨论了有动力学方程约束
f
(
x
,
x
˙
,
t
)
=
0
\boldsymbol f(\boldsymbol x,\dot {\boldsymbol x},t)=0
f(x,x˙,t)=0的动态系统，若无其他约束，这个系统的最优轨线遵循以下必要条件
H
x
−
d
d
t
H
x
˙
=
0
f
(
x
(
t
)
,
x
(
t
)
˙
,
t
)
=
0
\begin{aligned} H_x-\frac{\text d}{\text d t}H_{\dot x}=0\\ \boldsymbol f(\boldsymbol x(t),\dot{\boldsymbol x(t)},t)=0 \end{aligned}
Hx−dtdHx˙=0f(x(t),x(t)˙,t)=0
其中的Hamilton函数
H
(
x
,
x
˙
,
λ
,
t
)
=
L
(
x
,
x
˙
,
t
)
−
λ
T
f
(
x
,
x
˙
,
t
)
H(x,\dot x,\lambda,t)=L(x,\dot x,t)-\lambda^{\mathrm T}f(x,\dot x,t)
H(x,x˙,λ,t)=L(x,x˙,t)−λTf(x,x˙,t)。控制系统中更常见的一阶非线性系统方程，问题是这样的：
t
f
t_f
tf给定，终端状态未知或已知（仅边界条件不同），除状态方程外没有约束，且
x
˙
=
f
[
x
(
t
)
,
u
(
t
)
,
t
]
;
x
(
t
o
)
=
x
0
t
o
≤
t
≤
t
f
min
u
(
t
)
J
=
φ
[
x
(
t
f
)
,
t
f
]
+
∫
t
o
t
f
L
[
x
(
t
)
,
u
(
t
)
,
t
]
d
t
(1)
\dot{x}=f[x(t), u(t), t] ; \quad x\left(t_{o}\right)=x_0\quad t_{o} \leq t \leq t_{f}\\ \min_{u(t)}J=\varphi\left[x\left(t_{f}\right), t_{f}\right]+\int_{t_{o}}^{t_{f}} L[x(t), u(t), t] d t \tag 1
x˙=f[x(t),u(t),t];x(to)=x0to≤t≤tfu(t)minJ=φ[x(tf),tf]+∫totfL[x(t),u(t),t]dt(1)
式中包括了控制项
u
(
t
)
u(t)
u(t)。这样的问题仍按照上一章的方法来考虑，对动力学方程约束引入Lagrange乘子
J
[
x
(
t
)
,
x
˙
(
t
)
,
u
(
t
)
,
t
]
=
φ
(
0
)
+
∫
t
t
f
{
L
+
λ
T
[
f
(
x
,
u
,
t
)
−
x
˙
]
+
d
φ
(
x
,
t
)
d
t
}
d
t
=
φ
(
0
)
+
∫
t
t
f
H
ˉ
(
x
,
x
˙
,
λ
,
u
,
t
)
d
t
(†)
J[x(t),\dot x(t),u(t),t]=\varphi(0)+\int_{t}^{t_{f}}\left\{L+\lambda^{\mathrm T}[f(x, u, t)-\dot{x}]+\frac{\text d\varphi(x, t)}{\text d t}\right\} d t\\ =\varphi(0)+\int_{t}^{t_{f}}\bar H(x,\dot x,\lambda,u,t)\text d t\tag{\dag}
J[x(t),x˙(t),u(t),t]=φ(0)+∫ttf{L+λT[f(x,u,t)−x˙]+dtdφ(x,t)}dt=φ(0)+∫ttfHˉ(x,x˙,λ,u,t)dt(†)
对 x ( t ) x(t) x(t)、 u ( t ) u(t) u(t)和 λ ( t ) \lambda(t) λ(t)都考虑最优性必要条件，即 H ˉ x − d d t H ˉ x ˙ = 0 \bar H_x-\frac{\text d}{\text d t}\bar H_{\dot x}=0 Hˉx−dtdHˉx˙=0以及 H ˉ u = 0 \bar H_u=0 Hˉu=0， H ˉ λ = 0 \bar H_\lambda=0 Hˉλ=0，可以得到Euler方程：
∂ L ∂ x + ∂ f T ∂ x λ ( t ) + λ ˙ ( t ) = 0 f ( x , u , t ) − x ˙ = 0 ∂ L ∂ u + ∂ f T ∂ u λ λ ( t ) = 0 \begin{equation}\begin{aligned} \frac{\partial L}{\partial x}+\frac{\partial f^{\mathrm{T}}}{\partial x} \lambda(t)+\dot{\lambda}(t)=0 \\ f(x,u,t)-\dot x=0\\ \frac{\partial L}{\partial u}+\frac{\partial f^{\mathrm{T}}}{\partial u^{\lambda}} \lambda(t)=0 \tag{EL1} \end{aligned}\end{equation} ∂x∂L+∂x∂fTλ(t)+λ˙(t)=0f(x,u,t)−x˙=0∂u∂L+∂uλ∂fTλ(t)=0(EL1)
此外还有状态方程和边界条件：终端固定 x ( t f ) = x f x(t_f)=x_f x(tf)=xf或终端自由 H ˉ x ˙ ( t f ) = 0 \bar H_{\dot x}(t_f)=0 Hˉx˙(tf)=0.
2. Hamilton函数法
上面这个方程的形式不是很好，我们重新定义一个哈密尔顿函数：
H [ x ( t ) , u ( t ) , λ ( t ) , t ] ≜ L [ x ( t ) , u ( t ) , t ] + λ T ( t ) f [ x ( t ) , u ( t ) , t ] (2) H[x(t), u(t), \lambda(t), t]\triangleq L[x(t), u(t), t]+\lambda^{\mathrm T}(t) f[x(t), u(t), t] \tag 2 H[x(t),u(t),λ(t),t]≜L[x(t),u(t),t]+λT(t)f[x(t),u(t),t](2)
那么性能指标化为：
J
[
x
(
t
)
,
x
˙
(
t
)
,
u
(
t
)
,
t
]
=
φ
(
0
)
−
λ
T
x
∣
0
t
f
+
∫
t
o
t
f
(
H
+
λ
˙
T
x
)
d
t
\begin{aligned} J[x(t),\dot x(t),u(t),t] &=\varphi(0)-\lambda^{\mathrm T} x\big|_0^{t_f} &+\int_{t_{o}}^{t_{f}}(H+\dot{\lambda}^{\mathrm T} x) d t \end{aligned}
J[x(t),x˙(t),u(t),t]=φ(0)−λTx∣
∣0tf+∫totf(H+λ˙Tx)dt
此时，仿照上面式
(EL1)
\text{(EL1)}
(EL1)，可以重写Euler方程为以下几个等式：
λ
˙
=
−
∂
H
∂
x
=
−
∂
L
∂
x
−
λ
T
∂
f
∂
x
x
˙
=
∂
H
∂
λ
=
f
(
x
,
u
,
t
)
0
=
∂
H
∂
u
=
∂
L
∂
u
+
(
∂
f
∂
u
)
T
λ
(3,EL2)
\begin{aligned} \dot{\lambda}&=-\frac{\partial H}{\partial x}=-\frac{\partial L}{\partial x}-\lambda^{\mathrm T}\frac{\partial f}{\partial x}\tag{3,EL2}\\ \dot{x}&=\frac{\partial H}{\partial\lambda}=f(x,u,t)\\ 0&=\frac{\partial H}{\partial u}=\frac{\partial L}{\partial u}+\left(\frac{\partial f}{\partial u}\right)^{\mathrm T} \lambda \end{aligned}
λ˙x˙0=−∂x∂H=−∂x∂L−λT∂x∂f=∂λ∂H=f(x,u,t)=∂u∂H=∂u∂L+(∂u∂f)Tλ(3,EL2)
这样，最优控制问题被规范化为3个Euler方程，按公式 ( 3 ) (3) (3)，依次是协态方程、状态方程和控制方程，方程按照Hamilton函数的偏导数的形式非常简洁。由于式 (EL2) \text{(EL2)} (EL2)由2个微分方程构成，实际上Hamilton函数法构造了两点边值问题，只要同时得到协态、状态，就相当于获得了最优控制。
按照公式 ( 1 ) − ( 3 ) (1)-(3) (1)−(3)的过程进行展开来求解，哈密尔顿函数法求解最优控制问题的具体过程如下：
- 首先写出哈密尔顿函数 H = L + λ T f H=L+\lambda^Tf H=L+λTf
- 依次列写协态方程 λ ˙ = − ∂ H ∂ x \dot\lambda=-\frac{\partial H}{\partial x} λ˙=−∂x∂H、控制方程 ∂ H ∂ u = 0 \frac{\partial H}{\partial u}=0 ∂u∂H=0
- 将最优控制代入状态方程 x ˙ = f ( x , u , t ) \dot x=f(x,u,t) x˙=f(x,u,t)
- 写出边界条件和横截条件如 x ( t f ) , λ ( t f ) , H ( ∗ , t f ) x(t_f),\lambda(t_f),H(*,t_f) x(tf),λ(tf),H(∗,tf)
- 求解整个Hamilton系统
这个方法在通用性很强，可以解决大多数无约束问题、以及带有终端约束的最优控制问题。
2.1 Hamilton函数的性质
2.1.1 哈密尔顿系统
哈密尔顿函数的引入，使得无约束最优控制的Euler-Lagrange方程可以用一阶常微分方程组描述。两者之间的关系如下，一阶必要条件Euler-Lagrange为：
∂
L
[
x
(
t
)
,
x
˙
(
t
)
,
t
]
∂
x
−
d
d
t
∂
L
[
x
(
t
)
,
x
˙
(
t
)
,
t
]
∂
x
˙
=
0
\frac{\partial L[x(t),\dot x(t),t]}{\partial x}-\frac{\text d}{\text d t}\frac{\partial L[x(t),\dot x (t),t]}{\partial {\dot x}}=0
∂x∂L[x(t),x˙(t),t]−dtd∂x˙∂L[x(t),x˙(t),t]=0
定义协态变量
λ
(
t
)
=
−
∂
L
[
x
(
t
)
,
x
˙
(
t
)
,
t
]
∂
x
˙
\lambda(t)=-\frac{\partial L[x(t),\dot x(t),t]}{\partial {\dot x}}
λ(t)=−∂x˙∂L[x(t),x˙(t),t]把E-L方程代入得到
λ
˙
=
−
L
x
\dot\lambda=-L_x
λ˙=−Lx. 定义哈密尔顿函数
H
[
x
(
t
)
,
λ
(
t
)
,
t
]
=
L
[
x
(
t
)
,
x
˙
(
t
)
,
t
]
+
λ
T
(
t
)
f
[
x
(
t
)
,
u
(
t
)
,
t
]
H[x(t),\lambda(t),t]=L[x(t),\dot x(t),t]+\lambda^\mathrm T (t)f[x(t),u(t),t]
H[x(t),λ(t),t]=L[x(t),x˙(t),t]+λT(t)f[x(t),u(t),t]，并求它的梯度
∂
H
∂
x
=
λ
T
∂
f
∂
x
+
∂
L
∂
x
+
∂
L
∂
x
˙
∂
x
˙
∂
x
=
∂
L
∂
x
=
−
λ
˙
\begin{aligned}\frac{\partial H}{\partial x}=&\lambda^\mathrm T\frac{\partial f}{\partial x}+\frac{\partial L}{\partial x}+\frac{\partial L}{\partial \dot x}\frac{\partial \dot x}{\partial x}=\frac{\partial L}{\partial x}=-\dot\lambda \end{aligned}
∂x∂H=λT∂x∂f+∂x∂L+∂x˙∂L∂x∂x˙=∂x∂L=−λ˙
另外，
x
˙
=
f
=
∂
H
∂
λ
\dot x=f=\frac{\partial H}{\partial \lambda}
x˙=f=∂λ∂H，于是得到以下Hamilton系统
x
˙
=
∂
H
∂
λ
λ
˙
=
−
∂
H
∂
x
\begin{aligned} \dot x=& \frac{\partial H}{\partial \lambda}\\ \dot\lambda=&-\frac{\partial H}{\partial x} \end{aligned}
x˙=λ˙=∂λ∂H−∂x∂H
以上这个结论在求解BVP问题时会用到，它的性质对于数值求解的影响是负面性的，Diehl【3】有以下描述：
2.1.2 时间不变性
沿最优轨线
x
∗
(
t
)
x^*(t)
x∗(t)，Hamilton函数对时间的全导数等于其对时间的偏导数，即
d
H
d
t
=
∂
H
∂
t
(1)
\frac{\mathrm{d} H}{\mathrm{d} t}=\frac{\partial H}{\partial t}\tag{1}
dtdH=∂t∂H(1)
证明：对Hamiltonian按照链式求导法则全导数：
d
H
d
t
=
∂
H
T
∂
x
x
˙
+
∂
H
T
∂
λ
λ
˙
+
∂
H
T
∂
u
U
˙
+
∂
H
∂
t
\frac{\mathrm{d} H}{\mathrm{d} t}=\frac{\partial H^{\mathrm{T}}}{\partial x} \dot{x}+\frac{\partial H^{\mathrm{T}}}{\partial \lambda} \dot{\lambda}+\frac{\partial H^{\mathrm{T}}}{\partial u} \dot{U}+\frac{\partial H}{\partial t}
dtdH=∂x∂HTx˙+∂λ∂HTλ˙+∂u∂HTU˙+∂t∂H
考虑到最优轨线附近满足
−
∂
H
∂
x
=
λ
˙
∂
H
∂
u
=
0
∂
H
∂
λ
=
f
=
x
˙
\begin{aligned} -\frac{\partial H}{\partial x}&=\dot\lambda\\ \frac{\partial H}{\partial u}&=0\\ \frac{\partial H}{\partial \lambda}&=f=\dot x \end{aligned}
−∂x∂H∂u∂H∂λ∂H=λ˙=0=f=x˙
代入则公式
(
4
)
(4)
(4)可证。
□
\square
□
此外，若Hamiltonian不显含时间t，则显然有
d
H
d
t
=
∂
H
∂
t
=
0
\frac{\mathrm{d} H}{\mathrm{d} t}=\frac{\partial H}{\partial t}=0
dtdH=∂t∂H=0
于是可得 H ( x ∗ ( t ) , u ∗ ( t ) , λ ∗ ( t ) , t ) = Const H(x^*(t),u^*(t),\lambda^*(t),t)=\text{Const} H(x∗(t),u∗(t),λ∗(t),t)=Const，若进一步考虑边界条件，对于 t f t_f tf固定的，则 H ( ∗ , t ) = H ( 0 ) H(*,t)=H(0) H(∗,t)=H(0)； 对于 t f t_f tf自由的问题，查表，有终端约束 H ( ∗ , t f ) = 0 H(*,t_f)=0 H(∗,tf)=0，则 H ( ∗ , t ) = 0 H(*,t)=0 H(∗,t)=0，也就是说终端时间自由的最优控制必然有哈密尔顿函数始终为0.
2.2 Hamilton函数的边界条件和横截条件
除了Euler方程，还要考虑边界条件和定解条件才能实际求解。
方程中
x
(
t
)
,
λ
(
t
)
∈
R
n
,
u
(
t
)
∈
R
m
x(t),\lambda(t)\in \Reals^n,u(t)\in\Reals^m
x(t),λ(t)∈Rn,u(t)∈Rm总共有
2
n
+
m
2n+m
2n+m个未知的时变参数。协态方程和状态方程
x
(
t
)
,
λ
(
t
)
x(t),\lambda(t)
x(t),λ(t)是一阶常微分方程组，需要知道
2
n
2n
2n个边界条件才能求解；控制方程
u
(
t
)
u(t)
u(t)是代数方程，由
x
(
t
)
x(t)
x(t)和
λ
(
t
)
\lambda(t)
λ(t)直接得到。
下面给出几种常用的边界条件和横截条件，式中变量均按照问题
(
1
)
(1)
(1)的表述。
| 问题描述 | 未知变量个数 | 边界条件 | 横截条件 |
|---|---|---|---|
| t f , x f t_f,x_f tf,xf均给定 | 2 n 2n 2n | x ( t 0 ) = x 0 , x ( t f ) = x f x(t_0)=x_0,x(t_f)=x_f x(t0)=x0,x(tf)=xf | \ |
| t f t_f tf给定， x f x_f xf自由 | 2 n 2n 2n | x ( t 0 ) = x 0 x(t_0)=x_0 x(t0)=x0 | λ ( t f ) = ∂ φ ( ⋅ ∗ , t f ) ∂ x \begin{aligned}\lambda(t_f)=\frac{\partial \varphi(\cdot^*,t_f)}{\partial x}\end{aligned} λ(tf)=∂x∂φ(⋅∗,tf) |
| t f t_f tf自由， x f x_f xf给定 | 2 n + 1 2n+1 2n+1 | x ( t 0 ) = x 0 , x ( t f ) = x f x(t_0)=x_0,x(t_f)=x_f x(t0)=x0,x(tf)=xf | H ( ⋅ ∗ , t f ) + ∂ φ ( ⋅ ∗ , t f ) ∂ t = 0 \begin{aligned}H(\cdot^*,t_f)+\frac{\partial \varphi(\cdot^*,t_f)}{\partial t}=0\end{aligned} H(⋅∗,tf)+∂t∂φ(⋅∗,tf)=0 |
| t f ， x f t_f，x_f tf，xf均自由 | 2 n + 1 2n+1 2n+1 | x ( t 0 ) = x 0 x(t_0)=x_0 x(t0)=x0 | λ ( t f ) = ∂ φ ∂ x ; H ( ⋅ ∗ , t f ) + ∂ φ ( ⋅ ∗ , t f ) ∂ t = 0 \begin{aligned}&\lambda(t_f)=\frac{\partial \varphi}{\partial x};\\&H(\cdot^*,t_f)+\frac{\partial \varphi(\cdot^*,t_f)}{\partial t}=0\end{aligned} λ(tf)=∂x∂φ;H(⋅∗,tf)+∂t∂φ(⋅∗,tf)=0 |
- 上面问题 ( 1 ) (1) (1)中的性能指标若为Meyer型，即 φ ( x ( t f ) , t f ) ) ≡ 0 \varphi(x(t_f),t_f))\equiv0 φ(x(tf),tf))≡0，则横截条件中出现相应的项为0。
- 另外，表1中的形式是简写，还要把哈密尔顿函数展开。如对于第四行
t
f
,
x
f
t_f,x_f
tf,xf均自由时，可以把横截条件代入Hamiltonian，得
λ ( t f ) = ∂ φ ∂ x H ( ⋅ ∗ , t f ) + ∂ φ ( ⋅ ∗ , t f ) ∂ t = [ L + ∂ φ ∂ x f + ∂ φ ∂ t ] t f = 0 \begin{aligned}&\lambda(t_f)=\frac{\partial \varphi}{\partial x}\\ &H(\cdot^*,t_f)+\frac{\partial \varphi(\cdot^*,t_f)}{\partial t}= \left[L+\frac{\partial \varphi}{\partial x}f+\frac{\partial \varphi}{\partial t}\right]_{t_f}=0\end{aligned} λ(tf)=∂x∂φH(⋅∗,tf)+∂t∂φ(⋅∗,tf)=[L+∂x∂φf+∂t∂φ]tf=0
3. 终端约束时的横截条件
设终端时刻
t
f
t_f
tf自由或给定，终端状态
x
f
x_f
xf自由但满足代数约束，两者之间的关系为
ψ
(
x
(
t
f
)
,
t
f
)
=
0
，
ψ
∈
R
q
,
q
<
n
(5)
\psi(x(t_f),t_f)=0，\psi\in\Reals^q,q\lt n\tag 5
ψ(x(tf),tf)=0，ψ∈Rq,q<n(5)
有q个终端约束，仍考虑表达式 ( 1 ) (1) (1)所述的性能指标。这样的终端约束可以表达以下两种关系，如：
- x f x_f xf的部分状态量 x i ( t f ) = x f ( i ) , i = 1 , 2 , … , q < n x_i(t_f)=x_{f}^{(i)},i=1,2,\dots,q<n xi(tf)=xf(i),i=1,2,…,q<n给定，其他状态量自由；
- x f x_f xf互相之间存在代数等式约束关系；
参考文献[2]，按照Lagrange乘数法，设一个常数向量
μ
∈
R
q
\mu\in\Reals^{q}
μ∈Rq对终端约束函数进行相乘，则性能指标变成
J
=
[
φ
+
μ
T
ψ
]
t
f
+
∫
0
t
f
{
L
(
x
,
u
,
t
)
+
λ
T
[
f
(
x
,
u
,
t
)
−
x
˙
]
}
d
t
=
Φ
t
f
+
∫
0
t
f
(
H
−
λ
T
x
˙
)
d
t
\begin{aligned} J&=\left[\varphi+\mu^{\mathrm T} \psi\right]_{t_{f}}+\int_{0}^{t_{f}}\left\{L(x, u, t)+\lambda^{\mathrm T}[f(x, u, t)-\dot{x}]\right\} d t\\ &=\Phi_{t_f}+\int_{0}^{t_{f}}(H-\lambda^{\mathrm T}\dot x)dt \end{aligned}
J=[φ+μTψ]tf+∫0tf{L(x,u,t)+λT[f(x,u,t)−x˙]}dt=Φtf+∫0tf(H−λTx˙)dt
上式仍然定义相同的Hamilton函数
H
≜
L
+
λ
T
f
H\triangleq L+\lambda^{\mathrm T}f
H≜L+λTf，以及一个新定义的标量函数
Φ
(
x
(
t
f
)
,
t
f
)
≜
φ
+
μ
T
ψ
(‡)
\color{blue}\Phi(\mathbf x(t_f),t_f)\triangleq \varphi+\mu^{\mathrm T} \psi\tag\ddag
Φ(x(tf),tf)≜φ+μTψ(‡)用于解决终端约束。接下来对终端时刻的性能指标求全微分：
d
J
=
(
(
∂
Φ
∂
t
+
L
)
d
t
+
∂
Φ
∂
x
d
x
)
t
f
+
∫
0
t
f
(
∂
H
∂
x
δ
x
+
∂
H
∂
u
δ
u
−
λ
T
δ
x
˙
)
d
t
\begin{aligned} d J=\left(\left(\frac{\partial \Phi}{\partial t}+L\right) d t+\frac{\partial \Phi}{\partial x} d x\right) _{t_f} &+\int_{0}^{t_{f}}\left(\frac{\partial H}{\partial x} \delta x+\frac{\partial H}{\partial u} \delta u-\lambda^{\mathrm T} \delta \dot{x}\right) d t \end{aligned}
dJ=((∂t∂Φ+L)dt+∂x∂Φdx)tf+∫0tf(∂x∂Hδx+∂u∂Hδu−λTδx˙)dt
并考虑
δ
x
(
t
)
=
d
x
(
t
)
−
x
˙
(
t
)
d
t
\delta x(t)=\text d x(t)-\dot x(t)\text d t
δx(t)=dx(t)−x˙(t)dt，上式可变换为
d
J
=
(
∂
Φ
∂
t
+
L
+
λ
T
x
˙
)
t
f
d
t
f
+
[
(
∂
Φ
∂
x
−
λ
T
)
d
x
]
t
f
+
(
λ
T
δ
x
)
t
0
+
∫
0
t
f
[
(
∂
H
∂
x
+
λ
˙
T
)
δ
x
+
∂
H
∂
u
δ
u
]
d
t
(6)
\text d J=\left(\frac{\partial \Phi}{\partial t}+L+\lambda^{\mathrm T} \dot{x}\right)_{t_{f}}\text d t_{f}+\left[\left(\frac{\partial \Phi}{\partial x}-\lambda^{\mathrm T}\right)\text d x\right]_{t_{f}}+\\\left(\lambda^{\mathrm T} \delta x\right)_{t_0} +\int_{0}^{t_{f}}\left[\left(\frac{\partial H}{\partial x}+\dot{\lambda}^{\mathrm T}\right) \delta x+\frac{\partial H}{\partial u} \delta u\right]\text d t \tag 6
dJ=(∂t∂Φ+L+λTx˙)tfdtf+[(∂x∂Φ−λT)dx]tf+(λTδx)t0+∫0tf[(∂x∂H+λ˙T)δx+∂u∂Hδu]dt(6)
针对上公式的结果进行分析。
3.1 性能指标变分的推导结果
公式
(
6
)
(6)
(6)中，按照最优性的必要条件，令每一项的系数都为0，可以得到Euler方程、
λ
˙
=
−
∂
H
∂
x
=
−
∂
L
∂
x
−
λ
T
∂
f
∂
x
∂
H
∂
u
=
0
=
∂
L
∂
u
+
(
∂
f
∂
u
)
T
λ
(7)
\begin{aligned} \dot{\lambda}&=-\frac{\partial H}{\partial x}=-\frac{\partial L}{\partial x}-\lambda^{\mathrm T}\frac{\partial f}{\partial x}\\ \frac{\partial H}{\partial u}&=0=\frac{\partial L}{\partial u}+\left(\frac{\partial f}{\partial u}\right)^{\mathrm T} \lambda \end{aligned}\tag{7}
λ˙∂u∂H=−∂x∂H=−∂x∂L−λT∂x∂f=0=∂u∂L+(∂u∂f)Tλ(7)
协态变量的定解条件、
λ
T
(
t
f
)
=
∂
Φ
(
⋅
∗
,
t
f
)
∂
x
=
∂
φ
(
x
f
,
t
f
)
∂
x
+
μ
T
∂
ψ
(
x
f
,
t
f
)
∂
x
(8)
\begin{aligned} \lambda^{\mathrm T}\left(t_{f}\right)=\frac{\partial \Phi(\cdot^*,t_f)}{\partial x}=\frac{\partial \varphi(x_f,t_f)}{\partial x}+\mu^{\mathrm T} \frac{\partial \psi(x_f,t_f)}{\partial x} \end{aligned}\tag{8}
λT(tf)=∂x∂Φ(⋅∗,tf)=∂x∂φ(xf,tf)+μT∂x∂ψ(xf,tf)(8)
以及时间
t
f
t_f
tf不固定的定解条件：
(
∂
Φ
∂
t
+
λ
T
x
˙
+
L
)
t
=
t
f
≡
(
d
Φ
d
t
+
L
)
t
=
t
f
=
0
(9)
\left(\frac{\partial \Phi}{\partial t}+\lambda^{\mathrm T} \dot{x}+L\right)_{t=t_{f}}\equiv \left(\frac{\text d \Phi}{\text d t}+L\right)_{t=t_{f}}=0 \tag{9}
(∂t∂Φ+λTx˙+L)t=tf≡(dtdΦ+L)t=tf=0(9)
可见文献[2]按照公式 ( 6 ) (6) (6)推导的结果和变分法得到的结果完全一致。
3.2 初值部分未知的处理
公式
(
6
)
(6)
(6)中，如果初始状态完全给定，即
δ
x
(
t
0
)
=
0
\delta \mathbf{x}(t_0)=0
δx(t0)=0，则总和
λ
T
δ
x
∣
t
0
=
0
{\lambda^{\mathrm T} \delta \mathbf x}|_{t_0}=0
λTδx∣t0=0；如果初始状态部分给定，而其余状态未知，即
x
j
(
t
0
)
,
j
=
1
,
2
,
.
.
.
,
k
x_j(t_0),j=1,2,...,k
xj(t0),j=1,2,...,k，也就是任意的变分需要遵循：
δ
x
j
(
t
0
)
{
≠
0
,
j
=
1
,
2
,
.
.
.
,
k
=
0
,
j
=
k
+
1
,
.
.
.
,
n
\delta x_j(t_0) \left\{\begin{matrix} \neq0&,j=1,2,...,k\\ =0& ,j=k+1,...,n \end{matrix} \right.
δxj(t0){=0=0,j=1,2,...,k,j=k+1,...,n
要使得 λ T δ x ∣ t 0 = 0 {\lambda^{\mathrm T} \delta \mathbf{x}}|_{t_0}=0 λTδx∣t0=0，则要求对应拉格朗日乘子为0，即 λ j ( t 0 ) = 0 , j = 1 , 2 , . . . , k \lambda_j(t_0)=0,j=1,2,...,k λj(t0)=0,j=1,2,...,k。于是有这样的结论：初值未给定状态的协态变量初值为0，即 λ j ( t 0 ) = 0 , 若 x j ( t 0 ) 未知 , j = 1 , 2 , . . . k 。 (10) \lambda_j(t_0)=0,\ 若x_j(t_0)未知\tag{10}\ ,j=1,2,...k。 λj(t0)=0, 若xj(t0)未知 ,j=1,2,...k。(10)在BVP问题中， x j ( t 0 ) x_j(t_0) xj(t0)未知而 λ j ( t 0 ) \lambda_j(t_0) λj(t0)已知，构成的问题仍然可解。
3.3 表格总结
实际上，公式 ( 6 ) (6) (6)就是把所有的边界条件写进一个式子里的表达形式，虽然推导麻烦，但是很凝练。其中 x ( t ) , λ ( t ) x(t),\lambda(t) x(t),λ(t)以及Lagrange乘数 μ \mu μ未知，以下再给出表格总结：
| 问题描述 | 未知变量个数 | 边界条件 | 横截条件 |
|---|---|---|---|
| t f , x f t_f,x_f tf,xf均给定 | 2 n 2n 2n | x ( t 0 ) = x 0 ; x ( t f ) = x f x(t_0)=x_0;\ x(t_f)=x_f x(t0)=x0; x(tf)=xf | \ |
| t f t_f tf给定， x f x_f xf自由 | 2 n → ( x , λ ) 2n\ \rightarrow( x,\lambda) 2n →(x,λ) | x ( t 0 ) = x 0 x(t_0)=x_0 x(t0)=x0 | λ ( t f ) = ∂ φ ∂ x \begin{aligned}\lambda(t_f)=\frac{\partial \varphi}{\partial x}\end{aligned} λ(tf)=∂x∂φ |
| t f , x f t_f,x_f tf,xf均自由 | 2 n + 1 → ( x , λ , t f ) 2n+1\ \rightarrow( x,\lambda,t_f) 2n+1 →(x,λ,tf) | x ( t 0 ) = x 0 x(t_0)=x_0 x(t0)=x0 | λ ( t f ) = ∂ φ ∂ x ; { L + ∂ φ ∂ x f + ∂ φ ∂ t } t f = 0 \begin{aligned}&\lambda(t_f)=\frac{\partial \varphi}{\partial x};\\&\left\{L+\frac{\partial \varphi}{\partial x}f+\frac{\partial \varphi}{\partial t}\right\}_{t_f}=0\end{aligned} λ(tf)=∂x∂φ;{L+∂x∂φf+∂t∂φ}tf=0 |
| t f t_f tf给定， x f x_f xf自由，且有终端约束 ψ ( x f , t f ) = 0 \psi(x_f,t_f)=0 ψ(xf,tf)=0 | 2 n + q → ( x , λ , μ ) 2n+q\ \rightarrow( x,\lambda,\mu) 2n+q →(x,λ,μ) | x ( t 0 ) = x 0 ; ψ ( x f , t f ) = 0 x(t_0)=x_0;\\ \psi(x_f,t_f)=0 x(t0)=x0;ψ(xf,tf)=0 | λ ( t f ) = ∂ Φ ∂ x \begin{aligned}\lambda(t_f)=\frac{\partial \Phi}{\partial x}\end{aligned} λ(tf)=∂x∂Φ |
| t f ， x f t_f，x_f tf，xf均自由，且有终端约束 ψ ( x f , t f ) = 0 \psi(x_f,t_f)=0 ψ(xf,tf)=0 | 2 n + q + 1 → ( x , λ , μ , t f ) 2n+q+1\ \rightarrow( x,\lambda,\mu,t_f) 2n+q+1 →(x,λ,μ,tf) | x ( t 0 ) = x 0 ; ψ ( x f , t f ) = 0 x(t_0)=x_0;\\ \psi(x_f,t_f)=0 x(t0)=x0;ψ(xf,tf)=0 | λ ( t f ) = ∂ Φ ∂ x ; { L + ∂ Φ ∂ x f + ∂ Φ ∂ t } t f = 0 \begin{aligned}&\lambda(t_f)=\frac{\partial \Phi}{\partial x};\\ &\left\{L+\frac{\partial \Phi}{\partial x}f+\frac{\partial \Phi}{\partial t}\right\}_{t_f}=0\end{aligned} λ(tf)=∂x∂Φ;{L+∂x∂Φf+∂t∂Φ}tf=0 |
1）. 第4、5行说明
上表第4、5行因为引入了拉格朗日乘子
μ
\mu
μ而与第2、3行形式一致，显得简洁，其中的标量函数
Φ
(
x
(
t
f
)
,
t
f
)
\Phi(x(t_f),t_f)
Φ(x(tf),tf)由公式
(
‡
)
(\ddag)
(‡)定义。终端状态未知所产生的横截条件展开为：
λ
(
t
f
)
=
∂
Φ
∂
x
∣
t
=
t
f
≡
∂
φ
∂
x
+
μ
T
∂
ψ
∂
x
(11)
\lambda(t_f)=\frac{\partial \Phi}{\partial x}\Big|_{t=t_f}\equiv\frac{\partial \varphi}{\partial x}+\mu^{\mathrm T}\frac{\partial\psi}{\partial x}\tag{11}
λ(tf)=∂x∂Φ∣
∣t=tf≡∂x∂φ+μT∂x∂ψ(11)
第5行比第4行多了一个条件，这个条件用于确定自由的
t
f
t_f
tf。定解条件展开后为：
[
L
+
(
∂
φ
∂
t
+
μ
T
∂
ψ
∂
t
)
+
(
∂
φ
∂
x
+
μ
T
∂
ψ
∂
x
)
f
]
t
f
=
0
(12)
\left[ L+\left(\frac{\partial \varphi}{\partial t}+\mu^{\mathrm T} \frac{\partial \psi}{\partial t}\right)+\left(\frac{\partial \varphi}{\partial x}+\mu^{\mathrm T} \frac{\partial \psi}{\partial x}\right) f\right]_{t_f}=0 \tag{12}
[L+(∂t∂φ+μT∂t∂ψ)+(∂x∂φ+μT∂x∂ψ)f]tf=0(12)
事实上，上式
(
12
)
(12)
(12)可以这样记忆。
H
∗
(
⋅
∗
,
t
f
)
+
∂
Φ
(
⋅
∗
,
t
f
)
∂
t
=
0
Φ
(
x
(
t
f
)
,
t
f
)
≡
φ
+
μ
T
ψ
H
∗
[
x
(
t
)
,
λ
(
t
)
,
x
(
t
f
)
,
λ
(
t
f
)
,
t
,
t
f
]
≡
L
+
λ
T
f
+
μ
T
ψ
\begin{aligned} &H^*(\cdot^*,t_f)+\frac{\partial \Phi(\cdot^*,t_f)}{\partial t}=0\\ &\Phi(\mathbf x(t_f),t_f)\equiv \varphi+\mu^{\mathrm T} \psi\\ &H^*[x(t),\lambda(t),x(t_f),\lambda(t_f),t,t_f]\equiv L+\lambda^\mathrm Tf+\mu^\mathrm T\psi \end{aligned}
H∗(⋅∗,tf)+∂t∂Φ(⋅∗,tf)=0Φ(x(tf),tf)≡φ+μTψH∗[x(t),λ(t),x(tf),λ(tf),t,tf]≡L+λTf+μTψ
2）. 其他型性能指标
表2中的内容假设性能指标为Bolza型（即包含终端项与积分项），但是如果性能指标只包含积分项（Lagrange型）或终端函数项（Meyer型），只需要把不存在的函数项设为0即可。
3）. 终端状态部分给定的处理
另外，作为终端约束
(
5
)
(5)
(5)的一种特殊形式，如果部分终端状态已知，即
ψ
[
x
(
t
f
)
,
t
f
]
=
x
i
(
t
f
)
−
x
i
f
=
0
，
i
=
1
,
2
,
.
.
.
,
q
,
q
<
n
\psi[x(t_f),t_f]=x_i(t_f)-x_{if}=0，i=1,2,...,q,q\lt n
ψ[x(tf),tf]=xi(tf)−xif=0，i=1,2,...,q,q<n
这种情况仍然可以套用表格中的边界条件与横截条件
x
i
(
t
f
)
=
x
i
f
λ
(
t
f
)
=
{
0
+
μ
i
,
i
=
1
,
2
,
.
.
.
q
∂
φ
∂
x
i
,
i
=
q
+
1
,
.
.
.
,
n
\begin{aligned} x_i(t_f)&=x_{if}\\ \lambda(t_f)&=\left\{ \begin{matrix} 0+\mu_i&\ ,i=1,2,...q\\ \frac{\partial \varphi}{\partial x_i}&\ ,i=q+1,...,n \end{matrix}\right. \end{aligned}
xi(tf)λ(tf)=xif={0+μi∂xi∂φ ,i=1,2,...q ,i=q+1,...,n
4. 应用举例
4.1 倒立摆问题
倒立摆按照方程
I
θ
¨
+
b
θ
˙
−
m
g
l
sin
θ
=
u
I\ddot\theta+b\dot\theta-mgl\sin\theta=u
Iθ¨+bθ˙−mglsinθ=u,初始状态
x
0
=
[
θ
,
ω
]
T
=
[
π
,
0
]
x_0=[\theta,\omega]^{\mathrm T}=[\pi,0]
x0=[θ,ω]T=[π,0]，控制目标
[
θ
f
,
ω
f
]
T
=
[
0
,
0
]
[\theta_f,\omega_f]^{\mathrm T}=[0,0]
[θf,ωf]T=[0,0]，终端时刻
t
f
t_f
tf自由，二次型性能指标
min
u
(
t
)
=
1
2
x
f
T
Q
x
f
+
1
2
∫
0
t
f
R
u
2
\min_{u(t)}=\frac1 2\mathbf x_f^{\mathrm T}\text Q\mathbf x_f+\frac1 2\int_0^{t_f}\text R\mathbf u^2
u(t)min=21xfTQxf+21∫0tfRu2
首先写出一阶非线性微分方程组
[
θ
˙
ω
˙
]
=
f
(
x
)
=
[
ω
1
/
I
(
−
m
g
l
sin
θ
+
b
ω
+
u
)
]
\begin{bmatrix}\dot\theta\\\dot\omega\end{bmatrix}=\mathbf f(\mathbf x)=\begin{bmatrix}\omega\\1/I(-mgl\sin\theta+b\omega+u)\end{bmatrix}
[θ˙ω˙]=f(x)=[ω1/I(−mglsinθ+bω+u)]
写出Hamilton函数
H
=
1
2
R
u
2
+
λ
1
ω
+
λ
2
ω
˙
H=\frac1 2\text R u^2+\lambda_1\omega+\lambda_2\dot\omega
H=21Ru2+λ1ω+λ2ω˙
协态方程
λ
˙
1
=
−
∂
H
∂
θ
=
λ
2
m
g
l
cos
θ
/
I
λ
˙
2
=
−
∂
H
∂
θ
=
−
λ
1
\begin{aligned} \dot\lambda_1&=-\frac{\partial H}{\partial \theta}=\lambda_2mgl\cos\theta/I\\ \dot\lambda_2&=-\frac{\partial H}{\partial \theta}=-\lambda_1 \end{aligned}
λ˙1λ˙2=−∂θ∂H=λ2mglcosθ/I=−∂θ∂H=−λ1
最优控制
∂
H
∂
u
=
R
u
+
λ
2
/
I
=
0
⟹
u
=
−
λ
2
/
R
I
\frac{\partial{H}}{\partial u}=Ru+\lambda_2/I=0\implies u=-\lambda_2/{RI}
∂u∂H=Ru+λ2/I=0⟹u=−λ2/RI
横截条件
H
(
t
f
)
+
∂
φ
(
x
f
)
∂
t
=
1
2
R
u
2
+
λ
2
u
/
I
=
0
⟹
λ
2
(
t
f
)
=
0
H(t_f)+\frac{\partial \varphi(x_f)}{\partial t}=\frac1 2\text R u^2+\lambda_2u/I=0\implies \lambda_2(t_f)=0
H(tf)+∂t∂φ(xf)=21Ru2+λ2u/I=0⟹λ2(tf)=0
代入数据，调用MATLAB中的 sol = bvp4c(odefun,bcfun,solinit,options) \texttt{sol = bvp4c(odefun,bcfun,solinit,options)} sol = bvp4c(odefun,bcfun,solinit,options)求解这个问题，得到结果。
4.2 连续推力轨道转移问题
轨道动力学方程
r
¨
=
−
μ
r
3
r
+
a
\mathbf{\ddot r}=\mathbf -\frac\mu{r^3}\mathbf r+\mathbf a
r¨=−r3μr+a，初始状态已知
r
(
t
0
)
=
r
0
,
v
(
t
0
)
=
v
0
r(t_0)=r_0,v(t_0)=v_0
r(t0)=r0,v(t0)=v0，终端时刻
t
f
t_f
tf给定，终端状态约束
ψ
(
r
(
t
f
)
,
v
(
t
f
)
)
=
r
T
v
=
0
\psi(\mathbf r(t_f),\mathbf v(t_f))=\mathbf r^{\mathrm T}\mathbf v=0
ψ(r(tf),v(tf))=rTv=0。最小能量问题
min
a
(
t
)
J
=
1
2
∫
t
0
t
f
a
T
a
d
t
\min_{a(t)}J=\frac1 2\int_{t_0}^{t_f}\mathbf a^{\mathrm T}\mathbf a\text d t
mina(t)J=21∫t0tfaTadt。
套用Hamilton函数法，状态变量
x
=
[
r
v
]
f
(
x
)
=
[
r
−
μ
r
3
r
+
a
]
\mathbf x=\begin{bmatrix}\mathbf r\\ \mathbf v\end{bmatrix}\\ \mathbf f(\mathbf x)=\begin{bmatrix}\mathbf r\\ -\frac\mu{r^3}\mathbf r+\mathbf a\end{bmatrix}
x=[rv]f(x)=[r−r3μr+a]
则Hamilton函数为
H
=
1
2
a
T
a
+
λ
r
T
v
+
λ
v
T
(
g
(
r
)
+
a
)
H=\frac 1 2\mathbf a^{\mathrm T}\mathbf a+\mathbf\lambda_r^{\mathrm T}\mathbf v+\mathbf\lambda_v^{\mathrm T}(\mathbf g(\mathbf r)+\mathbf a)
H=21aTa+λrTv+λvT(g(r)+a)
协态方程
λ
˙
r
T
=
−
∂
H
∂
r
=
−
λ
v
T
∂
g
(
r
)
∂
r
λ
˙
v
T
=
−
∂
H
∂
v
=
−
λ
r
T
\begin{aligned} \dot{\lambda}_{r}^{\mathrm T}&=-\frac{\partial H}{\partial \mathbf{r}}=-\lambda_{\mathrm{v}}^{\mathrm T} \frac{\partial \mathbf{g}(\mathbf{r})}{\partial \mathbf{r}}\\ \dot{\lambda}_{\mathrm{v}}^{\mathrm T}&=-\frac{\partial H}{\partial \mathbf{v}}=-\lambda_{r}^{\mathrm T} \end{aligned}
λ˙rTλ˙vT=−∂r∂H=−λvT∂r∂g(r)=−∂v∂H=−λrT
对终端约束引入Lagrange乘数
μ
∈
R
1
\mu\in\Reals^1
μ∈R1，查表得到横截条件
λ
r
(
t
f
)
=
μ
T
∂
ψ
∂
r
(
t
f
)
=
μ
v
f
λ
v
(
t
f
)
=
μ
T
∂
ψ
∂
v
(
t
f
)
=
μ
r
f
\lambda_{r}\left(t_{f}\right)=\mu^{\mathrm T}\frac{\partial\psi}{\partial \mathbf{r}\left(t_{f}\right)}=\mu\mathbf v_f \\ \lambda_{\mathrm{v}}\left(t_{f}\right)=\mu^{\mathrm T}\frac{\partial\psi}{\partial \mathbf{v}\left(t_{f}\right)}=\mu\mathbf r_f
λr(tf)=μT∂r(tf)∂ψ=μvfλv(tf)=μT∂v(tf)∂ψ=μrf
最优控制
∂
H
∂
a
=
a
+
λ
v
=
0
(13)
\frac{\partial{H}}{\partial \mathbf a}=\mathbf a+\lambda_v=0 \tag {13}
∂a∂H=a+λv=0(13)
则最优控制的控制律为
a
=
−
λ
v
\mathbf a=-\lambda_{\mathbf v}
a=−λv，公式
(
13
)
(13)
(13)最早由Lawden提出，被称为主矢量理论，它的主要含义是：最优推力方向总是与协态变量
λ
v
\lambda_{\mathbf v}
λv反向，只要求出协态变量就可以确定最优推力。代入最优控制的结果，构成以下12个未知变量的哈密尔顿系统，
r
˙
=
r
v
˙
=
−
μ
r
3
r
−
λ
v
λ
˙
r
T
=
−
λ
v
T
∂
g
(
r
)
∂
r
λ
˙
v
T
=
−
λ
r
T
\begin{aligned} \dot{\mathbf r}&=\mathbf r\\ \dot{\mathbf v}&= -\frac\mu{r^3}\mathbf {r}-\lambda_v\\ \dot{\lambda}_{r}^{\mathrm T}&=-\lambda_{\mathrm{v}}^{\mathrm T} \frac{\partial \mathbf{g}(\mathbf{r})}{\partial \mathbf{r}}\\ \dot{\lambda}_{\mathrm{v}}^{\mathrm T}&=-\lambda_{r}^{\mathrm T} \end{aligned}
r˙v˙λ˙rTλ˙vT=r=−r3μr−λv=−λvT∂r∂g(r)=−λrT
1个未知常数 μ \mu μ，边界条件共有13个，可以通过求解两点边值问题求解最优轨迹。这个问题的主要难点是确定协态变量初值，只要得到它就可以用数值积分方法求解，然而难点在于如何满足终端约束。
参考文献
[1] 邢继祥. 最优控制应用基础[M]. 科学出版社, 2003.
[2] Bryson A E , Ho Y C ,Applied optimal control : optimization, estimation, and control[J]. IEEE Transactions on Systems Man & Cybernetics, 1975
[3] Moritz Diehl, Numerical Optimal Control (draft), 2011
[4] 还有一些不重要的内容被我放到另一篇博客里了： 最优控制理论 二+、哈密尔顿函数法的补充
462
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
