典型问题

根据系统的建模可以划分为：

• 线性系统： $\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u}$
• 非线性系统 $\dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}(t),\mathbf{u}(t),t)$
  • 可以通过一些手段进行线性化

根据约束类型可以划分为：

• 控制约束：$\mathbf{u}(t)\in \mathbb{U}\subseteq {\mathbb{R}}^m$
• 状态约束： $\mathbf{f}(\mathbf{x})\le 0$ or $\mathbf{g}(\mathbf{x})=0$

其实无论是线性系统/非线性系统，控制约束/状态约束，使用最优化方法解决时都是约束，都需要使用庞特里亚金极小值原理解决

一个最优控制的典型问题形式如下：

• 被控系统(动态约束)：$$\dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}(t),\mathbf{u}(t),t),\mathbf{x}(t_0)={\mathbf{x}}_0$$
• 目标函数：$$J=\Phi \left[\mathbf{x}(t_f)\right]+{\int }_{t_0}^{t_f}L\left[\mathbf{x}(t),\mathbf{u}(t),t\right]\mathrm{d}t$$
  • 寻找最优控制 ${\mathbf{u}}^{\ast }(t)$ 和最优轨迹${\mathbf{x}}^{\ast }(t)$，使 $J$ 最小化（或最大化），同时满足状态方程和控制约束 $\mathbf{u}(t)\in \mathbb{U}$。
• 约束条件：
  • 输入约束：$\mathbf{u}(t)\in \mathbb{U}\subseteq {\mathbb{R}}^m$
  • 状态约束，末尾状态也可以设定为一个约束

典型方法

变分法

将问题理解为求代价泛函$J[y]$极小值函数$y$,一般用于无约束问题，如果遇见约束问题，需要构建拉格朗日函数，将原问题的约束条件 “嵌入” 到新构造的增广代价泛函中，从而将有约束的优化问题转化为无约束的极值问题。
变分法基础
变分法研究泛函的极值，泛函以函数为自变量（如曲线长度、物理作用量）。通过变分$\delta y(x)$，函数的微小扰动），将泛函极值转化为微分方程求解。

• 泛函示例：$$J[y]={\int }_{x_1}^{x_2}F(x,y,y^{\prime })dx$$其中 $ y(x)$ 是函数，$y^{\prime }=dy/dx$，求函数$y^{\ast }(x)$使得泛函取得最极值。
• 变分：泛函的自变量（函数）的微小变化称为变分，记作$\delta y$，
• 变分与微分/积分满足可交换性：$\delta \left(\frac{dy}{dx}\right)=\frac{d}{dx}(\delta y),\delta \int Ldx=\int \delta Ldx$
• 变分满足链式法则， $\delta F=\frac{\partial F}{\partial x}\delta x+\frac{\partial F}{\partial y}\delta y+\frac{\partial F}{\partial \dot{y}}\delta \dot{y}$，其中根据变分与微分的交换律： $\delta \dot{y}=\frac{d}{dx}(\delta y)$ ,自变量x的变分$\delta x=0$

欧拉-拉格朗日（E-L）方程

E-L 方程描述了极值曲线在函数空间中的 “驻点”，即沿任何变分方向的泛函变化率为零，类似于函数极值的一阶导数为零。
E-L方程是泛函极值的必要条件，推导如下：

