非线性规划中的差分动态规划（DDP）算法进阶题
题目描述
考虑一个离散时间的确定性最优控制问题，可表述为一个非线性规划问题。假设系统动力学由非线性状态方程描述：
\[x_{k+1} = f(x_k, u_k), \quad k=0,1,\dots,N-1 \]
其中 \(x_k \in \mathbb{R}^{n_x}\) 是状态向量，\(u_k \in \mathbb{R}^{n_u}\) 是控制向量，\(f: \mathbb{R}^{n_x} \times \mathbb{R}^{n_u} \to \mathbb{R}^{n_x}\) 是二次连续可微的函数。给定初始状态 \(x_0\)。目标是最小化总代价函数：
\[J = \phi(x_N) + \sum_{k=0}^{N-1} L(x_k, u_k) \]
其中 \(\phi: \mathbb{R}^{n_x} \to \mathbb{R}\) 是终端代价函数，\(L: \mathbb{R}^{n_x} \times \mathbb{R}^{n_u} \to \mathbb{R}\) 是阶段代价函数，两者均二次连续可微。控制序列 \(\{u_0, u_1, \dots, u_{N-1}\}\) 是决策变量，无显式约束。
本题要求：使用差分动态规划（DDP） 算法求解此最优控制问题。解释DDP算法的核心思想，推导其关键步骤（包括回溯法和前向法中的二次展开、局部优化、反馈策略更新等），并阐述如何通过迭代改进控制序列。提供一个简单的数值例子（例如，小车倒立摆的摆动问题或简单二次型问题）展示算法步骤。
解题过程详解
差分动态规划（DDP）是一种基于动态规划和二阶近似的迭代方法，用于求解离散时间最优控制问题。其核心思想是：在给定当前控制序列（标称轨迹）附近，对动力学和代价函数进行二阶展开，通过求解局部二次规划（即线性二次型调节器LQR问题）得到改进的控制策略，并沿新轨迹迭代更新。下面循序渐进讲解：
第一步：问题重述与DDP基本原理
-
问题形式：我们有一个 \(N\) 阶段的控制问题，决策变量是整个控制序列 \(U = [u_0, u_1, \dots, u_{N-1}]\)，状态序列 \(X = [x_1, \dots, x_N]\) 由动力学方程决定。这是一个高维无约束非线性规划问题，维度为 \(N \cdot n_u\)。
-
动态规划视角：定义值函数（最优代价-动作函数）：
\[ V_k(x) = \min_{u_k, u_{k+1}, \dots} \left[ \phi(x_N) + \sum_{j=k}^{N-1} L(x_j, u_j) \right] \]
满足 \(x_k = x\) 和系统动力学。根据贝尔曼最优性原理：
\[ V_k(x) = \min_{u} \left[ L(x, u) + V_{k+1}(f(x, u)) \right] \]
- DDP思路：从一个初始猜测控制序列 \(U^{(0)}\) 出发，通过前向模拟得到标称状态轨迹 \(\bar{X}\)。然后，在标称轨迹附近，对上述贝尔曼方程中的右端项进行二阶泰勒展开，得到局部LQR问题，其解给出一个反馈控制律（包括开环修正和线性反馈项）。利用这个反馈律，在下一个迭代中沿轨迹施加改进的控制，并更新标称轨迹。迭代直至收敛。
第二步：DDP算法核心推导
假设当前迭代中，我们有标称轨迹 \(\{\bar{x}_k, \bar{u}_k\}_{k=0}^{N-1}\)，满足 \(\bar{x}_{k+1} = f(\bar{x}_k, \bar{u}_k)\)。
2.1 反向传播（Backward Pass）
从最后一步 \(k=N\) 开始，向后计算到 \(k=0\)。定义状态和控制的偏差：
\[\delta x_k = x_k - \bar{x}_k, \quad \delta u_k = u_k - \bar{u}_k \]
- 展开代价函数：
定义哈密顿函数（仅用于推导）：
\[ H_k(x, u) = L(x, u) + V_{k+1}(f(x, u)) \]
在标称点 \((\bar{x}_k, \bar{u}_k)\) 处，计算其一阶和二阶导数：
- 梯度：
\[ \begin{aligned} q_k &= \nabla_x H_k = L_x + f_x^T V'_{k+1} \\ r_k &= \nabla_u H_k = L_u + f_u^T V'_{k+1} \end{aligned} \]
其中 $ L_x = \nabla_x L(\bar{x}_k, \bar{u}_k) $，$ f_x = D_x f(\bar{x}_k, \bar{u}_k) $ 等，$ V'_{k+1} = \nabla_x V_{k+1}(\bar{x}_{k+1}) $。
- 海森矩阵块：
\[ \begin{aligned} Q_{xx,k} &= H_{xx} = L_{xx} + f_x^T V''_{k+1} f_x + \sum_{i} (V'_{k+1})_i f_{xx}^i \\ Q_{uu,k} &= H_{uu} = L_{uu} + f_u^T V''_{k+1} f_u + \sum_{i} (V'_{k+1})_i f_{uu}^i \\ Q_{ux,k} &= H_{ux} = L_{ux} + f_u^T V''_{k+1} f_x + \sum_{i} (V'_{k+1})_i f_{ux}^i \end{aligned} \]
这里 $ f_{xx}^i $ 是 $ f $ 的第 $ i $ 分量对 $ x $ 的海森矩阵，通常为了计算效率，DDP常忽略二阶动力学项（即 $ f_{xx}, f_{uu}, f_{ux} $），这导出了**iLQR（迭代LQR）** 算法。但经典DDP保留这些项。为简化，以下推导采用iLQR假设（忽略二阶动力学项），但会指出完整DDP的区别。
- 二次逼近的贝尔曼方程：
对 \(H_k\) 在标称点做二阶展开：
\[ \Delta H_k \approx q_k^T \delta x_k + r_k^T \delta u_k + \frac12 \delta x_k^T Q_{xx,k} \delta x_k + \frac12 \delta u_k^T Q_{uu,k} \delta u_k + \delta u_k^T Q_{ux,k} \delta x_k \]
值函数的增量：
\[ V_k(\bar{x}_k + \delta x_k) \approx V_k(\bar{x}_k) + \Delta V_k \]
根据贝尔曼方程，\(V_k(\bar{x}_k + \delta x_k) = \min_{\delta u_k} \left[ H_k(\bar{x}_k + \delta x_k, \bar{u}_k + \delta u_k) \right]\)。
- 求解局部最优控制：
对 \(\delta u_k\) 最小化 \(\Delta H_k\)（一个二次函数）：- 最优条件：\(\frac{\partial \Delta H_k}{\partial \delta u_k} = r_k + Q_{uu,k} \delta u_k + Q_{ux,k} \delta x_k = 0\)
- 解得最优控制修正：
\[ \delta u_k^* = -Q_{uu,k}^{-1} (r_k + Q_{ux,k} \delta x_k) = k_k + K_k \delta x_k \]
其中：
\[ k_k = -Q_{uu,k}^{-1} r_k, \quad K_k = -Q_{uu,k}^{-1} Q_{ux,k} \]
$ k_k $ 是开环修正项，$ K_k $ 是反馈增益矩阵。
- 更新值函数导数的递归：
将 \(\delta u_k^*\) 代回 \(\Delta H_k\) 表达式，得到 \(\Delta V_k\) 的二次型：
\[ \Delta V_k = \frac12 \delta x_k^T S_k \delta x_k + s_k^T \delta x_k + \text{常数} \]
通过匹配系数，得到值函数梯度 \(V'_k\) 和海森 \(V''_k\) 的递推（在 \(\delta x_k=0\) 处）：
\[ \begin{aligned} V'_k &= q_k + K_k^T r_k + K_k^T Q_{uu,k} k_k + Q_{ux,k}^T k_k \\ V''_k &= Q_{xx,k} + K_k^T Q_{uu,k} K_k + Q_{ux,k}^T K_k + K_k^T Q_{ux,k} \end{aligned} \]
实际上，更紧凑的形式（经过代数整理，并利用 \(k_k, K_k\) 定义）是：
\[ \begin{aligned} V'_k &= q_k - K_k^T Q_{uu,k} k_k + Q_{ux,k}^T k_k \\ V''_k &= (Q_{xx,k} - Q_{ux,k}^T Q_{uu,k}^{-1} Q_{ux,k}) + K_k^T Q_{uu,k} K_k \end{aligned} \]
常用计算步骤：
- 首先计算 \(Q_{uu,k}^{-1}\)（需确保正定性，否则正则化）。
- 然后计算 \(k_k, K_k\)。
- 更新：
\[ \begin{aligned} V'_k &\leftarrow q_k + K_k^T r_k + K_k^T Q_{uu,k} k_k + Q_{ux,k}^T k_k \\ V''_k &\leftarrow Q_{xx,k} + K_k^T Q_{uu,k} K_k + Q_{ux,k}^T K_k + K_k^T Q_{ux,k} \end{aligned} \]
初始条件：\(V'_N = \nabla \phi(\bar{x}_N)\), \(V''_N = \nabla^2 \phi(\bar{x}_N)\)。
2.2 前向滚动（Forward Pass）
从 \(k=0\) 到 \(N-1\)，使用反向传播计算出的反馈律，沿新状态轨迹更新控制：
- 给定初始状态 \(x_0^{(new)} = x_0\)。
- 对于 \(k=0,1,\dots,N-1\)：
\[ \begin{aligned} u_k^{(new)} &= \bar{u}_k + k_k + K_k (x_k^{(new)} - \bar{x}_k) \\ x_{k+1}^{(new)} &= f(x_k^{(new)}, u_k^{(new)}) \end{aligned} \]
- 计算新轨迹下的总代价 \(J^{(new)}\)。
第三步：算法流程与迭代改进
-
初始化：猜测初始控制序列 \(U^{(0)}\)，设定正则化参数 \(\mu > 0\)（用于保证 \(Q_{uu}\) 可逆），容忍误差 \(\epsilon\)，最大迭代次数。
-
迭代循环：
a. 前向模拟：用当前控制序列 \(U^{(i)}\) 积分动力学，得到标称状态轨迹 \(\bar{X}\) 和代价 \(J^{(i)}\)。
b. 反向传播：从 \(k=N\) 到 \(0\) 计算 \(q_k, r_k, Q_{xx}, Q_{uu}, Q_{ux}\) 及 \(k_k, K_k\)。如果 \(Q_{uu}\) 不正定，增加 \(\mu\)（即用 \(Q_{uu} + \mu I\) 代替）以保证下降方向。
c. 前向滚动：用反馈律 \(u_k = \bar{u}_k + k_k + K_k (x_k - \bar{x}_k)\) 从 \(x_0\) 生成新轨迹 \(X^{(new)}, U^{(new)}\) 和代价 \(J^{(new)}\)。
d. 接受判断：如果 \(J^{(new)} < J^{(i)}\)，接受新轨迹，更新 \(\bar{X} \leftarrow X^{(new)}\), \(\bar{U} \leftarrow U^{(new)}\)，并减少 \(\mu\)（如 \(\mu \leftarrow \mu/2\)）。否则，增加 \(\mu\) 以更保守，重新反向传播（或减少步长）。
e. 检查收敛：如果代价下降量 \(|J^{(i)} - J^{(new)}| < \epsilon\) 或达到最大迭代，停止。
第四步：简单数值例子
考虑一个简单的二次型问题，以便手工验证：
- 系统：\(x_{k+1} = x_k + u_k\)，标量 \(x, u\)。
- 代价：\(J = \frac12 x_N^2 + \frac12 \sum_{k=0}^{N-1} (x_k^2 + u_k^2)\)。
- 参数：\(N=2, x_0=1\)。
步骤：
- 初始猜测：\(\bar{u}_0 = \bar{u}_1 = 0\)，则 \(\bar{x}_0=1, \bar{x}_1=1, \bar{x}_2=1\)，初始代价 \(J = 0.5*(1+1+1) = 1.5\)。
- 反向传播（iLQR，忽略二阶动力学）：
- \(k=2\)：\(V'_2 = \phi_x = x_2 = 1\), \(V''_2 = 1\)。
- \(k=1\)：\(f_x=1, f_u=1, L_x=x_1=1, L_u=u_1=0, L_{xx}=1, L_{uu}=1, L_{ux}=0\)。
- \(q_1 = L_x + f_x V'_2 = 1 + 1*1 = 2\)
- \(r_1 = L_u + f_u V'_2 = 0 + 1*1 = 1\)
- \(Q_{xx} = L_{xx} + f_x^2 V''_2 = 1 + 1 = 2\)
- \(Q_{uu} = L_{uu} + f_u^2 V''_2 = 1 + 1 = 2\)
- \(Q_{ux} = L_{ux} + f_u f_x V''_2 = 0 + 1 = 1\)
- \(k_1 = -Q_{uu}^{-1} r_1 = -0.5*1 = -0.5\)
- \(K_1 = -Q_{uu}^{-1} Q_{ux} = -0.5*1 = -0.5\)
- 更新 \(V'_1 = q_1 + K_1 r_1 + K_1^2 Q_{uu} + Q_{ux} k_1 = 2 + (-0.5)*1 + 0.25*2 + 1*(-0.5) = 2 -0.5+0.5-0.5=1.5\)
\(V''_1 = Q_{xx} + K_1^2 Q_{uu} + 2 K_1 Q_{ux} = 2 + 0.25*2 + 2*(-0.5)*1 = 2+0.5-1=1.5\)
- \(k=0\)：类似计算（略），得到 \(k_0, K_0\)。
- 前向滚动：用反馈律从 \(x_0=1\) 计算新控制。例如，若得 \(k_0=-0.6, K_0=-0.6\)，则 \(u_0 = 0 + (-0.6) + (-0.6)*(1-1) = -0.6\)，\(x_1 = 1 + (-0.6) = 0.4\)；再用 \(k_1, K_1\) 计算 \(u_1\)，得新轨迹。
- 新代价应下降，迭代至收敛（解析解可通过Riccati方程得到）。
要点小结：
- DDP利用动态规划框架，在标称轨迹附近二次展开，得到局部LQR反馈策略。
- 反向传播计算反馈增益和开环修正，前向滚动执行新控制。
- 通过迭代更新轨迹，通常具有二次收敛速度。
- 与SQP相比，DDP利用问题的时间结构，适合较长时域控制问题。
这个算法进阶题要求理解展开推导、反馈律构造、以及迭代更新机制。在实际应用中，还需处理正则化、步长调整、约束处理（可通过内点法或惩罚函数）等细节。