最优控制理论二、哈密尔顿函数法

Hamilton函数方法是变分法应用在控制系统上的标准化方法，即使不懂变分法，简单套用表格中的公式也可以列写出方程，这个方法是最优控制理论用的最多的方法。
本篇博客是本系列最长也是最重要的一篇，持续更新，欢迎同学和朋友们提出修改建议。

Hamiltonian 目录

1. 规范化的最优控制问题
2. Hamilton函数法
- 2.1 Hamilton函数的性质
- - 2.1.1 哈密尔顿系统
  - 2.1.2 时间不变性
- 2.2 Hamilton函数的边界条件和横截条件
3. 终端约束时的横截条件
- 3.1 性能指标变分的推导结果
- 3.2 初值部分未知的处理
- 3.3 表格总结
- - 1）. 第4、5行说明
  - 2）. 其他型性能指标
  - 3）. 终端状态部分给定的处理
4. 应用举例
- 4.1 倒立摆问题
- 4.2 连续推力轨道转移问题
参考文献

1. 规范化的最优控制问题

按照第一章最优控制理论一、变分法和泛函极值问题，我们已经讨论了有动力学方程约束f(x,x˙,t)=0\boldsymbol f(\boldsymbol x,\dot {\boldsymbol x},t)=0f(x,x˙,t)=0的动态系统，若无其他约束，这个系统的最优轨线遵循以下必要条件
Hx−ddtHx˙=0f(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)−λTf(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)。控制系统中更常见的一阶非线性系统方程，问题是这样的：tft_ftf给定，终端状态未知或已知（仅边界条件不同），除状态方程外没有约束，且
x˙=f[x(t),u(t),t];x(to)=x0to≤t≤tfmin⁡u(t)J=φ[x(tf),tf]+∫totfL[x(t),u(t),t]dt(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)+∫ttf{L+λT[f(x,u,t)−x˙]+dφ(x,t)dt}dt=φ(0)+∫ttfHˉ(x,x˙,λ,u,t)dt(†)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−ddtHˉx˙=0\bar H_x-\frac{\text d}{\text d t}\bar H_{\dot x}=0Hˉx−dtdHˉx˙=0以及Hˉu=0\bar H_u=0Hˉu=0，Hˉλ=0\bar H_\lambda=0Hˉλ=0，可以得到Euler方程：

∂L∂x+∂fT∂xλ(t)+λ˙(t)=0f(x,u,t)−x˙=0∂L∂u+∂fT∂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(tf)=xfx(t_f)=x_fx(tf)=xf或终端自由Hˉx˙(tf)=0\bar H_{\dot x}(t_f)=0Hˉ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 2H[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)−λTx∣0tf+∫totf(H+λ˙Tx)dt\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∂xx˙=∂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+λTfH=L+\lambda^TfH=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(tf),λ(tf),H(∗,tf)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−ddt∂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方程代入得到λ˙=−Lx\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函数对时间的全导数等于其对时间的偏导数，即
dHdt=∂H∂t(1)\frac{\mathrm{d} H}{\mathrm{d} t}=\frac{\partial H}{\partial t}\tag{1} dtdH=∂t∂H(1)

证明：对Hamiltonian按照链式求导法则全导数：
dHdt=∂HT∂xx˙+∂HT∂λλ˙+∂HT∂uU˙+∂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，则显然有
dHdt=∂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)=ConstH(x^*(t),u^*(t),\lambda^*(t),t)=\text{Const}H(x∗(t),u∗(t),λ∗(t),t)=Const，若进一步考虑边界条件，对于tft_ftf固定的，则H(∗,t)=H(0)H(*,t)=H(0)H(∗,t)=H(0)；对于tft_ftf自由的问题，查表，有终端约束H(∗,tf)=0H(*,t_f)=0H(∗,tf)=0，则H(∗,t)=0H(*,t)=0H(∗,t)=0，也就是说终端时间自由的最优控制必然有哈密尔顿函数始终为0.

2.2 Hamilton函数的边界条件和横截条件

