最優控制：從變分法到龐特里亞金原理

典型問題

根據系統的建模可以劃分為：

線性系統： $\mathbf{\dot{x}} = \boldsymbol{A}\mathbf{x}+\boldsymbol{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), \quad \mathbf{x}(t_0) = \mathbf{x}_0$
目標函數： $\varPhi \left[ \boldsymbol{x}(t_f) \right] + \int_{t_0}^{t_f} L \left[ \boldsymbol{x}(t), \boldsymbol{u}(t), t \right] \mathrm{d}t$
- 尋找最優控制 $\mathbf{u}^*(t)$ 和最優軌跡 $\mathbf{x}^*(t)$ ，使 $J$ 最小化（或最大化），同時滿足狀態方程和控制約束 $\mathbf{u}(t) \in \mathbb{U}$ 。
約束條件：
- 輸入約束： $\mathbf{u}(t) \in \mathbb{U} \subseteq \mathbb{R}^m$
- 狀態約束，末尾狀態也可以設定為一個約束

典型方法

變分法

將問題理解為求代價泛函 $J [y]$ 極小值函數 $y$ ,一般用于無約束問題，如果遇見約束問題，需要構建拉格朗日函數，將原問題的約束條件 “嵌入” 到新構造的增廣代價泛函中，從而將有約束的優化問題轉化為無約束的極值問題。
變分法基礎
變分法研究泛函的極值，泛函以函數為自變量（如曲線長度、物理作用量）。通過變分 $\delta y(x)$ ，函數的微小擾動），將泛函極值轉化為微分方程求解。