• 邻近曲线：设 $y(x)=y(x)+ϵ\delta y(x)$，$\delta$就是变分，需要任意可微，为了固定边界确保端点不变还需要满足$\delta (x_1)=\delta (x_2)=0$，$ϵ\to 0$。
• 此时泛函可写做：$$J(y)={\int }_{x_1}^{x_2}F(x,y+ϵ\delta y,y^{\prime }+ϵ\delta y^{\prime })dx$$
• 泛函导数：对 $J[ϵ]$ 求导并令 $ϵ=0$，分部积分后得：$$\int \delta \left(\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)\right)dx=0$$
  • 推导：利用莱布尼兹法则，对 $ϵ$求导保留线性部分：$$\delta J(ϵ)={\int }_{x_1}^{x_2}\left(\frac{\partial F}{\partial (y+ϵ\delta y)}\cdot \delta y+\frac{\partial F}{\partial (y^{\prime }+ϵ\delta y^{\prime })}\cdot \delta y^{\prime }\right)dx$$
  • 根据变分定义，令$ϵ=0$，得到泛函在$y(x)$处沿 $\eta (x)$ 方向的变分导数： $$\delta J={\int }_{x_1}^{x_2}\left(\frac{\partial F}{\partial y}\delta +\frac{\partial F}{\partial y^{\prime }}\delta y^{\prime }\right)dx$$
  • 然后取出后半部分进行分部积分以消除$\delta y^{\prime }$：$${\int }_{x_1}^{x_2}\frac{\partial F}{\partial y^{\prime }}\delta y^{\prime }dx={\left[\frac{\partial F}{\partial y^{\prime }}\delta y\right]}_{x_1}^{x_2}-{\int }_{x_1}^{x_2}\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)\delta ydx$$
  • 由于$\delta y(x_1)=\delta y(x_2)=0$,带入后化简可得：$$\delta J={\int }_{x_1}^{x_2}\left(\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)\right)\delta ydx=0$$
• 由于$\delta$任意，根据变分引理则被积函数必恒为零。故 $E-L方程$：$$\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)=0$$
  • 变分引理：若连续函数$G(x)$对于任意$\eta \in C^2([x_1,x_2])$满足$\int G\eta dx=0$，则$G(x)\equiv 0$

求解E-L方程就可以得到极值处的函数 $y^{\ast }(x)$，这个过程可以结合导数的概念来理解，,函数空间中一个点$y(x)$，临近曲线 $y(x)+ϵ\delta y$就是在这个点附近的点，且$ϵ\to 0$，泛函对$ϵ$求导得到的就是$y(x)$这个点在$\delta y(x)$这个方向的导数，满足朝向任何方向$\delta y(x)$上的导数都为0的点自然就是极值点。最优控制的庞特里亚金原理中的哈密顿函数，本质是变分法在动态系统的扩展（含协态变量的E-L方程推广）。

• 如果F含高阶导数，则E-L方程也要扩展，$$\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)+\frac{d^2}{d{x}^2}\left(\frac{\partial F}{\partial y^{\prime \prime }}\right)-\cdots =0$$
• 如果存在约束，则需要使用拉格朗日乘子构造增广泛函

拉格朗日函数（Lagrange Function）