除了Euler方程，还要考虑边界条件和定解条件才能实际求解。
方程中x(t),λ(t)∈Rn,u(t)∈Rmx(t),\lambda(t)\in \Reals^n,u(t)\in\Reals^mx(t),λ(t)∈Rn,u(t)∈Rm总共有2n+m2n+m2n+m个未知的时变参数。协态方程和状态方程x(t),λ(t)x(t),\lambda(t)x(t),λ(t)是一阶常微分方程组，需要知道2n2n2n个边界条件才能求解；控制方程u(t)u(t)u(t)是代数方程，由x(t)x(t)x(t)和λ(t)\lambda(t)λ(t)直接得到。
下面给出几种常用的边界条件和横截条件，式中变量均按照问题(1)(1)(1)的表述。

表1. 不含终端约束时的定解条件

问题描述	未知变量个数	边界条件	横截条件
tf,xft_f,x_ftf,xf均给定	2n2n2n	x(t0)=x0,x(tf)=xfx(t_0)=x_0,x(t_f)=x_fx(t0)=x0,x(tf)=xf	\
tft_ftf给定，xfx_fxf自由	2n2n2n	x(t0)=x0x(t_0)=x_0x(t0)=x0	λ(tf)=∂φ(⋅∗,tf)∂x\begin{aligned}\lambda(t_f)=\frac{\partial \varphi(\cdot^*,t_f)}{\partial x}\end{aligned}λ(tf)=∂x∂φ(⋅∗,tf)
tft_ftf自由，xfx_fxf给定	2n+12n+12n+1	x(t0)=x0,x(tf)=xfx(t_0)=x_0,x(t_f)=x_fx(t0)=x0,x(tf)=xf	H(⋅∗,tf)+∂φ(⋅∗,tf)∂t=0\begin{aligned}H(\cdot^,t_f)+\frac{\partial \varphi(\cdot^,t_f)}{\partial t}=0\end{aligned}H(⋅∗,tf)+∂t∂φ(⋅∗,tf)=0
tf，xft_f，x_ftf，xf均自由	2n+12n+12n+1	x(t0)=x0x(t_0)=x_0x(t0)=x0	λ(tf)=∂φ∂x;H(⋅∗,tf)+∂φ(⋅∗,tf)∂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(tf),tf))≡0\varphi(x(t_f),t_f))\equiv0φ(x(tf),tf))≡0，则横截条件中出现相应的项为0。
另外，表1中的形式是简写，还要把哈密尔顿函数展开。如对于第四行tf,xft_f,x_ftf,xf均自由时，可以把横截条件代入Hamiltonian，得
λ(tf)=∂φ∂xH(⋅∗,tf)+∂φ(⋅∗,tf)∂t=[L+∂φ∂xf+∂φ∂t]tf=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. 终端约束时的横截条件

设终端时刻tft_ftf自由或给定，终端状态xfx_fxf自由但满足代数约束，两者之间的关系为
ψ(x(tf),tf)=0，ψ∈Rq,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)所述的性能指标。这样的终端约束可以表达以下两种关系，如：

xfx_fxf的部分状态量xi(tf)=xf(i),i=1,2,…,q<nx_i(t_f)=x_{f}^{(i)},i=1,2,\dots,q<nxi(tf)=xf(i),i=1,2,…,q<n给定，其他状态量自由；
xfx_fxf互相之间存在代数等式约束关系；

