📁 Sunhao's Log
← 返回笔记列表 / 微分动态规划
微分动态规划 ✍️ Second Brain ★ 读者评价

微分动态规划 (Differential Dynamic Programming, DDP) 全景探秘

A (Long) Peek into Differential Dynamic Programming & Iterative LQR

#optimal-control #trajectory-optimization #ddp #ilqr #robotics #math-heavy #reinforcement-learning

微分动态规划 (Differential Dynamic Programming, DDP) 全景探秘

A (Long) Peek into Differential Dynamic Programming & Iterative LQR
参考 Lilian Weng 经典技术长文范式整理 · 经典控制论 ✕ 现代强化学习核心底座


💡 导语
在波士顿动力双足人形机器人 Atlas 令人惊叹的后空翻、四足机器狗在崎岖路面的奔跑,以及 SpaceX 猎鹰九号火箭垂直回收着陆的背后,有一个共同的“数学发动机”——微分动态规划 (Differential Dynamic Programming, DDP) 及其近代变种 迭代线性二次型调节器 (Iterative LQR, iLQR)。
本文从第一性原理推导出发,彻底打通从经典贝尔曼动态规划、LQR Riccati 差分方程,到局部二阶泰勒展开、反馈/前馈增益解耦、LM 阻尼正规化、Box-QP 约束处理,再到它与现代模型强化学习(Model-Based RL,如 Guided Policy Search)的内在血脉关系。


目录


1. 问题的缘起:从贝尔曼“维度灾难”到局部轨迹优化

理查德·贝尔曼(Richard Bellman)在 1950 年代提出了著名的动态规划(Dynamic Programming, DP)与最优性原理。DP 的宏伟蓝图在于求出整个连续状态空间上的全局最优闭环策略 。然而在实际工程中,当状态维度 超过 4 维时,空间网格化离散所需要的计算复杂度以指数爆炸 攀升,这就是臭名昭著的“维度灾难 (Curse of Dimensionality)”。

面对连续、高维、非线性的机器人多刚体动力学系统(例如一个 12 自由度四足机器人的状态维度高达 24 到 36 维),我们真的需要求解全状态空间每一个角落的最优策略吗?

答案是否定的。机器人在执行具体机动时(如从椅子上站立、越过障碍物、倒车入库),其状态演化往往紧密围绕在某条**局部时空管道(Local Trajectory Tube)**之内。这就催生了控制理论与现代 AI 的两大核心流派:

对比维度全局策略搜索 (Global Policy Search, 如 PPO/SAC)局部轨迹优化 (Local Trajectory Optimization, 如 DDP)
样本需求量极高(需数百万甚至数千万次物理交互采样)极低(通常仅需数十次动力学导数计算即可收敛)
动力学利用把系统视为黑盒环境,不利用精确物理方程深度融合分析动力学(牛顿-欧拉或拉格朗日方程的导数)
泛化能力覆盖全状态域,训练后推理极快,适应广泛初始状态强依赖初值猜测(Initial Guess),容易陷入局部极小

在局部轨迹优化的世界中,通常有三大主流方法:

  1. 单打靶法 (Single Shooting):直接将控制序列 作为决策变量,状态通过动力学前向积分隐式确定。缺点是对非线性长时域系统极度敏感,微小的控制扰动会导致末端状态巨幅发散。
  2. 直接配置法 (Direct Collocation):将状态序列和控制序列同时离散化为决策变量,动力学校验作为等式约束。转化为一个极其庞大但稀疏的大规模非线性规划(NLP,常借助 IPOPT/SNOPT 求解)。优点是易处理复杂约束,缺点是每次迭代需分解巨型稀疏 KKT 矩阵,且不直接生成闭环反馈增益。
  3. 微分动态规划 (DDP / iLQR):两者的黄金平衡点。它将贝尔曼最优性原理与二阶牛顿优化巧妙结合。沿时间反向递归时,由于贝尔曼马尔可夫性,时间复杂度关于时域长度 保持严格线性 (而非稠密非线性规划的 )。更关键的是,它在反向递推的同时天然输出了时变局部反馈增益矩阵 ,赋予了生成的轨迹对抗扰动的鲁棒跟踪能力!