拉格朗日函数（Lagrange Function）是数学优化理论中的核心工具，用于求解带有约束条件的极值问题。它通过引入拉格朗日乘数（Lagrange Multiplier），将原问题转化为无约束的极值问题，从而利用微积分中的求导方法求解。以下是详细介绍：
假设我们需要优化（最大化或最小化）一个目标函数 $$f(\mathbf{x})$$，同时满足若干等式约束 $$g_i(\mathbf{x})=0（i=1,2,\dots ,m）$$。
拉格朗日函数的核心思想是：通过引入拉格朗日乘数 ${\lambda }_i$，将约束条件与目标函数“合并”为一个新的函数，使得原问题的极值可以通过求解新函数的无约束极值得到。
设目标函数为 $f(\mathbf{x})$，约束条件为 $m$ 个等式约束 $g_i(\mathbf{x})=0,(\mathbf{x}\in {\mathbb{R}}^n)$，且 $m<n$。
拉格朗日函数定义为：
$$L(\mathbf{x},{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_m)=f(\mathbf{x})-\sum _{i=1}^m{\lambda }_i{g}_i(\mathbf{x})$$
其中：

• $\mathbf{x}=(x_1,x_2,\dots ,x_n)$是优化变量，
• ${\lambda }_i$是拉格朗日乘数（可以是常数，也可以是函数，取决于优化变量，如果优化变量是一个函数，那拉格朗日乘数就也是个变量，此时的拉格朗日乘数被称为协态，拉格朗日函数称作动态拉格朗日函数，也叫哈密顿函数）。
• 拉格朗日函数基于对偶原理—— 将原始问题的约束转化为对偶变量（乘子 / 协态）的优化条件，从而在增广空间中消除显式约束。

求解原问题的极值等价于求解拉格朗日函数$L$的无约束极值。根据微积分理论，极值点处的梯度（各偏导数）必须为零，即：

• 对优化变量 $x_i$ 求偏导并令其为零：
  $$\frac{\partial L}{\partial x_j}=\frac{\partial f}{\partial x_j}-\sum _{i=1}^m{\lambda }_i\frac{\partial g_i}{\partial x_j}=0(j=1,2,\dots ,n)$$
• 对拉格朗日乘数 ${\lambda }_i$ 求偏导并令其为零（自动满足约束条件）：
  $$\frac{\partial L}{\partial {\lambda }_i}=-g_i(\mathbf{x})=0(i=1,2,\dots ,m)$$

这组方程（共 $m+n$个）称为KKT条件（对等式约束情形，KKT条件退化为拉格朗日乘数法的条件），其解 $({\mathbf{x}}^{\ast },{\lambda }^{\ast })$即为原问题的候选极值点。
若约束条件包含不等式约束 $h_i(\mathbf{x})\le 0$，则需要使用扩展的拉格朗日函数和KKT条件，此时拉格朗日函数为：
$$L(\mathbf{x},\lambda ,\mu )=f(\mathbf{x})-\sum _{i=1}^m{\lambda }_i{g}_i(\mathbf{x})-\sum _{j=1}^k{\mu }_j{h}_j(\mathbf{x})$$
其中 ${\mu }_j\ge 0$ 是对应不等式约束的拉格朗日乘数，且需满足互补松弛条件：
$${\mu }_j{h}_j(\mathbf{x})=0$$(即约束 active 时 $h_j=0$，否则 ${\mu }_j=0$)

Pontryagin 极小值原理（庞特里亚金极小值原理）

庞特里亚金极小值原理（Pontryagin’s Minimum Principle）是最优控制理论中的核心定理之一，由苏联数学家列夫·庞特里亚金（Lev Pontryagin）及其团队在20世纪50年代提出。用于求解带有约束的动态系统最优控制问题，尤其适用于控制变量存在不等式约束（如幅值限制）的场景，是变分法的扩展和推广，也被称为现代变分法。它通过构造哈密顿函数，将最优控制问题转化为一组微分方程（正则方程）的求解，并给出控制变量的最优条件。

核心推导

问题描述
代价泛函：$J_a={\int }_{t_0}^{t_f}L(x,u,t)dt+\Phi (x(t_f))$
系统约束： $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u},t)$
首末约束： $\mathbf{x}(t_0)={\mathbf{x}}_0$

1. 构造哈密顿函数（Hamiltonian）
   引入协态变量（伴随变量） $\mathbf{\lambda }(t)\in {\mathbb{R}}^n$（与状态变量维数相同），定义哈密顿函数：$$H(\mathbf{x},\mathbf{u},\mathbf{\lambda },t)=L(\mathbf{x},\mathbf{u},t)+{\mathbf{\lambda }}^T\mathbf{f}(\mathbf{x},\mathbf{u},t)$$
   • 协态变量 $\mathbf{\lambda }(t)$ 的物理意义可理解为状态变量的“影子价格”或“边际代价”。
   • 如果有状态约束，还需要再加一项，变成：$$H(\mathbf{x},\mathbf{u},\mathbf{\lambda },t)=L(\mathbf{x},\mathbf{u},t)+{\mathbf{\lambda }}^T\mathbf{f}(\mathbf{x},\mathbf{u},t)+{\mathbf{\eta }}^T(\mathbf{x}-{\mathbf{x}}_{\mathbf{m}})$$