参考文献[2]，按照Lagrange乘数法，设一个常数向量μ∈Rq\mu\in\Reals^{q}μ∈Rq对终端约束函数进行相乘，则性能指标变成
J=[φ+μTψ]tf+∫0tf{L(x,u,t)+λT[f(x,u,t)−x˙]}dt=Φtf+∫0tf(H−λTx˙)dt\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+λTfH\triangleq L+\lambda^{\mathrm T}fH≜L+λTf，以及一个新定义的标量函数Φ(x(tf),tf)≜φ+μTψ(‡)\color{blue}\Phi(\mathbf x(t_f),t_f)\triangleq \varphi+\mu^{\mathrm T} \psi\tag\ddagΦ(x(tf),tf)≜φ+μTψ(‡)用于解决终端约束。接下来对终端时刻的性能指标求全微分：
dJ=((∂Φ∂t+L)dt+∂Φ∂xdx)tf+∫0tf(∂H∂xδx+∂H∂uδu−λTδx˙)dt\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)=dx(t)−x˙(t)dt\delta x(t)=\text d x(t)-\dot x(t)\text d tδx(t)=dx(t)−x˙(t)dt，上式可变换为
dJ=(∂Φ∂t+L+λTx˙)tfdtf+[(∂Φ∂x−λT)dx]tf+(λTδx)t0+∫0tf[(∂H∂x+λ˙T)δx+∂H∂uδu]dt(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(tf)=∂Φ(⋅∗,tf)∂x=∂φ(xf,tf)∂x+μT∂ψ(xf,tf)∂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)

以及时间tft_ftf不固定的定解条件：
(∂Φ∂t+λTx˙+L)t=tf≡(dΦdt+L)t=tf=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(t0)=0\delta \mathbf{x}(t_0)=0δx(t0)=0，则总和λTδx∣t0=0{\lambda^{\mathrm T} \delta \mathbf x}|_{t_0}=0λTδx∣t0=0；如果初始状态部分给定，而其余状态未知，即xj(t0),j=1,2,...,kx_j(t_0),j=1,2,...,kxj(t0),j=1,2,...,k，也就是任意的变分需要遵循：
δxj(t0){≠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∣t0=0{\lambda^{\mathrm T} \delta \mathbf{x}}|_{t_0}=0λTδx∣t0=0，则要求对应拉格朗日乘子为0，即λj(t0)=0,j=1,2,...,k\lambda_j(t_0)=0,j=1,2,...,kλj(t0)=0,j=1,2,...,k。于是有这样的结论：初值未给定状态的协态变量初值为0，即λj(t0)=0,若xj(t0)未知 ,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问题中，xj(t0)x_j(t_0)xj(t0)未知而λj(t0)\lambda_j(t_0)λj(t0)已知，构成的问题仍然可解。

3.3 表格总结

实际上，公式(6)(6)(6)就是把所有的边界条件写进一个式子里的表达形式，虽然推导麻烦，但是很凝练。其中x(t),λ(t)x(t),\lambda(t)x(t),λ(t)以及Lagrange乘数μ\muμ未知，以下再给出表格总结：

表2. 不同问题的定解条件

问题描述	未知变量个数	边界条件	横截条件
tf,xft_f,x_ftf,xf均给定	2n2n2n	x(t0)=x0;x(tf)=xfx(t_0)=x_0;\ x(t_f)=x_fx(t0)=x0; x(tf)=xf	\
tft_ftf给定，xfx_fxf自由	2n→(x,λ)2n\ \rightarrow( x,\lambda)2n →(x,λ)	x(t0)=x0x(t_0)=x_0x(t0)=x0	λ(tf)=∂φ∂x\begin{aligned}\lambda(t_f)=\frac{\partial \varphi}{\partial x}\end{aligned}λ(tf)=∂x∂φ
tf,xft_f,x_ftf,xf均自由	2n+1→(x,λ,tf)2n+1\ \rightarrow( x,\lambda,t_f)2n+1 →(x,λ,tf)	x(t0)=x0x(t_0)=x_0x(t0)=x0	λ(tf)=∂φ∂x;{L+∂φ∂xf+∂φ∂t}tf=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
tft_ftf给定，xfx_fxf自由，且有终端约束ψ(xf,tf)=0\psi(x_f,t_f)=0ψ(xf,tf)=0	2n+q→(x,λ,μ)2n+q\ \rightarrow( x,\lambda,\mu)2n+q →(x,λ,μ)	x(t0)=x0;ψ(xf,tf)=0x(t_0)=x_0;\\ \psi(x_f,t_f)=0x(t0)=x0;ψ(xf,tf)=0	λ(tf)=∂Φ∂x\begin{aligned}\lambda(t_f)=\frac{\partial \Phi}{\partial x}\end{aligned}λ(tf)=∂x∂Φ
tf，xft_f，x_ftf，xf均自由，且有终端约束ψ(xf,tf)=0\psi(x_f,t_f)=0ψ(xf,tf)=0	2n+q+1→(x,λ,μ,tf)2n+q+1\ \rightarrow( x,\lambda,\mu,t_f)2n+q+1 →(x,λ,μ,tf)	x(t0)=x0;ψ(xf,tf)=0x(t_0)=x_0;\\ \psi(x_f,t_f)=0x(t0)=x0;ψ(xf,tf)=0	λ(tf)=∂Φ∂x;{L+∂Φ∂xf+∂Φ∂t}tf=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(tf),tf)\Phi(x(t_f),t_f)Φ(x(tf),tf)由公式(‡)(\ddag)(‡)定义。终端状态未知所产生的横截条件展开为：
λ(tf)=∂Φ∂x∣t=tf≡∂φ∂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行多了一个条件，这个条件用于确定自由的tft_ftf。定解条件展开后为：
[L+(∂φ∂t+μT∂ψ∂t)+(∂φ∂x+μT∂ψ∂x)f]tf=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∗(⋅∗,tf)+∂Φ(⋅∗,tf)∂t=0Φ(x(tf),tf)≡φ+μTψH∗[x(t),λ(t),x(tf),λ(tf),t,tf]≡L+λTf+μ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(tf),tf]=xi(tf)−xif=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

