变分法求一个最优控制问题

终端时刻固定、终端状态固定

$x^{\prime }=-x+u$，$x(0)=2$，$x(1)=1$，$J=0.5{\int }_0^1(x^2+u^2)\text{d}t$，求使$J$最小的解u*(t)。

经典变分定义法

设$x(t)=x_0(t)+ϵ\eta (t)$，则
$$\begin{array}{rl}J& =0.5{\int }_0^{1}f(x,x^{\prime },u)\text{d}t=0.5{\int }_0^1((x_0+ϵ\eta {)}^2+u^2)\text{d}t\\ \frac{\text{d}J}{\text{d}ϵ}& ={\int }_0^1(\frac{\partial f}{\partial x}\frac{\partial x}{\partial ϵ}+\frac{\partial f}{\partial u}\frac{\partial u}{\partial ϵ})\text{d}t\\ & ={\int }_0^1(x\eta +u({\eta }^{\prime }+\eta ))\text{d}t\\ & ={\int }_0^1(x\eta +(x+x^{\prime })({\eta }^{\prime }+\eta ))\text{d}t\\ & ={\int }_0^1((2x+x^{\prime })\eta +(x+x^{\prime }){\eta }^{\prime })\text{d}t\\ & ={\int }_0^1(2x+x^{\prime })\eta \text{d}t+{\int }_0^1(x+x^{\prime })\text{d}\eta \\ & ={\int }_0^1(2x-x^{\prime \prime })\eta \text{d}t+{\left[(x+x^{\prime })\eta \right]}_0^1\end{array}$$
因为$\eta (0)=\eta (1)$，所以$2x-x^{\prime \prime }=0$，$x(t)=C_1{\text{e}}^{\sqrt{2}t}+C_2{\text{e}}^{-\sqrt{2}t}$，代入起止条件解得
$$x(t)=0.13275162522663497{\text{e}}^{\sqrt{2}t}+1.867248374773365{\text{e}}^{-\sqrt{2}t}$$

经典变分法

$$\begin{array}{rl}& 2f=x^2+u^2=x^2+(x+x^{\prime }{)}^2=x^{\prime 2}+2x{x}^{\prime }+2x^2\\ & \frac{\partial f}{\partial x}=\frac{\text{d}}{\text{d}t}\frac{\partial f}{\partial x^{\prime }}\\ & x^{\prime }+2x=(x+x^{\prime }{)}^{\prime }\\ & x^{\prime \prime }=2x\end{array}$$

哈密顿方程组法

$$\begin{array}{rl}& H=f(x,x^{\prime },u)+p(x)g(x)=0.5(x^2+u^2)+p(-x+u)\\ & \frac{\partial H}{\partial u}=0=u+p\\ & \frac{\partial H}{\partial p}=x^{\prime }=-x+u\\ & \frac{\partial H}{\partial x}=-p^{\prime }=x-p\end{array}$$
解得$x^{\prime \prime }=2x$。

仿真结果

x: 1
x’: 0.13022268764
u: 1.13022268764
J: 1.01806017092

终端时刻固定、终端状态自由

在前面问题的基础上，$x(1)$为任意值。由$x(1)+x^{\prime }(1)=0$计算得到$x$。

仿真结果

x: 0.563939069277
x’: -0.563939069277
u: -3.10862446895e-14
J: 0.771637192373

相关代码

import numpy as np
from scipy import integrate
def f(t):
    tx = np.sqrt(2)
    ty = np.exp(tx)
    C2 = (2*ty-1)/(ty-1/ty)
    C1 = 2-C2
    x = C1*np.exp(tx*t) + C2*np.exp(-tx*t)
    x1 = C1*tx*np.exp(tx*t) - C2*tx*np.exp(-tx*t)
    u = x + x1
    return 0.5*(u*u + x*x)
v, err=integrate.quad(f, 0, 1)
print(v, err)

#include <iostream>
#include "simucpp.hpp"
using namespace simucpp;
using namespace std;

class Funcu: public UserFunc
{
public:
    virtual double Function(double t) const {
        double tx = sqrt(2);
        double ty = exp(tx);
        double C2 = (2*ty-1)/(ty-1/ty);
        double C1 = 2-C2;
        // double C1 = 2/(1+exp(2*sqrt(2)));
        // double C2 = 2 - C1;
        double x = C1*exp(tx*t) + C2*exp(-tx*t);
        double x1 = C1*tx*exp(tx*t) - C2*tx*exp(-tx*t);
        return x + x1;
    }
};
class FuncJ: public UserFunc
{
public:
    virtual double Function(double *u) const {
        return 0.5*(u[0]*u[0]+u[1]*u[1]);
    }
};

int main()
{
    Simulator sim1(1);
    FMInput(minu, &sim1);
    FMSum(msumx, &sim1);
    FMFcnMISO(mfcnJ, &sim1);
    FMIntegrator(mintx, &sim1);
    FMIntegrator(mintJ, &sim1);
    FMOutput(moutu, &sim1);
    FMOutput(moutx, &sim1);
    sim1.connect(minu, msumx);
    sim1.connect(msumx,mintx);
    sim1.connect(mintx, msumx);
    msumx->Set_InputGain(-1);
    sim1.connect(mintx, mfcnJ);
    sim1.connect(minu, mfcnJ);
    sim1.connect(mfcnJ, mintJ);
    sim1.connect(minu, moutu);
    sim1.connect(mintx, moutx);
    Funcu *funcu = new Funcu;
    FuncJ *funcJ = new FuncJ;
    minu->Set_Function(funcu);
    mfcnJ->Set_Function(funcJ);
    mintx->Set_InitialValue(2);
    sim1.Initialize();
    sim1.Simulate();
    cout.precision(12);
    cout << "x" <<": "<< mintx->Get_OutValue() <<endl;
    cout << "x'" <<": "<< msumx->Get_OutValue() <<endl;
    cout << "u" <<": "<< minu->Get_OutValue() <<endl;
    cout << "J" <<": "<< mintJ->Get_OutValue() <<endl;
    sim1.Plot();
    return 0;
}

参考

python3通过scipy实现定积分
张杰,王飞跃. 最优控制:数学理论与智能方法(上册).