2. 正则方程（协态方程与状态方程）
   根据极小值原理，最优解需满足以下正则方程：

• 状态方程（与原系统一致）：$\frac{\partial H}{\partial \mathbf{\lambda }}=\mathbf{f}(\mathbf{x},\mathbf{u},t)=\dot{\mathbf{x}}(t)$
  • 对$\lambda$求变分分析推导出来
• 协态方程（伴随方程）： $-\frac{\partial H}{\partial \mathbf{x}}=-\left[\frac{\partial L}{\partial \mathbf{x}}+{\left(\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\right)}^T\mathbf{\lambda }\right]=\dot{\mathbf{\lambda }}(t)$
  • 根据变分和积分的交换律，代价函数的变分$$\delta J={\int }_{t_0}^{t_f}\delta L(x,u,t)dt={\int }_{t_0}^{t_f}\left(\frac{\partial L}{\partial x}\delta x+\frac{\partial L}{\partial u}\delta u\right)dt$$
  • 构造增广代价泛函： $$J_a={\int }_{t_0}^{t_f}\left[L(\mathbf{x},\mathbf{u},t)+{\lambda }^T(t)\mathbf{f}(\mathbf{x},\mathbf{u},t)\right]dt+\Phi (\mathbf{x}(t_f))$$
  • 将其中的 ${\lambda }^T\dot{x}$部分进行分部积分得：$$J_a={\int }_{t_0}^{t_f}\left[L+{\mathbf{\lambda }}^{T}f+{\dot{\mathbf{\lambda }}}^{T}x\right]dt-{\mathbf{\lambda }}^T(t_f)x(t_f)+{\mathbf{\lambda }}^T(t_0)x(t_0)+\Phi (x(t_f))$$
  • 带入哈密顿函数得： $$J_a={\int }_{t_0}^{t_f}\left[H+{\dot{\mathbf{\lambda }}}^{T}x\right]dt-{\mathbf{\lambda }}^T(t_f)x(t_f)+{\mathbf{\lambda }}^T(t_0)x(t_0)+\Phi (x(t_f))$$
  • 对增广泛函进行变分：$$\delta J_a={\int }_{t_0}^{t_f}\left[\frac{\partial H}{\partial x}\delta x+\frac{\partial H}{\partial u}\delta u+{\dot{\lambda }}^T\delta x\right]dt+\left[-{\lambda }^T(t_f)\delta x(t_f)+\frac{\partial \Phi }{\partial x(t_f)}\delta x(t_f)\right]=0$$
  • 分离变分项
    • 状态变分 $\delta x$：$${\int }_{t_0}^{t_f}\left(\frac{\partial H}{\partial x}+{\dot{\lambda }}^T\right)\delta xdt+\left[\left(-{\lambda }^T(t_f)+\frac{\partial \Phi }{\partial x(t_f)}\right)\delta x(t_f)\right]=0$$
    • 控制变分 $\delta u$：${\int }_{t_0}^{t_f}\frac{\partial H}{\partial u}\delta udt=0$
  • 对任意 $\delta x$ ，变分 $\delta J_a=0$ 需满足以下条件：
    • 积分项为0：$\frac{\partial H}{\partial x}+{\dot{\lambda }}^T=0\Rightarrow \dot{\lambda }=-{\left(\frac{\partial H}{\partial x}\right)}^T$这就是协态方程
    • 剩余项为0：$-\lambda (t_f)+\frac{\partial \Phi }{\partial x(t_f)}=0\Rightarrow \lambda (t_f)=\frac{\partial \Phi }{\partial x(t_f)}$
  • 对任意 $\delta u$，变分 $\delta J_a=0$ 需满足以下条件：
    • $$\frac{\partial H}{\partial u}=0\text{（无约束时）}\text{或}{u}^{\ast }=\arg \underset{u\in U}{\min }H\text{（有约束时）}$$即：无约束时直接求导为0，有约束时必须是最优控制曲线， $u^{\ast }$为函数空间中的 “驻点”，即沿任何变分方向的泛函变化率为零。