2. 轨迹优化数学表述与 Q 算子

考虑离散时间有限时域(Finite-Horizon)非线性动态系统:

其中状态向量 ,控制输入向量 , 为可微非线性动力学映射。

系统在全时域内的累积标量代价值(Total Cost Functional)定义为各步运行代价(Running Cost)与最终终端代价(Terminal Cost)之和:

轨迹优化的目标是:给定已知初始状态 ,寻找一个最优控制输入序列 ,使得总体代价 达到极小:

价值函数 (Value Function / Cost-to-Go)

根据动态规划思想,我们定义从第 步状态 开始到时域结束所能取得的最优剩余代价为价值函数 :

在终止时刻 ,边界条件为:

贝尔曼最优性原理与连续空间的 Q 算子

贝尔曼最优性逆向递归方程式为:

定义连续状态-动作价值函数(Action-Value / Q-function):

于是贝尔曼最优性方程可以极简地表达为:


3. 基石演进:有限时域离散时间 LQR

假设动力学为线性时变系统,代价函数为严格二次型:

其中权矩阵满足 (半正定),(正定),。

Riccati 差分方程推导

在时刻 ,。猜想对任意时刻 ,价值函数均能保持为纯二次型结构:

代入第 步的 函数:

Q(x_t, u_t) &= \frac{1}{2} x_t^T Q_t x_t + \frac{1}{2} u_t^T R_t u_t + \frac{1}{2} (A_t x_t + B_t u_t)^T P_{t+1} (A_t x_t + B_t u_t) \\ &= \frac{1}{2} x_t^T \big( Q_t + A_t^T P_{t+1} A_t \big) x_t + \frac{1}{2} u_t^T \big( R_t + B_t^T P_{t+1} B_t \big) u_t + u_t^T \big( B_t^T P_{t+1} A_t \big) x_t \end{aligned}$$ 对控制输入 $u_t$ 求偏导并令其为 0: $$\nabla_{u_t} Q = (R_t + B_t^T P_{t+1} B_t) u_t + B_t^T P_{t+1} A_t x_t = 0$$ 解析解得**最优控制律为纯线性的状态反馈控制器**: $$u_t^* = - K_t x_t, \quad K_t = \big( R_t + B_t^T P_{t+1} B_t \big)^{-1} B_t^T P_{t+1} A_t$$ 代回可得著名的**离散时间 Riccati 差分方程 (Discrete-time Riccati Difference Equation)**: $$P_t = Q_t + A_t^T P_{t+1} A_t - A_t^T P_{t+1} B_t \big( R_t + B_t^T P_{t+1} B_t \big)^{-1} B_t^T P_{t+1} A_t$$ --- ## 4. 微分动态规划 (Full DDP) 的完整数学推导 Mayne 与 Jacobson 在 1966-1970 年间创立的 **DDP**,其核心思想是:**在当前已有标称轨迹周围建立局部二阶泰勒展开,并以此驱动贝尔曼反向更新**。 ### 4.1 扰动展开与二阶张量项 设系统当前拥有一条先验的标称轨迹 $(\bar{x}_t, \bar{u}_t)$。定义扰动变量: $$\delta x_t \triangleq x_t - \bar{x}_t, \quad \delta u_t \triangleq u_t - \bar{u}_t$$ 价值函数局部展开至二阶: $$V(\bar{x}_t + \delta x_t, t) \approx V(\bar{x}_t, t) + V_x^T \delta x_t + \frac{1}{2} \delta x_t^T V_{xx} \delta x_t$$ 将动作-价值函数 $Q(x, u)$ 进行二阶多元泰勒展开: $$Q(\bar{x}_t + \delta x_t, \bar{u}_t + \delta u_t) \approx Q(\bar{x}_t, \bar{u}_t) + \begin{bmatrix} Q_x \\ Q_u \end{bmatrix}^T \begin{bmatrix} \delta x_t \\ \delta u_t \end{bmatrix} + \frac{1}{2} \begin{bmatrix} \delta x_t \\ \delta u_t \end{bmatrix}^T \begin{bmatrix} Q_{xx} & Q_{xu} \\ Q_{ux} & Q_{uu} \end{bmatrix} \begin{bmatrix} \delta x_t \\ \delta u_t \end{bmatrix}$$ 利用链式法则求导(上标撇号表示在 $t+1$ 步处评估): **一阶偏导数(梯度项):** $$Q_x = l_x + f_x^T V_x'$$ $$Q_u = l_u + f_u^T V_x'$$ **二阶偏导数(包含动力学二阶张量项):** $$Q_{xx} = l_{xx} + f_x^T V_{xx}' f_x + \sum_{i=1}^n V_{x, i}' \cdot \nabla_{xx}^2 f^i$$ $$Q_{uu} = l_{uu} + f_u^T V_{xx}' f_u + \sum_{i=1}^n V_{x, i}' \cdot \nabla_{uu}^2 f^i$$ $$Q_{ux} = l_{ux} + f_u^T V_{xx}' f_x + \sum_{i=1}^n V_{x, i}' \cdot \nabla_{ux}^2 f^i \quad (Q_{xu} = Q_{ux}^T)$$ > *注:$\sum_{i=1}^n V_{x, i}' \cdot \nabla^2 f^i$ 常记作张量缩并积 $V_x' \cdot f_{xx}$,度量了动力学非线性曲率。* ### 4.2 反向传播 (Backward Pass) 对局部二次模型极小化: $$\nabla_{\delta u_t} Q = Q_u + Q_{uu} \delta u_t + Q_{ux} \delta x_t = 0$$ 若 $Q_{uu} \succ 0$,解得**最优控制扰动解**: $$\delta u_t^*(\delta x_t) = k_t + K_t \delta x_t$$ 其中: - **开环前馈步长 (Feedforward Gain)**:$k_t = - Q_{uu}^{-1} Q_u$ - **闭环状态反馈增益 (Feedback Gain)**:$K_t = - Q_{uu}^{-1} Q_{ux}$ 期望单步代价缩减量: $$\Delta V_t = -\frac{1}{2} Q_u^T Q_{uu}^{-1} Q_u$$ 将最优控制律代回,价值函数的一阶度与二阶度逆向递推: $$V_x(t) = Q_x - K_t^T Q_{uu} k_t = Q_x + K_t^T Q_u$$ $$V_{xx}(t) = Q_{xx} - K_t^T Q_{uu} K_t$$ 反向传播初始化边界条件: $$V_x(N) = \nabla_x l_f(\bar{x}_N), \qquad V_{xx}(N) = \nabla_{xx}^2 l_f(\bar{x}_N)$$ ### 4.3 前向仿真回放 (Forward Pass) 在完成反向传播后,在真实非线性动力学中应用带线搜索因子 $\alpha \in (0, 1]$ 的前向回放: $$\hat{x}_0 = x_0$$ $$\hat{u}_t = \bar{u}_t + \alpha k_t + K_t (\hat{x}_t - \bar{x}_t)$$ $$\hat{x}_{t+1} = f(\hat{x}_t, \hat{u}_t)$$ 采用 **Armijo 回溯准则** 校验充分下降: $$J(\hat{\mathbf{X}}, \hat{\mathbf{U}}) - J(\bar{\mathbf{X}}, \bar{\mathbf{U}}) \le c_1 \left( \alpha \sum_{t=0}^{N-1} k_t^T Q_u + \frac{1}{2} \alpha^2 \sum_{t=0}^{N-1} k_t^T Q_{uu} k_t \right)$$ --- ## 5. 迭代 LQR (iLQR / SLQ) 与高斯-牛顿近似 在多刚体机器人中,计算三阶张量 $f_{xx} \in \mathbb{R}^{n \times n \times n}$ 代价极其高昂。 **高斯-牛顿近似(Gauss-Newton Approximation)** 忽略动力学的所有二阶项($f_{xx} \approx 0, f_{uu} \approx 0, f_{ux} \approx 0$): $$Q_{xx} \approx l_{xx} + f_x^T V_{xx}' f_x$$ $$Q_{uu} \approx l_{uu} + f_u^T V_{xx}' f_u$$ $$Q_{ux} \approx l_{ux} + f_u^T V_{xx}' f_x$$ **对比速览:** - **Full DDP**:二阶牛顿法,二次收敛速度,计算复杂,对非凸 Hessian 敏感; - **iLQR**:高斯-牛顿法,超线性收敛,仅需一阶雅可比 $f_x, f_u$,速度快 3~10 倍,结构性半正定,是工业界主流。 --- ## 6. 工程落地与数值避坑技巧 ### 6.1 LM 阻尼与正规化 为确保 $Q_{uu} \succ 0$: - **控制空间阻尼**:$\tilde{Q}_{uu} = Q_{uu} + \mu I_m$($\mu \to \infty$ 时退化为最速梯度下降,$\mu \to 0$ 时为牛顿步); - **状态空间阻尼**:$\tilde{V}_{xx}' = V_{xx}' + \mu I_n$。 ### 6.2 控制输入约束与 Box-QP 电机存在物理扭矩限制 $u_{\min} \le u_t \le u_{\max}$。 Yuval Tassa 等人(IROS 2012 / 2014)提出 **Box-QP 投影牛顿法**,将反向每步建模为带上下界的 QP 问题,仅在自由变量空间求逆,3~5 次内即可求得精确解析解。 ### 6.3 不可行初值与多重打靶 FDDP 传统 DDP 要求初始轨迹动力学严格可行。 Carlos Mastalli 等人(ICRA 2020, Crocoddyl)提出 **FDDP (Feasibility-driven DDP)**,允许初始残差 $\bar{d}_t = f(\bar{x}_t, \bar{u}_t) - \bar{x}_{t+1} \neq 0$,并在迭代中闭合间隙,极大增强在复杂接触动作中的初值鲁棒性。 --- ## 7. 控制论与现代强化学习的交汇 微分动态规划并非与现代强化学习割裂的孤岛,二者在底层哲学与算法架构上深度交融。 ```mermaid flowchart TD Sensors["📷 高维原始感知观测 (RGB 相机图像 / 激光雷达点云 $o_t$)"] --> Teacher subgraph GPS ["🧭 引导策略搜索 (Guided Policy Search, GPS) 双层闭环架构"] direction TB Teacher["👨‍🏫 导师分支 (Teacher: DDP / iLQR 局部轨迹优化簇)<br/>轨迹 1: $(x_t^1, u_t^1, K_t^1)$<br/>轨迹 2: $(x_t^2, u_t^2, K_t^2)$<br/>轨迹 M: $(x_t^M, u_t^M, K_t^M)$"] Distill["🔄 监督蒸馏与策略对齐 (KL 散度 / ADMM 约束: $D_{KL}(p_i \parallel \pi_\\theta) \le \epsilon$)"] Student["🎓 学生分支 (Student: 深度神经网络全局策略 $\pi_\\theta(a|o)$)<br/>高维感知端到端直接映射至连续控制动作"] Teacher --> Distill --> Student end Student --> Output["⚡ 实时高频端到端毫秒级控制推理 (500Hz ~ 1kHz Real-time Control)"] ``` 1. **Guided Policy Search (GPS, Sergey Levine)**: 在多个局部初始点使用 iLQR 优化出优质轨迹与时变反馈增益序列,利用监督蒸馏(KL 散度 / ADMM 约束)训练单一深度神经网络全局策略 $\pi_\theta(a|o)$,同时兼顾了局部轨迹高精度优化与高维原始视觉感知的全局抗扰泛化能力。 2. **滚动时域模型预测控制 (DDP-MPC)**: 以 100Hz 运行局部短时域 iLQR 优化,只执行第 1 步控制量,并在下一时刻以最新传感器状态热启动(Warm-start)迭代,实现实时对抗外力冲击的动态自愈闭环。 --- ## 8. 完整 Python 算法架构实战 ```python import numpy as np class IterativeLQR: def __init__(self, dynamics, cost, n_x, n_u, horizon): self.dynamics = dynamics self.cost = cost self.n_x = n_x self.n_u = n_u self.N = horizon self.mu = 1e-4 self.mu_min = 1e-6 self.mu_max = 1e8 self.delta_mu = 2.0 def solve(self, x0, u_guess, max_iters=50, tol=1e-5): x_bar = np.zeros((self.N + 1, self.n_x)) u_bar = np.array(u_guess) x_bar[0] = x0 for t in range(self.N): x_bar[t+1] = self.dynamics.f(x_bar[t], u_bar[t]) current_cost = self._calc_total_cost(x_bar, u_bar) for iteration in range(max_iters): k_seq, K_seq, expected_red, success = self._backward_pass(x_bar, u_bar) if not success: self.mu = min(self.mu_max, self.mu * self.delta_mu) continue x_new, u_new, new_cost, accepted = self._forward_pass( x_bar, u_bar, k_seq, K_seq, current_cost, expected_red ) if accepted: cost_diff = current_cost - new_cost x_bar, u_bar, current_cost = x_new, u_new, new_cost self.mu = max(self.mu_min, self.mu / self.delta_mu) if cost_diff < tol: print(f"[iLQR] 收敛于第 {iteration+1} 次迭代,代价: {current_cost:.5f}") break else: self.mu = min(self.mu_max, self.mu * self.delta_mu) return x_bar, u_bar, K_seq def _backward_pass(self, x_bar, u_bar): k_seq = np.zeros((self.N, self.n_u)) K_seq = np.zeros((self.N, self.n_u, self.n_x)) expected_reduction = 0.0 V_x = self.cost.lfx(x_bar[self.N]) V_xx = self.cost.lfxx(x_bar[self.N]) for t in reversed(range(self.N)): xt, ut = x_bar[t], u_bar[t] fx, fu = self.dynamics.fx(xt, ut), self.dynamics.fu(xt, ut) Q_x = self.cost.lx(xt, ut) + fx.T @ V_x Q_u = self.cost.lu(xt, ut) + fu.T @ V_x Q_xx = self.cost.lxx(xt, ut) + fx.T @ V_xx @ fx Q_uu = self.cost.luu(xt, ut) + fu.T @ V_xx @ fu Q_ux = self.cost.lux(xt, ut) + fu.T @ V_xx @ fx Q_uu_damped = Q_uu + self.mu * np.eye(self.n_u) try: L = np.linalg.cholesky(Q_uu_damped) except np.linalg.LinAlgError: return None, None, 0.0, False k = -np.linalg.solve(L.T, np.linalg.solve(L, Q_u)) K = -np.linalg.solve(L.T, np.linalg.solve(L, Q_ux)) expected_reduction += (-k.T @ Q_u - 0.5 * k.T @ Q_uu @ k) V_x = Q_x + K.T @ Q_u + K.T @ Q_uu @ k + Q_ux.T @ k V_xx = Q_xx + K.T @ Q_uu @ K + K.T @ Q_ux + Q_ux.T @ K V_xx = 0.5 * (V_xx + V_xx.T) k_seq[t] = k K_seq[t] = K return k_seq, K_seq, expected_reduction, True def _forward_pass(self, x_bar, u_bar, k_seq, K_seq, current_cost, expected_reduction): alphas = [1.0, 0.5, 0.25, 0.125, 0.0625, 0.0] c1 = 1e-4 for alpha in alphas: x_new = np.zeros_like(x_bar) u_new = np.zeros_like(u_bar) x_new[0] = x_bar[0] for t in range(self.N): u_new[t] = u_bar[t] + alpha * k_seq[t] + K_seq[t] @ (x_new[t] - x_bar[t]) x_new[t+1] = self.dynamics.f(x_new[t], u_new[t]) new_cost = self._calc_total_cost(x_new, u_new) if current_cost - new_cost >= c1 * alpha * expected_reduction: return x_new, u_new, new_cost, True return None, None, current_cost, False def _calc_total_cost(self, x, u): total = self.cost.lf(x[self.N]) for t in range(self.N): total += self.cost.l(x[t], u[t]) return total ``` --- ## 9. 轨迹优化方法全景对比速查 | 算法名称 | 数学本质 | 时间复杂度 | 收敛速度 | 是否自带反馈增益 | 约束适应能力 | 工业级代表开源库 | | :--- | :--- | :--- | :--- | :--- | :--- | :--- | | **LQR** | 解析线性二次最优 | $\mathcal{O}(N)$ | 单步解析完成 | **是**(天然伴生 $K_t$) | 弱(无硬约束) | SciPy / Control-Toolbox | | **iLQR / SLQ** | 高斯-牛顿局部动态规划 | $\mathcal{O}(N (n+m)^3)$ | 超线性 (Superlinear) | **是**(开环前馈+闭环增益) | 良(支持 Box-QP) | Crocoddyl / TrajOpt | | **Full DDP** | 二阶牛顿动态规划 | $\mathcal{O}(N (n^3 + n^2 m))$ | 严格二次 (Quadratic) | **是**(天然伴生 $K_t$) | 良(投影/增广拉格朗日) | Crocoddyl / DDP-Matlab | | **直接配置法 (DIRCOL)** | 稀疏大规模 NLP (SQP) | 稀疏分解 $\mathcal{O}(N)$ | 二次收敛 (IPOPT) | **否**(仅产出开环轨迹) | **极优**(支持任意非线性隐式约束) | Drake / CasADi / SNOPT | | **无模型强化学习 (SAC/PPO)** | 黑盒随机策略梯度搜索 | 样本依赖 | 缓慢经验累积 | **是**(神经网络闭环输出) | 差(依赖惩罚项塑形) | Stable-Baselines3 / RLlib | --- ## 10. 经典文献与进阶阅读指南 1. **Mayne, David Q.** (1966). *"A second-order gradient method for determining optimal trajectories of non-linear discrete-time systems."* International Journal of Control, 3(1), 85-95. (DDP 开山鼻祖) 2. **Jacobson, David H., & Mayne, David Q.** (1970). *"Differential Dynamic Programming."* Elsevier, New York. (奠基巨著) 3. **Li, Weiwei, & Todorov, Emanuel.** (2004). *"Iterative linear quadratic regulator design for nonlinear systems."* ICINCO 2004. (iLQR 奠基论文) 4. **Tassa, Yuval, Erez, Tom, & Todorov, Emanuel.** (2012). *"Synthesis and stabilization of complex behaviors through online trajectory optimization."* IROS 2012. (Box-QP 控制限幅) 5. **Mastalli, Carlos, et al.** (2020). *"Crocoddyl: An efficient and versatile framework for multi-contact optimal control."* ICRA 2020. (FDDP 多重打靶与多接触运控) 6. **Levine, Sergey, & Koltun, Vladlen.** (2013). *"Guided policy search."* ICML 2013. (将 iLQR 优化作为导师引入深度学习全局策略)
🌟 读者学术评价与反馈 • PEER REVIEW & RATING

为本文打星评价

如果您觉得本文公式推导严谨、复现内容有启发,欢迎为文章打分支持!

0.0
☆☆☆☆☆
0 人已评价
👈 点击星星即可完成打分(1~5星)