泛函示例： $\int_{x_1}^{x_2} F(x, y, y') dx$ 其中 $ y(x)$ 是函數， $y^{'} = d y / d x$ ，求函數 $y^*(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方程是泛函極值的必要條件，推導如下：

鄰近曲線：設 $\epsilon\delta y(x)$ ， $\delta$ 就是變分，需要任意可微，為了固定邊界確保端點不變還需要滿足 $\delta(x_1)=\delta(x_2)=0$ ， $\epsilon \to 0$ 。
此時泛函可寫做： $\int_{x_1}^{x_2} F(x, y + \epsilon \delta y, y' + \epsilon \delta y') \, dx$
泛函導數：對 $J[\epsilon]$ 求導并令 $\epsilon=0$ ，分部積分后得： $\int \delta \left( \frac{\partial F}{\partial y} - \frac{d}{dx}\left(\frac{\partial F}{\partial y'}\right) \right) dx = 0$
- 推導：利用萊布尼茲法則，對 $\epsilon$ 求導保留線性部分： $\delta J(\epsilon) = \int_{x_1}^{x_2} \left( \frac{\partial F}{\partial (y + \epsilon \delta y)} \cdot \delta y + \frac{\partial F}{\partial (y' + \epsilon \delta y')} \cdot \delta y' \right) dx$
- 根據變分定義，令 $\epsilon = 0$ ，得到泛函在 $y (x)$ 處沿 $\eta(x)$ 方向的變分導數： $\delta J = \int_{x_1}^{x_2} \left( \frac{\partial F}{\partial y} \delta + \frac{\partial F}{\partial y'} \delta y' \right) dx$
- 然后取出后半部分進行分部積分以消除 $\delta y'$ ： $\int_{x_1}^{x_2} \frac{\partial F}{\partial y'} \delta y' \, dx = \left[ \frac{\partial F}{\partial y'} \delta y \right]_{x_1}^{x_2} - \int_{x_1}^{x_2} \frac{d}{dx}\left( \frac{\partial F}{\partial y'} \right) \delta y \, dx$
- 由于 $\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'} \right) \right) \delta y \, dx = 0$
由于 $\delta$ 任意，根據變分引理則被積函數必恒為零。故 $E ? L 方程$ ： $\frac{\partial F}{\partial y} - \frac{d}{dx}\left(\frac{\partial F}{\partial y'}\right) = 0$
- 變分引理：若連續函數 $G (x)$ 對于任意 $\eta \in C^2([x_1,x_2])$ 滿足 $\int G\eta dx = 0$ ，則 $\equiv 0$

求解E-L方程就可以得到極值處的函數 $y^*(x)$ ，這個過程可以結合導數的概念來理解，,函數空間中一個點 $y (x)$ ，臨近曲線 $\epsilon\delta y$ 就是在這個點附近的點，且 $\epsilon \to 0$ ，泛函對 $\epsilon$ 求導得到的就是 $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}{dx^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 \quad (j = 1, 2, \dots, n)$
對拉格朗日乘數 $\lambda_i$ 求偏導并令其為零（自動滿足約束條件）：
$\frac{\partial L}{\partial \lambda_i} = -g_i(\mathbf{x}) = 0 \quad (i = 1, 2, \dots, m)$

這組方程（共 $m + n$ 個）稱為KKT條件（對等式約束情形，KKT條件退化為拉格朗日乘數法的條件），其解 $(\mathbf{x}^*, \lambda^*)$ 即為原問題的候選極值點。
若約束條件包含不等式約束 $h_i(\mathbf{x}) \leq 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 \geq 0$ 是對應不等式約束的拉格朗日乘數，且需滿足互補松弛條件：
$\mu_j h_j(\mathbf{x}) = 0 \quad$ (即約束 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))$
系統約束： $\mathbf{\dot{x}} = \mathbf{f}(\mathbf{x}, \mathbf{u}, t)$
首末約束： $\mathbf{x}(t_0) = \mathbf{x_0}$

構造哈密頓函數（Hamiltonian）
引入協態變量（伴隨變量） $\boldsymbol{\lambda}(t) \in \mathbb{R}^n$ （與狀態變量維數相同），定義哈密頓函數： $H(\mathbf{x}, \mathbf{u}, \boldsymbol{\lambda}, t) = L(\mathbf{x}, \mathbf{u}, t) + \boldsymbol{\lambda}^T \mathbf{f}(\mathbf{x}, \mathbf{u}, t)$
- 協態變量 $\boldsymbol{\lambda}(t)$ 的物理意義可理解為狀態變量的“影子價格”或“邊際代價”。
- 如果有狀態約束，還需要再加一項，變成： $H(\mathbf{x}, \mathbf{u}, \boldsymbol{\lambda}, t) = L(\mathbf{x}, \mathbf{u}, t) + \boldsymbol{\lambda}^T \mathbf{f}(\mathbf{x}, \mathbf{u}, t)+\boldsymbol{\eta}^T(\mathbf{x}-\mathbf{x_m})$
正則方程（協態方程與狀態方程）
根據極小值原理，最優解需滿足以下正則方程：

狀態方程（與原系統一致）： $\frac{\partial H}{\partial \boldsymbol{\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 \boldsymbol{\lambda} \right] = \dot{\boldsymbol{\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 + \boldsymbol{\lambda}^T f + \dot{\boldsymbol{\lambda}}^T x \right] dt - \boldsymbol{\lambda}^T(t_f) x(t_f) + \boldsymbol{\lambda}^T(t_0) x(t_0)+ \Phi(x(t_f))$
- 帶入哈密頓函數得： $J_a = \int_{t_0}^{t_f} \left[ H + \dot{\boldsymbol{\lambda}}^T x \right] dt - \boldsymbol{\lambda}^T(t_f) x(t_f) + \boldsymbol{\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 x \, dt + \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 u \, dt = 0$
- 對任意 $\delta x$ ，變分 $\delta J_a = 0$ 需滿足以下條件：
  - 積分項為0： $\frac{\partial H}{\partial x} + \dot{\lambda}^T = 0 \quad \Rightarrow \quad \dot{\lambda} = -\left( \frac{\partial H}{\partial x} \right)^T$ 這就是協態方程
  - 剩余項為0： $-\lambda(t_f) + \frac{\partial \Phi}{\partial x(t_f)} = 0 \quad \Rightarrow \quad \lambda(t_f) = \frac{\partial \Phi}{\partial x(t_f)}$
- 對任意 $\delta u$ ，變分 $\delta J_a = 0$ 需滿足以下條件：
  - $\frac{\partial H}{\partial u} = 0 \quad \text{（無約束時）} \quad \text{或} \quad u^* = \arg \min_{u \in U} H \quad \text{（有約束時）}$ 即：無約束時直接求導為0，有約束時必須是最優控制曲線， $u^*$ 為函數空間中的 “駐點”，即沿任何變分方向的泛函變化率為零。

終端條件
根據終端狀態是否固定，終端條件分為：
- 終端狀態固定（ $\mathbf{x}(t_f) = \mathbf{x}_f$ 已知）： $\boldsymbol{\lambda}(t_f) = \frac{\partial \phi}{\partial \mathbf{x}(t_f)}$
  - 推導見上方協態方程部分
- 終端狀態自由：
  $\boldsymbol{\lambda}(t_f) = \frac{\partial \phi}{\partial \mathbf{x}(t_f)}, \quad H(t_f) = -\frac{\partial \phi}{\partial t_f} \quad (\text{若 } t_f \text{ 自由})$
極小值條件（核心條件）
對于所有 $\in [t_0, t_f]$ 和 $\mathbf{u}(t) \in \mathbb{U}$ ，最優控制 $\mathbf{u}^*(t)$ 需滿足： $H(\mathbf{x}^*(t), \mathbf{u}^*(t), \boldsymbol{\lambda}^*(t), t) = \min_{\mathbf{u} \in \mathbb{U}} H(\mathbf{x}^*(t), \mathbf{u}, \boldsymbol{\lambda}^*(t), t)$
即哈密頓函數在最優控制處取得極小值（若目標為最大化 $J$ ，則取極大值）。這一條件體現了控制約束 $\mathbf{u} \in \mathbb{U}$ 的作用，當控制無約束時，可通過求導 $\frac{\partial H}{\partial \mathbf{u}} = \mathbf{0}$ 直接求解。

步驟總結

構造哈密頓函數 $H$ ，包含狀態、控制、伴隨變量和時間。
建立正則方程：

解狀態方程 $\dot{\mathbf{x}} = \frac{\partial H}{\partial \boldsymbol{\lambda}}$ ，
解協態方程 $\dot{\boldsymbol{\lambda}} = -\frac{\partial H}{\partial \mathbf{x}}$ 。

確定終端條件：根據終端狀態和時間是否固定，設定 $\boldsymbol{\lambda}(t_f)$ 和 $t_f$ 的條件。
求解極小值條件：利用 $H$ 對 $\mathbf{u}$ 的極小化（或極大化），結合控制約束 $\mathbb{U}$ ，求出 $\mathbf{u}^*(t)$ 。
聯立方程求解：將 $\mathbf{u}^*(t)$ 代入正則方程，解兩點邊值問題（BVP）得到 $\mathbf{x}^*(t)$ 和 $\boldsymbol{\lambda}^*(t)$ 。

簡單示例：最速降線問題

以下將通過一個典型的雙積分器系統最短時間控制問題，詳細講解如何運用龐特里亞金極小值原理求解最優控制問題。我們將從問題建模、哈密頓函數構造、協態方程推導、控制律求解到最優軌跡分析逐步展開，并使用$$環境呈現公式。

問題設定：雙積分器系統的最短時間控制
考慮如下二階線性系統（雙積分器）： $\begin{cases} \dot{x}_1 = x_2, \\ \dot{x}_2 = u, \end{cases}$
其中：

$x_1(t), x_2(t)$ 為狀態變量（位置和速度），
$u (t)$ 為控制輸入，滿足幅值約束 $\leq 1$ ，
目標：從初始狀態 $x_{10}, x_{20})$ 轉移到終端狀態 $(0, 0)$ ，且時間最短。

系統狀態方程建模為：
$\begin{bmatrix} \dot{x_1}\\\dot{x_2} \end{bmatrix}=\begin{bmatrix}0 & 1\\0 & 0\end{bmatrix}\begin{bmatrix}x_1 \\x_2\end{bmatrix}+\begin{bmatrix}0 \\1\end{bmatrix}u$

性能指標為：
$\min_u \int_0^T 1 \, dt = T \quad \text{（最小化時間）}.$

第一步：構造哈密頓函數
根據龐特里亞金極小值原理，定義哈密頓函數 $H$ 為：
$\begin{aligned} H &= L + \boldsymbol{\lambda}^T\dot{x} \\ &= \lambda_1 \dot{x}_1 + \lambda_2 \dot{x}_2 + L \end{aligned}$
其中：

$\lambda_1(t), \lambda_2(t)$ 為協態變量（伴隨變量），
$L$ 為性能指標的被積函數（此處 $L = 1$ ）。

代入狀態方程 $\dot{x}_1 = x_2$ 和 $\dot{x}_2 = u$ ，得：
$\lambda_1 x_2 + \lambda_2 u + 1.$
第二步：推導協態方程
根據PMP，協態變量的導數滿足：
$\dot{\lambda}_i = -\frac{\partial H}{\partial x_i} \quad (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^*(t) = \arg\min_{|u| \leq 1} H = \arg\min_{|u| \leq 1} \left( \lambda_2 u + \text{其他} \right).$
由于 $\lambda_2$ 是線性函數，最大化 $H$ 等價于取 $u$ 與 $\lambda_2$ 同號。因此：
$u^*(t) = \text{sign}(\lambda_2(t)) = \begin{cases} +1, & \lambda_2(t) < 0, \\ -1, & \lambda_2(t) > 0. \end{cases}$
當 $\lambda_2(t) = 0$ 時，可能存在切換點，此時控制輸入發生翻轉。

第四步：分析協態變量與切換特性
由于 $\lambda_2(t) = -c_1 t + c_2$ 是線性函數，其符號最多改變一次，因此最優控制 $u^*(t)$ 最多切換一次，屬于Bang-Bang控制。分兩種情況討論：

無切換（直接Bang）：若 $\lambda_2(t)$ 始終非負或非正，則 $u^*(t)$ 保持 $+ 1$ 或 $? 1$ 不變；
一次切換（Bang-Bang）：若 $\lambda_2(t)$ 由正變負或由負變正，則存在唯一切換時刻 $\tau$ ，使得：
$u^*(t) = \begin{cases} +1, & t < \tau, \\ -1, & t \geq \tau \end{cases} \quad \text{或} \quad \begin{cases} -1, & t < \tau, \\ +1, & t \geq \tau. \end{cases}$

第五步：利用終端條件求解常數
終端狀態固定為 $(0, 0)$ ，終端時間 $T$ 自由。根據PMP，終端時刻的橫截條件為： $H (T) = 0.$ 代入 $H$ 的表達式：
$\lambda_1 x_2(T) + \lambda_2(T) u(T) + 1 = 0.$
由于 $x_2(T) = 0$ （終端速度為0），化簡得：
$\lambda_2(T) u(T) + 1 = 0.$
又因 $\text{sign}(\lambda_2(T))$ ，故：
$|\lambda_2(T)| = 1 \quad \Rightarrow \quad \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$ ，則：

當 $\in [0, \tau]$ 時， $u = + 1$ ，狀態方程積分得：
$\begin{cases} x_1(t) = x_{10} + x_{20} t + \frac{1}{2} t^2, \\ x_2(t) = x_{20} + t. \end{cases}$
當 $\in [\tau, T]$ 時， $u = ? 1$ ，狀態方程積分得：
$\begin{cases} 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{cases}$
終端條件 $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 \leq 0$ ），
$\Gamma^-$ （對應 $u = + 1$ 直接到達原點）： $x_1 = \frac{1}{2} x_2^2$ （當 $x_2 \geq 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*x22（上支）或x1 = -0.5*x22（下支）# 初始狀態在開關曲線上方時（x10 > 0.5*x202），切換策略為 -1 → +1（先制動后加速）if (x20 > 0 and x10 > -0.5 * x20**2) or (x20 < 0 and x10 > 0.5 * x20**2):# 推導方程：τ2 - 2x20τ + (0.5x202 - 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*x202），切換策略為 +1 → -1（先加速后制動）elif (x20 > 0 and x10 < -0.5 * x20**2) or (x20 < 0 and x10 < 0.5 * x20**2):# 推導方程：τ2 + 2x20τ + (x10 + 0.5x202) = 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.5x22）且速度非負時，直接用-1制動T = x20 + np.sqrt(x20**2 - 2*x10)  # 推導自x2(T)=0和x1(T)=0def u(t):return -1.0else:# 下開關曲線（x1=-0.5x22）且速度負時，直接用+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)