3. 终端条件
   根据终端状态是否固定，终端条件分为：
   • 终端状态固定（$\mathbf{x}(t_f)={\mathbf{x}}_f$ 已知）： $$\mathbf{\lambda }(t_f)=\frac{\partial \varphi }{\partial \mathbf{x}(t_f)}$$
     • 推导见上方协态方程部分
   • 终端状态自由：
     $$\mathbf{\lambda }(t_f)=\frac{\partial \varphi }{\partial \mathbf{x}(t_f)},H(t_f)=-\frac{\partial \varphi }{\partial t_f}(\text{若 }{t}_f\text{ 自由})$$
4. 极小值条件（核心条件）
   对于所有 $t\in [t_0,t_f]$ 和 $\mathbf{u}(t)\in \mathbb{U}$，最优控制 ${\mathbf{u}}^{\ast }(t)$ 需满足： $$H({\mathbf{x}}^{\ast }(t),{\mathbf{u}}^{\ast }(t),{\mathbf{\lambda }}^{\ast }(t),t)=\underset{\mathbf{u}\in \mathbb{U}}{\min }H({\mathbf{x}}^{\ast }(t),\mathbf{u},{\mathbf{\lambda }}^{\ast }(t),t)$$
   即哈密顿函数在最优控制处取得极小值（若目标为最大化 $J$，则取极大值）。这一条件体现了控制约束 $\mathbf{u}\in \mathbb{U}$ 的作用，当控制无约束时，可通过求导 $\frac{\partial H}{\partial \mathbf{u}}=0$ 直接求解。

步骤总结

1. 构造哈密顿函数 $H$，包含状态、控制、伴随变量和时间。
2. 建立正则方程：

• 解状态方程 $\dot{\mathbf{x}}=\frac{\partial H}{\partial \mathbf{\lambda }}$，
• 解协态方程 $\dot{\mathbf{\lambda }}=-\frac{\partial H}{\partial \mathbf{x}}$。

3. 确定终端条件：根据终端状态和时间是否固定，设定 $\mathbf{\lambda }(t_f)$ 和 $t_f$ 的条件。
4. 求解极小值条件：利用 $H$ 对 $\mathbf{u}$ 的极小化（或极大化），结合控制约束 $\mathbb{U}$，求出 ${\mathbf{u}}^{\ast }(t)$。
5. 联立方程求解：将 ${\mathbf{u}}^{\ast }(t)$ 代入正则方程，解两点边值问题（BVP）得到 ${\mathbf{x}}^{\ast }(t)$ 和 ${\mathbf{\lambda }}^{\ast }(t)$。

简单示例：最速降线问题

以下将通过一个典型的双积分器系统最短时间控制问题，详细讲解如何运用庞特里亚金极小值原理求解最优控制问题。我们将从问题建模、哈密顿函数构造、协态方程推导、控制律求解到最优轨迹分析逐步展开，并使用$$环境呈现公式。

