最优化问题求解中的汉密尔顿方程是最优控制方法解决动态优化问题的一阶必要条件

汉密尔顿方程（Hamilton Equation）

理论推导

首先看需要解决的问题

max∫t0t1f(x,u)dtx˙=g(x,u)x(t0)=x0x(t1):freemax\int_{t_0}^{t_1}f(x,u)dt\\\\ \dot{x}=g(x,u)\\\\ x(t_0)=x_0\\\\ x(t_1):free

引入一个函数 λ(t)\lambda (t) 表示一个定义在 t0≤t≤t1t_0\leq t\leq t_1 上的连续的可微函数，那么对于任何满足上述约束的 x(t),u(t)x(t),u(t) 有

∫t0t1f(x,u)dt=∫t0t1[f(x,u)+λg(x,u)−λx˙]dt\int_{t_0}^{t_1}f(x,u)dt=\int_{t_0}^{t_1}[f(x,u)+\lambda g(x,u)-\lambda \dot{x}]dt

对上式最右侧右侧部分进行分部积分可以得到

−∫t0t1λx˙dt=−λ(t1)x(t1)+λ(t0)x(t0)+∫t0t1xλ˙dt-\int_{t_0}^{t_1}\lambda \dot{x}dt=-\lambda(t_1)x(t_1)+\lambda (t_0)x(t_0)+\int_{t_0}^{t_1}x\dot{\lambda}dt

将上式带回原来的方程

∫t0t1f(x,u)dt=∫t0t1[f(x,u)+λg(x,u)+λ˙x]dt−λ(t1)x(t1)+λ(t0)x(t0)\int_{t_0}^{t_1}f(x,u)dt=\int_{t_0}^{t_1}[f(x,u)+\lambda g(x,u)+\dot{\lambda} x]dt-\lambda(t_1)x(t_1)+\lambda (t_0)x(t_0)

在最优控制方法中，可以发现，只要确定了控制，系统状态自然就被确定了，这个是由函数 g(x,u)g(x,u) 确定，所以问题变为求解最优的控制 $u^\star $

引入一个扰动项函数 h(t)h(t) 是任意的但是人为给定，令

u=u⋆+ahu=u^\star +ah

并且定义函数 y(t,a)y(t,a) 为关于 u⋆+ahu^\star +ah 的状态变量，并且满足

y(t,0)=x⋆y(t0,a)=x0y(t,0)=x^\star \\\\ y(t_0,a)=x_0

之前的目标函数可以写作为

∫t0t1f(y(t,a),u⋆+ah)dt\int_{t_0}^{t_1}f(y(t,a),u^\star +ah)dt

并且得到

J(a)=∫t0t1f(y,u⋆+ah)dt=∫t0t1[f(y,u⋆+ah)+λg(y,u⋆+ah)+λ˙y]dt−λ(t1)y(t1,a)+λ(t0)y(t0,a)J(a)=\int_{t_0}^{t_1}f(y,u^\star +ah)dt=\int_{t_0}^{t_1}[f(y,u^\star +ah)+\lambda g(y,u^\star +ah)+\dot{\lambda} y]dt-\lambda(t_1)y(t_1,a)+\lambda (t_0)y(t_0,a)

由于 $u^\star $ 是最优控制变量，所以该函数在 a=0a=0 处取得最大值，也就是 J˙(0)=0\dot{J}(0)=0，也就是

J˙(0)=∫t0t1[(fx+λgx+λ˙)ya+(fu+λgu)h]dt−λ(t1)ya(t1,0)\dot{J}(0)=\int_{t_0}^{t_1}[(f_x+\lambda g_x+\dot{\lambda})y_a+(f_u+\lambda g_u)h]dt-\lambda (t_1)y_a(t_1,0)

由于 y(t0,a)y(t_0,a) 是常数，所以 ya(t0,0)=0y_a(t_0,0)=0

上式中难以确定的是 yay_a ，因为不知道 y(t,a)y(t,a) 的具体形式，但是可以通过构造 λ\lambda 来消去含有 yy 的项，定义

λ˙=−[fx(x⋆,u⋆)+λgx(x⋆,u⋆)]λ(t1)=0\dot\lambda=-[f_x(x^\star ,u^\star )+\lambda g_x(x^\star ,u^\star )]\\\lambda(t_1)=0

所以因此可以得到