这种情况仍然可以套用表格中的边界条件与横截条件
xi(tf)=xifλ(tf)={0+μi,i=1,2,...q∂φ∂xi,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θ˙−mglsin⁡θ=uI\ddot\theta+b\dot\theta-mgl\sin\theta=uIθ¨+bθ˙−mglsinθ=u,初始状态x0=[θ,ω]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]，终端时刻tft_ftf自由，二次型性能指标
min⁡u(t)=12xfTQxf+12∫0tfRu2\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^2u(t)min=21xfTQxf+21∫0tfRu2

首先写出一阶非线性微分方程组
[θ˙ω˙]=f(x)=[ω1/I(−mglsin⁡θ+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=12Ru2+λ1ω+λ2ω˙H=\frac1 2\text R u^2+\lambda_1\omega+\lambda_2\dot\omegaH=21Ru2+λ1ω+λ2ω˙

协态方程
λ˙1=−∂H∂θ=λ2mglcos⁡θ/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=Ru+λ2/I=0⟹u=−λ2/RI\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(tf)+∂φ(xf)∂t=12Ru2+λ2u/I=0⟹λ2(tf)=0H(t_f)+\frac{\partial \varphi(x_f)}{\partial t}=\frac1 2\text R u^2+\lambda_2u/I=0\implies \lambda_2(t_f)=0H(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¨=−μr3r+a\mathbf{\ddot r}=\mathbf -\frac\mu{r^3}\mathbf r+\mathbf ar¨=−r3μr+a，初始状态已知r(t0)=r0,v(t0)=v0r(t_0)=r_0,v(t_0)=v_0r(t0)=r0,v(t0)=v0，终端时刻tft_ftf给定，终端状态约束ψ(r(tf),v(tf))=rTv=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=12∫t0tfaTadt\min_{a(t)}J=\frac1 2\int_{t_0}^{t_f}\mathbf a^{\mathrm T}\mathbf a\text d tmina(t)J=21∫t0tfaTadt。
套用Hamilton函数法，状态变量
x=[rv]f(x)=[r−μr3r+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=12aTa+λrTv+λvT(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)

协态方程
λ˙rT=−∂H∂r=−λvT∂g(r)∂rλ˙vT=−∂H∂v=−λrT\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乘数μ∈R1\mu\in\Reals^1μ∈R1，查表得到横截条件
λr(tf)=μT∂ψ∂r(tf)=μvfλv(tf)=μT∂ψ∂v(tf)=μrf\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˙=rv˙=−μr3r−λvλ˙rT=−λvT∂g(r)∂rλ˙vT=−λrT\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] 还有一些不重要的内容被我放到另一篇博客里了：最优控制理论二+、哈密尔顿函数法的补充