问题设定：双积分器系统的最短时间控制
考虑如下二阶线性系统（双积分器）： $\{\begin{array}{l}{\dot{x}}_1=x_2,\\ {\dot{x}}_2=u,\end{array}$
其中：

• $x_1(t),x_2(t)$ 为状态变量（位置和速度），
• $u(t)$ 为控制输入，满足幅值约束 $|u(t)|\le 1$，
• 目标：从初始状态 $(x_{10},x_{20})$ 转移到终端状态 $(0,0)$，且时间最短。

系统状态方程建模为：
$$\left[\begin{array}{c}\dot{x_1}\\ \dot{x_2}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]+\left[\begin{array}{c}0\\ 1\end{array}\right]u$$

性能指标为：
$$J=\underset{u}{\min }{\int }_0^{T}1dt=T\text{（最小化时间）}.$$

第一步：构造哈密顿函数
根据庞特里亚金极小值原理，定义哈密顿函数 $H$ 为：
$$\begin{array}{rl}H& =L+{\mathbf{\lambda }}^T\dot{x}\\ & ={\lambda }_1{\dot{x}}_1+{\lambda }_2{\dot{x}}_2+L\end{array}$$
其中：

• ${\lambda }_1(t),{\lambda }_2(t)$ 为协态变量（伴随变量），
• $L$ 为性能指标的被积函数（此处 $L=1$）。

代入状态方程 ${\dot{x}}_1=x_2$ 和 ${\dot{x}}_2=u$，得：
$$H={\lambda }_1{x}_2+{\lambda }_{2}u+1.$$
第二步：推导协态方程
根据PMP，协态变量的导数满足：
$${\dot{\lambda }}_i=-\frac{\partial H}{\partial x_i}(i=1,2).$$
计算偏导数：

• 对 $x_1$：$\frac{\partial H}{\partial x_1}=0$，故 ${\dot{\lambda }}_1=0$，即 ${\lambda }_1$ 为常数，记为 ${\lambda }_1=c_1$；
• 对 $x_2$：$\frac{\partial H}{\partial x_2}={\lambda }_1$，故 ${\dot{\lambda }}_2=-{\lambda }_1=-c_1$。

因此，协态变量 ${\lambda }_2(t)$ 是关于时间的线性函数：
$${\lambda }_2(t)=-c_{1}t+c_2,$$
其中 $c_1,c_2$ 为常数，由边界条件确定。
第三步：利用极小值条件求解控制律
PMP要求控制输入 $u(t)$ 在每时刻 $t$ 最大化哈密顿函数 $H$，即：
$$u^{\ast }(t)=\arg \underset{|u|\le 1}{\min }H=\arg \underset{|u|\le 1}{\min }\left({\lambda }_{2}u+\text{其他}\right).$$
由于 ${\lambda }_2$ 是线性函数，最大化 $H$ 等价于取 $u$ 与 ${\lambda }_2$ 同号。因此：
$$u^{\ast }(t)=\text{sign}({\lambda }_2(t))=\{\begin{array}{ll}+1,& {\lambda }_2(t)<0,\\ -1,& {\lambda }_2(t)>0.\end{array}$$
当 ${\lambda }_2(t)=0$ 时，可能存在切换点，此时控制输入发生翻转。

第四步：分析协态变量与切换特性
由于 ${\lambda }_2(t)=-c_{1}t+c_2$ 是线性函数，其符号最多改变一次，因此最优控制 $u^{\ast }(t)$ 最多切换一次，属于Bang-Bang控制。分两种情况讨论：

1. 无切换（直接Bang）：若 ${\lambda }_2(t)$ 始终非负或非正，则 $u^{\ast }(t)$ 保持 $+1$ 或 $-1$ 不变；
2. 一次切换（Bang-Bang）：若 ${\lambda }_2(t)$ 由正变负或由负变正，则存在唯一切换时刻 $\tau$，使得：
   $$u^{\ast }(t)=\{\begin{array}{ll}+1,& t<\tau ,\\ -1,& t\ge \tau \end{array}\text{或}\{\begin{array}{ll}-1,& t<\tau ,\\ +1,& t\ge \tau .\end{array}$$

第五步：利用终端条件求解常数
终端状态固定为 $(0,0)$，终端时间 $T$ 自由。根据PMP，终端时刻的横截条件为：$H(T)=0.$代入 $H$ 的表达式：
$$H(T)={\lambda }_1{x}_2(T)+{\lambda }_2(T)u(T)+1=0.$$
由于 $x_2(T)=0$（终端速度为0），化简得：
$${\lambda }_2(T)u(T)+1=0.$$
又因 $u(T)=\text{sign}({\lambda }_2(T))$，故：
$$|{\lambda }_2(T)|=1\Rightarrow {\lambda }_2(T)=\pm 1.$$
结合 ${\lambda }_2(t)=-c_{1}t+c_2$，可得：
$${\lambda }_2(T)=-c_{1}T+c_2=\pm 1.$$

第六步：求解状态轨迹与最优控制
以先正后负切换（+1→-1为例，假设切换时刻为 $\tau$，则：

• 当 $t\in [0,\tau ]$ 时，$u=+1$，状态方程积分得：
  $$\{\begin{array}{l}{x}_1(t)=x_{10}+x_{20}t+\frac{1}{2}{t}^2,\\ x_2(t)=x_{20}+t.\end{array}$$
• 当 $t\in [\tau ,T]$ 时，$u=-1$，状态方程积分得：
  $$\{\begin{array}{l}{x}_1(t)=x_1(\tau )+x_2(\tau )(t-\tau )-\frac{1}{2}(t-\tau {)}^2,\\ x_2(t)=x_2(\tau )-(t-\tau ).\end{array}$$
  终端条件 $x_1(T)=0,x_2(T)=0$ 给出方程组，可解出 $\tau$ 和 $T$。

第七步：相平面分析与开关曲线
最优轨迹在相平面 $(x_1,x_2)$ 上的形状由开关曲线决定。开关曲线 $\Gamma$ 分为上下两支：
$$\Gamma ={\Gamma }^{+}\cup {\Gamma }^{-},$$
其中：

• ${\Gamma }^{+}$（对应 $u=-1$ 直接到达原点）：$x_1=-\frac{1}{2}{x}_2^2$（当 $x_2\le 0$），
• ${\Gamma }^{-}$（对应 $u=+1$ 直接到达原点）：$x_1=\frac{1}{2}{x}_2^2$（当 $x_2\ge 0$）。

若初始状态在 $\Gamma$ 上方，则最优控制为先 $+1$ 后 $-1$；若在下方，则为先 $-1$ 后 $+1$；若在 $\Gamma$ 上，则无需切换。

具体实现代码：

import numpy as np
import matplotlib.pyplot as pltplt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = Falsedef solve_min_time_control(x10, x20):"""求解双积分器系统的最短时间控制参数（切换时间τ和总时间T）参数：x10: 初始位置x20: 初始速度返回：τ: 切换时间T: 总时间u_seq: 控制序列（分段函数）"""# 开关曲线判断：x1 = 0.5*x2²（上支）或x1 = -0.5*x2²（下支）# 初始状态在开关曲线上方时（x10 > 0.5*x20²），切换策略为 -1 → +1（先制动后加速）if (x20 > 0 and x10 > -0.5 * x20**2) or (x20 < 0 and x10 > 0.5 * x20**2):# 推导方程：τ² - 2x20τ + (0.5x20² - x10) = 0（详细推导见说明）a = 1b = -2 * x20c = 0.5 * x20**2 - x10discriminant = b**2 - 4*a*cif discriminant < 0:raise ValueError("无实数解，初始状态无法到达原点")τ = ( -b + np.sqrt(discriminant) ) / (2*a)  # 取较大的根（保证物理意义）T = 2*τ - x20  # 由x2(T)=0推导得到# 定义分段控制律：前τ时间用-1，之后用+1def u(t):if t < τ:return -1.0else:return 1.0return τ, T, u# 初始状态在开关曲线下方时（x10 < -0.5*x20²），切换策略为 +1 → -1（先加速后制动）elif (x20 > 0 and x10 < -0.5 * x20**2) or (x20 < 0 and x10 < 0.5 * x20**2):# 推导方程：τ² + 2x20τ + (x10 + 0.5x20²) = 0a = 1b = 2 * x20c = x10 + 0.5 * x20**2discriminant = b**2 - 4*a*cif discriminant < 0:raise ValueError("无实数解，初始状态无法到达原点")τ = (-b - np.sqrt(discriminant)) / (2*a)  # 取较小的根T = 2*τ + x20  # 由x2(T)=0推导得到def u(t):if t < τ:return 1.0else:return -1.0return τ, T, u# 初始状态在开关曲线上（无需切换）else:if x20 >= 0:# 上开关曲线（x1=0.5x2²）且速度非负时，直接用-1制动T = x20 + np.sqrt(x20**2 - 2*x10)  # 推导自x2(T)=0和x1(T)=0def u(t):return -1.0else:# 下开关曲线（x1=-0.5x2²）且速度负时，直接用+1加速T = -x20 + np.sqrt(x20**2 + 2*x10)def u(t):return 1.0return 0, T, udef simulate_trajectory(x10, x20, τ, T, u):"""模拟状态轨迹（分切换前后两段积分）"""t1 = np.linspace(0, τ, 100)t2 = np.linspace(τ, T, 100)t_total = np.concatenate([t1, t2])# 切换前控制（u=-1或+1）if u(0) == -1:  # 上开关曲线场景：先-1制动x1_1 = x10 + x20*t1 - 0.5*t1**2x2_1 = x20 - t1else:  # 下开关曲线场景：先+1加速x1_1 = x10 + x20*t1 + 0.5*t1**2x2_1 = x20 + t1# 切换后控制（u=+1或-1）x1_τ = x1_1[-1]  # 切换时刻的x1x2_τ = x2_1[-1]  # 切换时刻的x2dt = t2 - τif u(τ+1e-3) == 1:  # 切换后为+1x1_2 = x1_τ + x2_τ*dt + 0.5*dt**2x2_2 = x2_τ + dtelse:  # 切换后为-1x1_2 = x1_τ + x2_τ*dt - 0.5*dt**2x2_2 = x2_τ - dtx1_total = np.concatenate([x1_1, x1_2])x2_total = np.concatenate([x2_1, x2_2])return t_total, x1_total, x2_totaldef visualize_results(x10, x20, τ, T, t_total, x1_total, x2_total, u):"""可视化相平面轨迹和时间响应"""plt.figure(figsize=(12, 6))# 相平面轨迹（x1 vs x2）plt.subplot(1, 2, 1)plt.plot(x2_total, x1_total, 'b-', label='最优轨迹')plt.scatter(x20, x10, c='r', marker='o', label='初始状态')plt.scatter(0, 0, c='g', marker='s', label='终点（原点）')x2_upper = np.linspace(-5, 0, 100)x2_lower = np.linspace(0, 5, 100)x1_switch_upper = 0.5 * x2_upper**2x1_switch_lower = -0.5 * x2_lower**2plt.plot(x2_upper, x1_switch_upper, 'b--', label='上切换曲线u=1')plt.plot(x2_lower, x1_switch_lower, 'r--', label='下切换曲线u=-1')plt.xlabel('速度 $x_2$')plt.ylabel('位置 $x_1$')plt.title('相平面轨迹')plt.legend()plt.grid(True)# 时间响应（x1(t), x2(t), u(t)）plt.subplot(1, 2, 2)plt.plot(t_total, x1_total, 'r-', label='$x_1(t)$')plt.plot(t_total, x2_total, 'g-', label='$x_2(t)$')u_values = [u(t) for t in t_total]plt.step(t_total, u_values, 'b--', label='$u(t)$', where='post')plt.xlabel('时间 $t$')plt.ylabel('状态/控制')plt.title('时间响应曲线')plt.legend()plt.grid(True)plt.tight_layout()plt.show()# 主程序运行（示例初始状态）
if __name__ == "__main__":x10, x20 = 2.0, 4.0  # 初始状态（位置=2，速度=1）τ, T, u = solve_min_time_control(x10, x20)t_total, x1_total, x2_total = simulate_trajectory(x10, x20, τ, T, u)print(f"切换时间 τ = {τ:.2f}s，总时间 T = {T:.2f}s")visualize_results(x10, x20, τ, T, t_total, x1_total, x2_total, u)