∫t0t1−[fu(x⋆,u⋆)+λgu(x⋆,u⋆)]hdt=0\int_{t_0}^{t_1}-[f_u(x^\star ,u^\star )+\lambda g_u(x^\star ,u^\star )]hdt=0

并且对于任意连续的 hh 都成立，所以令

h=fu(x⋆,u⋆)+λgu(x⋆,u⋆)h=f_u(x^\star ,u^\star )+\lambda g_u(x^\star ,u^\star )

得

∫t0t1−[fu(x⋆,u⋆)+λgu(x⋆,u⋆)]2dt=0\int_{t_0}^{t_1}-[f_u(x^\star ,u^\star )+\lambda g_u(x^\star ,u^\star )]^2dt=0

即

fu(x⋆,u⋆)+λgu(x⋆,u⋆)=0  t0≤t≤t1f_u(x^\star ,u^\star )+\lambda g_u(x^\star ,u^\star )=0~~t_0\leq t\leq t_1

也就是说，如果 $u^\star ,x^\star $ 是最优解，存在一个连续函数 λ(t)\lambda(t) ，他们同时满足

State Equation:x˙=g(x,u)x(t0)=x0Multiplier Equation:λ˙=−[fx(x⋆,u⋆)+λgx(x⋆,u⋆)]λ(t1)=0Optimality Equation:fu(x⋆,u⋆)+λgu(x⋆,u⋆)=0State~Equation:\\\dot{x}=g(x,u)\\x(t_0)=x_0\\Multiplier~Equation:\\\dot\lambda=-[f_x(x^\star ,u^\star )+\lambda g_x(x^\star ,u^\star )]\\\lambda(t_1)=0\\Optimality~Equation:\\f_u(x^\star ,u^\star )+\lambda g_u(x^\star ,u^\star )=0

这就是汉密尔顿方程，普遍的表达形式为

H(t,x(t),u(t),λ(t))=f(t,x,u)+λg(t,x,u)H(t,x(t),u(t),\lambda(t))=f(t,x,u)+\lambda g(t,x,u)

并且方程满足

∂H∂u=fu+λgu=0−∂H∂x=−(fx+λgx)=λ˙∂H∂λ=g=x˙\frac{\partial H}{\partial u}=f_u+\lambda g_u=0\\\\ -\frac{\partial H}{\partial x}=-(f_x+\lambda g_x)=\dot\lambda\\\\ \frac{\partial H}{\partial \lambda}=g=\dot{x}

求解：

State Equation:x˙=g(x,u)x(t0)=x0Multiplier Equation:λ˙=−[fx(x⋆,u⋆)+λgx(x⋆,u⋆)]λ(t1)=0Optimality Equation:fu(x⋆,u⋆)+λgu(x⋆,u⋆)=0State~Equation:\\\\ \dot{x}=g(x,u)\\\\ x(t_0)=x_0\\\\ Multiplier~Equation:\\\\ \dot\lambda=-[f_x(x^\star ,u^\star )+\lambda g_x(x^\star ,u^\star )]\\\\ \lambda(t_1)=0\\\\ Optimality~Equation:\\\\ f_u(x^\star ,u^\star )+\lambda g_u(x^\star ,u^\star )=0

根据 Optimality EquationOptimality~Equation 可以得到 uu 关于 xx 和 λ\lambda 的方程，然后带入 Multiplier EquationMultiplier~Equation 和 State EquationState~Equation 中可以得到两个微分方程，并且注意两个边界条件，这就可以解出两个微分方程，从而求解

对于方程最大值问题，还要求解 Huu(t,x⋆,u⋆,λ)≤0H_{uu}(t,x^\star ,u^\star ,\lambda)\leq 0，对于最小值，需要保证 Huu(t,x⋆,u⋆,λ)≥0H_{uu}(t,x^\star,u^\star,\lambda)\geq 0

求解例子

对于一个系统

x˙=Ax+Buy=CxJ=12x(tf)TSx(tf)+12∫0tf[xTQx+uTRu]dt\dot{x}=Ax+Bu\\y=Cx\\J=\frac{1}{2} x(t_f)^TSx(t_f)+\frac{1}{2}\int_{0}^{t_f}{[x^TQx+u^TRu]}dt

求解 min Jmin~J 相当于求解 min ∫0tf[xTQx+uTRu]dtmin~\int_{0}^{t_f}{[x^TQx+u^TRu]}dt

所以引入 hamiltonhamilton 函数

H(t,x,u)=xTQx+uTRu+λT(Ax+Bu)H(t,x,u)=x^TQx+u^TRu+\lambda^T(Ax+Bu)

∂H∂u=(R+RT)u+BTλ=0⇓u=−(R+RT)−1BTλ\frac{\partial H}{\partial u}=(R+R^T)u+B^T\lambda=0\\\\ \Downarrow\\\\ u=-(R+R^T)^{-1}B^T\lambda

−∂H∂x=−(fx+λgx)=−[(Q+QT)x+ATλ]=λ˙∂H∂λ=g=Ax+Bu=Ax−B(R+RT)−1BTλ=x˙-\frac{\partial H}{\partial x}=-(f_x+\lambda g_x)=-[(Q+Q^T)x+A^T\lambda]=\dot\lambda\\\\ \frac{\partial H}{\partial \lambda}=g=Ax+Bu=Ax-B(R+R^T)^{-1}B^T\lambda=\dot{x}

由于我们设定 u=−Kxu=-Kx，所以可以知道 λ=Px\lambda=Px，并且 −(R+RT)−1BTP=K-(R+R^T)^{-1}B^TP=K

带入上述的方程得到

λ˙=−[(Q+QT)x+ATλ]=−(Q+QT+ATP)xx˙=Ax−B(R+RT)−1BTλ=[A−B(R+RT)−1BTP]x\dot\lambda=-[(Q+Q^T)x+A^T\lambda]=-(Q+Q^T+A^TP)x\\\\ \dot{x}=Ax-B(R+R^T)^{-1}B^T\lambda=[A-B(R+R^T)^{-1}B^TP]x

联立得到

λ˙=P˙x+Px˙=[P˙+P(A−B(R+RT)−1BTP)]x=−(Q+QT+ATP)x\dot{\lambda}=\dot{P}x+P\dot{x}=[\dot{P}+P(A-B(R+R^T)^{-1}B^TP)]x=-(Q+Q^T+A^TP)x

所以得到

P˙+PA−PB(R+RT)−1BTP+Q+QT+ATP=0\dot{P}+PA-PB(R+R^T)^{-1}B^TP+Q+Q^T+A^TP=0

令 Q1=Q+QT,Q2=R+RTQ_1=Q+Q^T,Q_2=R+R^T，得到

P˙+PA−PBQ2−1BTP+Q1+ATP=0\dot{P}+PA-PBQ_2^{-1}B^TP+Q_1+A^TP=0

最终得到一个 RiccatiRiccati 方程

Hamiltonian Matrix

首先引入一种特殊的矩阵

J=[0I−I0]  J∈R2n×2nJ=\begin{bmatrix}0&I\\-I&0\end{bmatrix}~~J\in R^{2n\times 2n}

这个矩阵具有如下性质

JT=−JJ−1=JTJTJ=I2nJTJT=−I2nJ2=−I2ndetJ=±1J^T=-J\\J^{-1}=J^T\\J^TJ=I_{2n}\\J^TJ^T=-I_{2n}\\J^2=-I_{2n}\\detJ=\pm1

定义 Hamiltonian matrixHamiltonian~ matrix

如果一个矩阵 A=JSA=JS，其中 JJ 是上面的特殊矩阵， S=STS=S^T，那么矩阵 AA 是一个 HamiltonianHamiltonian 矩阵

或者说，如果矩阵 AA 满足 (JA)T=JA(JA)^T=JA ， JJ 是上面的特殊矩阵，那么矩阵 AA 是一个 HamiltonianHamiltonian 矩阵

性质

  1. A,B∈Hn⇒A+B∈HnA,B\in H^n ⇒ A+B\in H^n
  2. A∈Hn⇒αA∈HnA\in H^n ⇒ \alpha A\in H^n
  3. [AB]∈Hn⇒det[AB]=AB−BA\begin{bmatrix}A&B\end{bmatrix}\in H^n⇒det\begin{bmatrix}A&B\end{bmatrix}=AB-BA
  4. 对于 A∈HnA\in H^n，定义 pA(x)pA(x) 为矩阵 A 的特征多项式
    • pA(x)=pA(−x)pA(x)=pA(-x)
    • if pA(c)=0,c∈R,then pA(−c)=pA(c‾)=pA(−c‾)=0if~pA(c)=0,c\in R,then~pA(-c)=pA(\overline{c})=pA(-\overline{c})=0