动态规划与最优控制理论及应用详解
简介：《动态规划与最优控制》是Bertsekas教授的经典教材，系统介绍了动态规划的基本理论与最优控制方法。作为优化理论与决策过程的重要工具，动态规划被广泛应用于工程、经济、计算机科学等多个领域。本书深入浅出地讲解了状态空间、Bellman方程、离散与连续时间动态规划、线性二次型控制（LQR）、数值求解方法、最优控制原理以及马尔科夫决策过程等内容。通过大量应用实例，如资源调度、机器人路径规划和通信网络优化，帮助读者掌握动态规划在实际问题中的建模与求解技巧，是控制理论与优化领域的重要参考书。
1. 动态规划的基本概念与核心思想
动态规划（Dynamic Programming, DP）是一种用于求解具有重叠子问题和最优子结构特性的多阶段决策问题的数学优化方法。其核心思想在于： 将一个复杂问题分解为一系列相互关联的子问题，并通过递推关系逐阶段求解，最终获得全局最优解 。
在动态规划中， 状态（State） 描述系统在某一阶段的特征， 决策（Decision） 表示在该状态下所采取的行动， 价值函数（Value Function） 则衡量在某一状态下采取最优策略所能获得的累计收益。通过这些基本要素的建模，可以构建出描述问题的动态系统。
最优性原理（Principle of Optimality） 是动态规划的理论基石，它指出： 无论初始状态和初始决策如何，剩余决策必须构成最优策略 。这一原理确保了动态规划可以通过递归或递推方式进行求解。
2. Bellman方程原理与动态规划建模
Bellman方程是动态规划理论的核心，它通过最优性原理（Principle of Optimality）建立了当前状态与后续状态之间的递推关系。本章将深入探讨Bellman方程的数学形式、推导逻辑以及在建模中的实际应用。我们将从基本形式入手，逐步过渡到离散时间系统下的建模方法，并通过经典案例（如库存控制与资源分配）来展示其应用价值。
2.1 Bellman方程的基本形式
Bellman方程的核心在于将最优策略分解为当前决策与未来最优决策的组合。它不仅适用于确定性系统，也广泛用于随机系统（如马尔科夫决策过程）。在本节中，我们将介绍Bellman方程的基本形式，并讨论其与最优性原理的关系。
2.1.1 最优性原理与递推关系
最优性原理由Richard Bellman提出，其核心思想是： 一个最优策略的任何子策略也是该子问题的最优策略 。这一原理允许我们将复杂问题分解为多个子问题，从而形成递推结构。
例如，假设我们有一个有限阶段决策问题，阶段数为 $ T $，每个阶段的状态为 $ s_t $，动作为 $ a_t $，即时成本为 $ c(s_t, a_t) $，状态转移函数为 $ s_{t+1} = f(s_t, a_t) $。那么最优性原理可以表示为：
V_t(s_t) = \min_{a_t} \left[ c(s_t, a_t) + V_{t+1}(f(s_t, a_t)) \right]
其中 $ V_t(s_t) $ 表示从状态 $ s_t $ 出发，从阶段 $ t $ 到 $ T $ 的最小总代价。该方程表明：当前阶段的最优价值函数等于当前动作的代价加上下一阶段的最优价值函数。
代码示例： 以下是一个简单的Python代码片段，用于计算有限阶段下的Bellman递推：
def bellman_recursion(costs, transitions, T):
V = [{} for _ in range(T+1)] # 存储每个阶段的价值函数
# 初始化终值函数
for s in states[T]:
V[T][s] = terminal_cost(s)
# 从 T-1 阶段逆推
for t in range(T-1, -1, -1):
for s in states[t]:
min_cost = float('inf')
for a in actions[t][s]:
next_s = transitions[t][s][a]
current_cost = costs[t][s][a] + V[t+1][next_s]
if current_cost < min_cost:
min_cost = current_cost
V[t][s] = min_cost
return V
代码分析：
-
costs：表示各阶段的即时代价函数。 -
transitions：状态转移函数，给出每个状态动作对的下一个状态。 -
T：阶段总数。 -
V：每个阶段的价值函数，存储为字典形式以适应状态空间。 - 代码从最后阶段开始逆推，利用Bellman递推公式计算每个阶段的最优价值函数。
2.1.2 状态转移函数与价值函数的定义
状态转移函数 $ f(s_t, a_t) $ 描述了系统在状态 $ s_t $ 下执行动作 $ a_t $ 后所进入的新状态 $ s_{t+1} $。对于确定性系统，它是一个确定性映射；对于随机系统，则是一个概率分布函数。
价值函数 $ V(s_t) $ 是从状态 $ s_t $ 开始，执行最优策略后所能获得的最小总回报（或最大总奖励）。它定义为：
V(s_t) = \min_{\pi} \mathbb{E} \left[ \sum_{k=t}^{T} \gamma^{k-t} r(s_k, a_k) \right]
其中 $ \gamma \in [0,1] $ 是折扣因子，用于衡量未来回报的重要性。
表格：状态转移与价值函数的对比
| 概念 | 确定性系统 | 随机系统（MDP） |
|---|---|---|
| 状态转移函数 | $ s_{t+1} = f(s_t, a_t) $ | $ P(s_{t+1} |
| 价值函数 | $ V(s_t) = \min_{a_t} \left[ c(s_t, a_t) + V(f(s_t, a_t)) \right] $ | $ V(s_t) = \min_{a_t} \left[ \mathbb{E} {s {t+1}}[c(s_t, a_t) + \gamma V(s_{t+1})] \right] $ |
2.2 离散时间动态规划建模
在实际问题中，大多数系统都可以建模为离散时间动态系统。本节将介绍如何构建状态空间、动作空间，并展示动态系统建模中的递归结构，最后通过库存控制与资源分配案例说明其应用。
2.2.1 状态空间与动作空间的表示
在离散时间动态规划中，状态空间 $ \mathcal{S} $ 和动作空间 $ \mathcal{A} $ 是建模的基础。
- 状态空间 ：系统可能处于的所有状态集合。例如，在库存控制问题中，状态可以是当前库存量。
- 动作空间 ：在每个状态下可执行的动作集合。例如，在库存控制中，动作可以是订购的数量。
状态和动作通常可以是离散的或连续的。对于离散系统，状态和动作空间可以表示为有限集合；对于连续系统，则需要进行离散化处理。
mermaid流程图：动态规划建模流程
graph TD
A[问题定义] --> B[定义状态空间]
B --> C[定义动作空间]
C --> D[定义状态转移函数]
D --> E[定义即时代价函数]
E --> F[Bellman方程建模]
F --> G[求解最优策略]
2.2.2 动态系统建模中的递归结构
动态系统建模的核心在于递归结构的设计。我们可以将系统的演化过程视为一个递归函数，其中每个阶段的决策依赖于当前状态和未来状态的价值。
例如，考虑一个简单的库存控制问题：
- 每个阶段 $ t $，库存量为 $ s_t $。
- 决策为订购数量 $ a_t $。
- 下一阶段库存量 $ s_{t+1} = s_t + a_t - d_t $，其中 $ d_t $ 为需求量。
- 即时成本包括库存持有成本和缺货成本。
我们可以将该问题建模为：
V_t(s_t) = \min_{a_t} \left[ h(s_t + a_t - d_t)^+ + p(d_t - s_t - a_t)^+ + V_{t+1}(s_t + a_t - d_t) \right]
其中 $ h $ 为单位库存持有成本，$ p $ 为单位缺货成本，$ x^+ = \max(x, 0) $。
代码示例：库存控制问题中的Bellman方程实现
def inventory_control(V_next, s_t, a_t, d_t, h, p):
holding_cost = h * max(s_t + a_t - d_t, 0)
shortage_cost = p * max(d_t - s_t - a_t, 0)
next_state = s_t + a_t - d_t
return holding_cost + shortage_cost + V_next.get(next_state, 0)
代码分析：
-
V_next：下一阶段的价值函数。 -
s_t：当前库存量。 -
a_t：订购数量。 -
d_t：当前阶段的需求量。 -
h：库存持有成本系数。 -
p：缺货惩罚系数。 - 该函数计算当前动作的总代价，包括持有成本、缺货成本和下一阶段的价值。
2.2.3 建模案例：库存控制与资源分配问题
库存控制建模
库存控制是动态规划的经典应用之一。我们考虑一个有限周期的库存控制问题，目标是通过订购策略使总成本最小。
- 状态：库存量 $ s_t $
- 动作：订购量 $ a_t $
- 状态转移：$ s_{t+1} = s_t + a_t - d_t $
- 成本函数：$ c(s_t, a_t) = h \cdot (s_t + a_t - d_t)^+ + p \cdot (d_t - s_t - a_t)^+ $
资源分配问题
资源分配问题中，我们需要将有限资源分配给多个项目以最大化总收益。例如，给定预算 $ B $，需在 $ n $ 个项目中分配资金 $ x_1, x_2, …, x_n $，使得总收益最大。
- 状态：当前阶段的剩余预算 $ b $
- 动作：在当前项目中投入的资金 $ x_i $
- 状态转移：$ b’ = b - x_i $
- 收益函数：$ R_i(x_i) $
该问题可以通过动态规划建模为：
V_i(b) = \max_{0 \leq x_i \leq b} \left[ R_i(x_i) + V_{i+1}(b - x_i) \right]
2.3 Bellman方程的数学推导与求解方法
Bellman方程的求解是动态规划的关键。本节将介绍价值迭代与策略迭代的基本思想，并分析其收敛性与误差估计。
2.3.1 价值迭代与策略迭代的基本思想
- 价值迭代（Value Iteration） ：通过不断更新价值函数来逼近最优解。其更新公式为：
V_{k+1}(s) = \min_a \left[ c(s,a) + \gamma \sum_{s’} P(s’ | s, a) V_k(s’) \right]
- 策略迭代（Policy Iteration） ：先固定一个策略 $ \pi $，评估其价值函数，再根据价值函数改进策略。其步骤包括：
- 策略评估（Policy Evaluation） ：给定策略 $ \pi $，求解对应的Bellman方程：
V^\pi(s) = c(s, \pi(s)) + \gamma \sum_{s’} P(s’ | s, \pi(s)) V^\pi(s’)
- 策略改进（Policy Improvement） ：根据当前价值函数更新策略：
\pi’(s) = \arg\min_a \left[ c(s,a) + \gamma \sum_{s’} P(s’ | s, a) V^\pi(s’) \right]
代码示例：价值迭代算法实现
def value_iteration(P, c, gamma, theta=1e-6):
V = {s: 0 for s in P}
while True:
delta = 0
for s in P:
v = V[s]
V[s] = min(sum(P[s][a][s_prime] * (c[s][a] + gamma * V[s_prime]) for s_prime in P) for a in P[s])
delta = max(delta, abs(v - V[s]))
if delta < theta:
break
return V
代码分析：
-
P：状态转移概率矩阵。 -
c：即时代价函数。 -
gamma：折扣因子。 -
theta：收敛阈值。 - 算法不断迭代更新价值函数，直到变化小于阈值为止。
2.3.2 收敛性分析与误差估计
价值迭代与策略迭代都具有全局收敛性，但在收敛速度和计算复杂度上有所不同。
| 方法 | 收敛性 | 时间复杂度 | 适用场景 |
|---|---|---|---|
| 价值迭代 | 线性收敛 | $ O( | S |
| 策略迭代 | 有限步收敛 | $ O( | S |
误差估计方面，设 $ V_k $ 为第 $ k $ 次迭代的价值函数，$ V^* $ 为最优价值函数，则有：
| V_k - V^ | \leq \frac{\gamma^k}{1 - \gamma} | V_0 - V^ |
这表明价值迭代的误差随迭代次数指数衰减，收敛速度受折扣因子 $ \gamma $ 影响。
通过本章内容，我们系统地介绍了Bellman方程的基本形式、离散时间动态规划建模方法以及价值迭代与策略迭代的核心思想与实现方式。这些内容为后续章节中连续时间动态规划、LQR控制及近似动态规划打下了坚实的理论基础。
3. 连续时间动态规划与微分动态规划（DDP）
在实际工程与控制系统中，许多问题本质上是连续时间系统，例如飞行器控制、机器人动力学、化工过程控制等。为了在这些系统中实现最优控制，我们需要将动态规划理论推广到连续时间域。本章将重点介绍连续时间系统下的Bellman方程——Hamilton-Jacobi-Bellman（HJB）方程，以及在此基础上发展出的微分动态规划（Differential Dynamic Programming, DDP）算法。我们将从理论推导到数值实现，全面解析连续时间动态规划的核心内容。
3.1 连续时间系统下的Bellman方程
3.1.1 Hamilton-Jacobi-Bellman（HJB）方程的推导
在连续时间系统中，系统的状态演化通常由微分方程描述：
\dot{x}(t) = f(x(t), u(t), t)
其中 $x(t) \in \mathbb{R}^n$ 是状态向量，$u(t) \in \mathbb{R}^m$ 是控制输入。我们的目标是最小化如下形式的代价函数：
J(x(t), u(\cdot)) = \int_{t}^{T} L(x(\tau), u(\tau), \tau) d\tau + \Phi(x(T))
其中 $L$ 是瞬时代价函数，$\Phi$ 是终端代价函数。
基于最优性原理（Principle of Optimality），我们可以定义价值函数（Value Function）为：
V(x, t) = \min_{u(\cdot)} \left{ \int_{t}^{T} L(x(\tau), u(\tau), \tau) d\tau + \Phi(x(T)) \right}
通过对价值函数进行微分并应用动态规划的递推思想，可以推导出著名的HJB方程：
-\frac{\partial V}{\partial t} = \min_{u} \left[ L(x, u, t) + \nabla_x V(x, t)^T f(x, u, t) \right]
这个偏微分方程是连续时间动态规划的核心，它将最优控制问题转化为一个非线性偏微分方程的求解问题。
HJB方程的直观解释
- 左边表示价值函数随时间的负变化率。
- 右边则由两部分组成：当前控制动作的代价 $L$ 和状态演化对价值函数的影响 $\nabla_x V^T f$。
- 通过最小化控制输入 $u$，我们可以得到最优控制律。
HJB方程的求解难点
- HJB方程是一个非线性偏微分方程，在高维空间中难以解析求解。
- 通常需要借助数值方法或近似方法（如DDP）来求解。
3.1.2 HJB方程与最优控制律的关系
一旦我们求解出HJB方程的解 $V(x, t)$，就可以从中提取出最优控制律。最优控制输入 $u^*(x, t)$ 满足：
u^*(x, t) = \arg\min_{u} \left[ L(x, u, t) + \nabla_x V(x, t)^T f(x, u, t) \right]
这意味着在每一个时刻，我们都可以通过求解这个最小化问题来确定当前的最优控制输入。这个控制律是闭环反馈控制，因为它依赖于当前状态 $x$。
最优控制律示例
考虑一个简单的线性系统：
\dot{x} = x + u
代价函数为：
L(x, u) = x^2 + u^2
我们可以构造HJB方程并尝试求解：
-\frac{\partial V}{\partial t} = \min_u \left( x^2 + u^2 + \frac{\partial V}{\partial x}(x + u) \right)
对 $u$ 求导并令其为0：
2u + \frac{\partial V}{\partial x} = 0 \Rightarrow u = -\frac{1}{2} \frac{\partial V}{\partial x}
将其代入HJB方程：
-\frac{\partial V}{\partial t} = x^2 + \left( -\frac{1}{2} \frac{\partial V}{\partial x} \right)^2 + \frac{\partial V}{\partial x} \left( x - \frac{1}{2} \frac{\partial V}{\partial x} \right)
这是一个非线性PDE，需要数值方法求解。
控制律提取的代码实现（Python示例）
import numpy as np
from scipy.optimize import minimize
def cost_L(x, u):
return x**2 + u**2
def system_f(x, u):
return x + u
def value_function_gradient(x, t):
# 假设我们已经有一个近似的价值函数梯度
return 2 * x # 示例梯度
def optimal_control(x, t):
def objective(u):
grad_V = value_function_gradient(x, t)
return cost_L(x, u) + grad_V * system_f(x, u)
res = minimize(objective, x0=0)
return res.x[0]
# 示例调用
x = 1.0
t = 0.0
u_opt = optimal_control(x, t)
print(f"Optimal control input u* at x={x}, t={t}: {u_opt}")
代码逻辑分析
-
cost_L：定义瞬时代价函数 $L(x, u)$。 -
system_f：定义系统动力学 $f(x, u)$。 -
value_function_gradient：模拟价值函数的梯度 $\nabla_x V$。 -
optimal_control：通过最小化目标函数来求解最优控制输入。 -
minimize：使用数值优化方法求解控制输入。
3.2 微分动态规划（DDP）算法
3.2.1 DDP的基本框架与二次近似
微分动态规划（Differential Dynamic Programming, DDP）是一种用于连续时间最优控制问题的数值优化方法。其核心思想是：在当前状态轨迹附近对动态系统和代价函数进行泰勒展开，并通过迭代优化轨迹和控制序列。
DDP的基本步骤
- 初始化 ：选择一个初始控制序列 $u_0(t)$，并计算对应的状态轨迹 $x_0(t)$。
- 前向模拟 ：根据当前控制输入模拟系统状态。
- 反向递推 ：
- 对代价函数和动态方程进行二次/线性近似。
- 计算价值函数的梯度和海森矩阵。
- 推导出控制更新律 $u_{k+1}(t) = u_k(t) + \delta u(t)$。 - 收敛判断 ：若更新量足够小，则停止；否则继续迭代。
二次近似过程
在时间步 $t$，我们对代价函数和动态方程进行二阶泰勒展开：
- 动态近似：$\delta x_{t+1} \approx A_t \delta x_t + B_t \delta u_t$
- 代价函数近似：$\delta J \approx \delta x_t^T Q_t \delta x_t + \delta u_t^T R_t \delta u_t + 2 \delta x_t^T Q_{tx} \delta u_t$
其中 $A_t = \frac{\partial f}{\partial x}$, $B_t = \frac{\partial f}{\partial u}$，$Q_t, R_t$ 是代价函数的二阶导数。
DDP算法流程图（mermaid格式）
graph TD
A[初始化控制序列 u0(t)] --> B[前向模拟状态轨迹 x0(t)]
B --> C[反向递推计算价值函数梯度]
C --> D[构造二次近似模型]
D --> E[求解控制更新律 δu(t)]
E --> F[更新控制序列 u_{k+1}(t) = u_k(t) + δu(t)]
F --> G{收敛判断}
G -- 否 --> B
G -- 是 --> H[输出最优控制序列与轨迹]
3.2.2 状态轨迹优化与控制输入更新
在DDP中，状态轨迹和控制输入的更新是通过以下两个关键步骤完成的：
- 反向递推计算价值函数的梯度与海森矩阵
- 前向更新控制输入与状态轨迹
反向递推过程
在时间步 $t$，我们定义：
- $V(x, t)$：价值函数
- $V_x = \nabla_x V$
- $V_{xx} = \nabla_{xx} V$
根据HJB方程的近似形式，我们可以递推地计算：
V_x(t) = Q_x + A_t^T V_x(t+1)
V_{xx}(t) = Q_{xx} + A_t^T V_{xx}(t+1) A_t
其中 $A_t = \frac{\partial f}{\partial x}$。
控制更新律
在每一步迭代中，我们计算控制输入的增量：
\delta u(t) = -K_t \delta x(t)
其中 $K_t$ 是反馈增益矩阵，可通过以下公式计算：
K_t = (R_{uu} + B_t^T V_{xx}(t+1) B_t)^{-1} (B_t^T V_{xx}(t+1) A_t + Q_{ux})
DDP算法实现（Python伪代码）
import numpy as np
def DDP(x0, T, dynamics, cost, max_iter=100, tol=1e-6):
N = len(x0)
u = np.zeros((T, N)) # 初始控制序列
x = simulate(dynamics, x0, u)
for iter in range(max_iter):
# 前向模拟
x = simulate(dynamics, x0, u)
# 反向递推
V_x = np.zeros((T, N))
V_xx = np.zeros((T, N, N))
for t in reversed(range(T)):
A, B = jacobian(dynamics, x[t], u[t])
Q, R, Qx, Qu, Qxx, Quu, Qux = cost(x[t], u[t])
# 反向更新价值函数梯度
if t == T - 1:
V_x[t] = Qx
V_xx[t] = Qxx
else:
V_x[t] = Qx + A.T @ V_x[t+1]
V_xx[t] = Qxx + A.T @ V_xx[t+1] @ A
# 控制更新增益矩阵
inv_term = np.linalg.inv(R + B.T @ V_xx[t+1] @ B)
K = inv_term @ (B.T @ V_xx[t+1] @ A + Qux.T)
# 控制更新
delta_u = -K @ (x[t] - x_desired[t]) # 假设有一个目标状态 x_desired
u[t] += delta_u
# 收敛判断
if np.linalg.norm(delta_u) < tol:
break
return u, x
def simulate(dynamics, x0, u):
T = len(u)
x = [x0]
for t in range(T):
x.append(dynamics(x[t], u[t]))
return np.array(x)
# 示例调用
def dynamics(x, u):
return x + u # 简单系统模型
def cost(x, u):
Q = 1.0
R = 1.0
Qx = 2 * x
Qu = 2 * u
Qxx = 2 * np.eye(1)
Quu = 2 * np.eye(1)
Qux = np.zeros((1, 1))
return Q*x**2 + R*u**2, Qx, Qu, Qxx, Quu, Qux
x0 = np.array([1.0])
T = 10
u_opt, x_opt = DDP(x0, T, dynamics, cost)
print("Optimal control sequence:")
print(u_opt)
代码逻辑分析
-
DDP函数是主循环，包含前向模拟和反向递推。 -
simulate函数用于前向模拟状态轨迹。 -
dynamics和cost是用户定义的系统模型和代价函数。 - 在每次迭代中，计算状态轨迹、价值函数梯度和海森矩阵。
- 根据反馈增益矩阵 $K$ 更新控制输入。
- 当控制更新量小于阈值时停止迭代。
DDP算法参数说明
| 参数 | 含义 | 示例值 |
|---|---|---|
x0 | 初始状态 | [1.0] |
T | 时间步数 | 10 |
dynamics | 系统动力学函数 | x + u |
cost | 代价函数 | 二次型 |
max_iter | 最大迭代次数 | 100 |
tol | 收敛阈值 | 1e-6 |
3.3 连续时间动态规划的数值实现
3.3.1 网格离散化与差分近似
由于HJB方程和DDP算法在连续时间域中难以直接求解，通常需要将问题离散化。常见的离散化方法包括：
- 时间离散化 ：将时间区间 $[t_0, T]$ 划分为 $N$ 个时间步。
- 状态空间离散化 ：将连续状态空间划分为网格点。
- 有限差分法 ：使用前向、后向或中心差分近似导数。
状态空间网格化示例（Python）
import numpy as np
# 定义状态空间范围与分辨率
x_min, x_max = -5, 5
num_points = 100
x_grid = np.linspace(x_min, x_max, num_points)
# 构建网格
X, U = np.meshgrid(x_grid, x_grid)
# 示例：计算价值函数在网格点上的近似值
def value_function_approx(x):
return x**2 # 假设价值函数为平方函数
V = value_function_approx(X)
# 差分近似导数
dV_dx = np.gradient(V, x_grid, axis=1)
差分近似表格说明
| 方法 | 公式 | 误差阶数 |
|---|---|---|
| 前向差分 | $\frac{f(x+h) - f(x)}{h}$ | $O(h)$ |
| 后向差分 | $\frac{f(x) - f(x-h)}{h}$ | $O(h)$ |
| 中心差分 | $\frac{f(x+h) - f(x-h)}{2h}$ | $O(h^2)$ |
3.3.2 实际应用中的稳定性与收敛性处理
在数值实现过程中，HJB方程和DDP算法可能会遇到以下问题：
- 数值不稳定 ：由于近似误差积累，可能导致算法发散。
- 收敛性问题 ：在高维状态空间中，收敛速度变慢甚至不收敛。
- 局部极小问题 ：DDP等优化方法可能陷入局部极小值。
解决方案
- 正则化 ：在优化目标中加入正则项，防止控制输入剧烈变化。
- 信任域方法 ：限制每次更新的步长，确保迭代过程稳定。
- 多起点优化 ：从多个初始点出发，提高找到全局最优的概率。
- 自适应网格划分 ：根据状态轨迹的分布动态调整网格密度。
稳定性增强的DDP代码改进（Python）
def regularized_DDP(...):
...
# 在控制更新中加入正则项
reg_term = 1e-4 * np.eye(len(u[t]))
inv_term = np.linalg.inv(R + B.T @ V_xx[t+1] @ B + reg_term)
...
通过引入正则化项，可以有效防止控制输入的剧烈波动，提高算法稳定性。
收敛性改进策略表格
| 策略 | 描述 | 优点 | 缺点 |
|---|---|---|---|
| 正则化 | 加入控制变化惩罚项 | 提高稳定性 | 可能牺牲最优性 |
| 信任域 | 限制更新步长 | 防止发散 | 增加计算复杂度 |
| 多起点 | 多个初始点启动 | 提高全局最优概率 | 增加运行时间 |
| 自适应网格 | 动态调整状态空间网格 | 提高精度 | 实现复杂度高 |
通过上述策略的综合应用，可以在实际系统中实现稳定、高效的连续时间动态规划控制。
4. 线性系统与二次成本函数下的LQR控制
在控制理论中，线性二次调节器（Linear Quadratic Regulator，简称 LQR）是一种广泛应用的最优控制方法，特别适用于线性系统和二次代价函数的场景。LQR 控制器通过动态规划原理求解最优反馈控制律，使得系统状态快速收敛到期望值，并且在扰动存在时仍能保持良好的鲁棒性。本章将从数学建模、动态规划解法以及扩展应用三个方面，系统阐述 LQR 控制的核心思想与实现方式。
4.1 LQR问题的数学建模
LQR 控制的核心在于对系统的动态模型和代价函数进行合理的数学建模。通过线性系统模型与二次代价函数的结合，LQR 能够将最优控制问题转化为一个可以通过 Riccati 方程求解的优化问题。
4.1.1 线性系统状态方程与二次代价函数
LQR 控制适用于线性系统，其状态方程通常表示为：
\dot{x}(t) = A x(t) + B u(t)
其中：
- $ x(t) \in \mathbb{R}^n $ 是系统状态向量；
- $ u(t) \in \mathbb{R}^m $ 是控制输入向量；
- $ A \in \mathbb{R}^{n \times n} $ 是状态矩阵；
- $ B \in \mathbb{R}^{n \times m} $ 是控制输入矩阵。
对于离散时间系统，状态方程可表示为：
x_{k+1} = A x_k + B u_k
LQR 的目标是设计一个状态反馈控制律 $ u(t) = -K x(t) $，使得系统状态趋于稳定，并最小化如下形式的二次代价函数（Quadratic Cost Function）：
J = \int_0^\infty \left( x^T(t) Q x(t) + u^T(t) R u(t) \right) dt
其中：
- $ Q \in \mathbb{R}^{n \times n} $ 是状态代价矩阵，通常为对称正定矩阵；
- $ R \in \mathbb{R}^{m \times m} $ 是控制代价矩阵，同样为对称正定矩阵；
- $ x^T Q x $ 表示状态偏离期望值的代价；
- $ u^T R u $ 表示控制输入的代价，用于防止控制信号过大。
代价函数的设计体现了对系统性能的权衡：增大 $ Q $ 表示更注重状态收敛速度，增大 $ R $ 则表示更关注控制输入的能量消耗。
4.1.2 代价函数的结构与优化目标
LQR 控制的优化目标是找到一个状态反馈增益矩阵 $ K $，使得代价函数 $ J $ 最小化。该问题可以被建模为一个无限时间范围的最优控制问题。
代价函数的结构具有良好的数学性质：
- 二次形式保证了代价函数是凸函数，因此存在唯一的最小值；
- 线性系统与二次代价函数的组合使得问题可以通过解析方法求解，而无需依赖复杂的数值优化算法；
- 系统稳定性可以通过选择合适的 $ Q $ 和 $ R $ 来保证。
代价函数的优化目标可以理解为：在尽可能小的控制能量下，使系统状态尽快趋近于零（或某个期望状态）。这种优化思想广泛应用于飞行器姿态控制、机器人运动规划、自动驾驶轨迹跟踪等领域。
为了更好地理解代价函数的作用，我们来看一个简单的例子。
代码示例：LQR代价函数的可视化
import numpy as np
import matplotlib.pyplot as plt
# 定义参数
Q = np.array([[1, 0], [0, 1]])
R = np.array([[0.1]])
# 构造网格
x1 = np.linspace(-5, 5, 100)
x2 = np.linspace(-5, 5, 100)
X1, X2 = np.meshgrid(x1, x2)
# 构造代价函数（假设 u = -Kx，K = [k1, k2]）
# 为了简化，我们固定 u = 0，观察代价函数分布
def cost(x1, x2, Q, R):
x = np.array([x1, x2])
return x.T @ Q @ x + 0 # 假设控制输入为0
# 计算代价函数值
J = np.zeros_like(X1)
for i in range(X1.shape[0]):
for j in range(X1.shape[1]):
J[i, j] = cost(X1[i, j], X2[i, j], Q, R)
# 可视化代价函数
plt.contourf(X1, X2, J, levels=50, cmap='viridis')
plt.colorbar(label='Cost')
plt.xlabel('x1')
plt.ylabel('x2')
plt.title('LQR Cost Function (Q = I, R = 0.1)')
plt.show()
代码逻辑分析：
- 本段代码构造了一个二维状态空间下的代价函数图像，用于展示代价函数的形状；
- 代价函数的值由 $ x^T Q x $ 决定，控制输入 $ u $ 被设为 0；
- 图像显示了状态偏离原点时代价函数的分布情况，中心最低点对应状态 $ x=0 $；
- 调整 $ Q $ 和 $ R $ 的值可以改变代价函数的“陡峭”程度，从而影响控制律的设计。
参数说明：
-
Q：状态代价矩阵，控制状态误差的权重； -
R：控制代价矩阵，控制控制输入的代价； -
x1, x2：状态空间的两个维度； -
X1, X2：状态网格； -
J：计算得到的代价函数值矩阵。
4.2 LQR控制的动态规划解法
LQR 控制问题可以通过动态规划方法求解，其核心在于推导并求解 Riccati 方程。该方程描述了最优控制律的结构，并提供了状态反馈增益矩阵 $ K $ 的解析解。
4.2.1 Riccati方程的推导与求解
连续时间系统下，LQR 控制的最优控制律可通过求解连续时间代数 Riccati 方程（CARE）获得：
A^T P + P A - P B R^{-1} B^T P + Q = 0
其中：
- $ P \in \mathbb{R}^{n \times n} $ 是对称正定矩阵；
- $ A, B, Q, R $ 如前所述；
- $ P $ 的求解决定了反馈增益矩阵 $ K $ 的值；
- 控制律为 $ u(t) = -K x(t) $，其中 $ K = R^{-1} B^T P $。
Riccati 方程是一个非线性矩阵方程，通常通过数值方法进行求解，如 Newton-Raphson 法、Schur 分解法等。
对于离散时间系统，对应的 Riccati 方程为：
P = A^T P A - A^T P B (R + B^T P B)^{-1} B^T P A + Q
该方程的求解同样可通过迭代方法实现。
代码示例：使用 Python 求解连续时间 Riccati 方程
from scipy.linalg import solve_continuous_are
import numpy as np
# 定义系统矩阵
A = np.array([[0, 1],
[0, 0]])
B = np.array([[0],
[1]])
Q = np.eye(2)
R = np.array([[1]])
# 求解 Riccati 方程
P = solve_continuous_are(A, B, Q, R)
# 计算反馈增益矩阵 K
K = np.linalg.inv(R) @ B.T @ P
print("Riccati 矩阵 P：")
print(P)
print("反馈增益矩阵 K：")
print(K)
代码逻辑分析：
- 使用
scipy.linalg.solve_continuous_are函数求解连续时间代数 Riccati 方程； - 系统矩阵 $ A $、$ B $ 定义了一个简单的二阶系统（如质量-弹簧系统）；
- $ Q $ 和 $ R $ 设定代价函数的权重；
- 求解得到的 $ P $ 用于计算反馈增益矩阵 $ K $；
- 得到的 $ K $ 可用于构建控制器 $ u = -Kx $。
参数说明：
-
A：状态矩阵； -
B：控制矩阵； -
Q：状态代价矩阵； -
R：控制代价矩阵； -
P：Riccati 方程的解； -
K：反馈增益矩阵。
4.2.2 状态反馈控制律的设计
在 LQR 控制中，状态反馈控制律的设计是关键环节。反馈增益矩阵 $ K $ 决定了系统响应的速度与稳定性。
控制律形式为：
u(t) = -K x(t)
将该控制律代入系统状态方程，得到闭环系统：
\dot{x}(t) = (A - B K) x(t)
闭环系统的稳定性由矩阵 $ A - B K $ 的特征值决定。如果所有特征值的实部小于零，则系统是渐近稳定的。
mermaid流程图：LQR控制律设计流程
graph TD
A[定义系统模型 A, B] --> B[设定代价函数矩阵 Q, R]
B --> C[求解 Riccati 方程]
C --> D[计算反馈增益矩阵 K]
D --> E[构建闭环控制系统]
E --> F[验证系统稳定性与性能]
4.3 LQR控制的扩展与应用
LQR 控制在理论和应用上都具有良好的扩展性。在实际工程中，常将 LQR 与非线性系统、时变系统等结合，形成扩展 LQR（Extended LQR）、离散时间 LQR（DLQR）等变种。
4.3.1 离散时间LQR的实现方法
离散时间 LQR（DLQR）适用于数字控制系统。其控制律形式为：
u_k = -K_d x_k
其中 $ K_d $ 通过离散 Riccati 方程求解得到。DLQR 的实现流程如下：
- 将连续时间系统离散化；
- 设置离散代价函数；
- 求解离散 Riccati 方程；
- 设计反馈控制律。
代码示例：使用 Python 求解离散时间 LQR
from scipy.linalg import solve_discrete_are
import numpy as np
# 离散化系统矩阵
A_d = np.array([[1, 0.1],
[0, 1]])
B_d = np.array([[0.005],
[0.1]])
Q_d = np.eye(2)
R_d = np.array([[1]])
# 求解离散 Riccati 方程
P_d = solve_discrete_are(A_d, B_d, Q_d, R_d)
# 计算反馈增益矩阵 K_d
K_d = np.linalg.inv(R_d + B_d.T @ P_d @ B_d) @ B_d.T @ P_d @ A_d
print("离散 Riccati 矩阵 P_d：")
print(P_d)
print("离散反馈增益矩阵 K_d：")
print(K_d)
代码逻辑分析：
- 使用
solve_discrete_are函数求解离散时间代数 Riccati 方程； - 系统矩阵 $ A_d $ 和 $ B_d $ 为离散化后的状态转移矩阵；
- 代价函数矩阵 $ Q_d $ 和 $ R_d $ 用于控制状态和控制输入的权重；
- 反馈增益矩阵 $ K_d $ 用于构建离散控制律；
- 该控制器适用于嵌入式系统、机器人控制、飞行器导航等应用场景。
参数说明：
-
A_d：离散状态矩阵； -
B_d：离散控制矩阵； -
Q_d：离散状态代价矩阵； -
R_d：离散控制代价矩阵； -
P_d：离散 Riccati 方程的解； -
K_d：离散反馈增益矩阵。
4.3.2 LQR在飞行器控制中的典型应用
LQR 控制广泛应用于飞行器姿态控制中，例如四旋翼无人机的姿态调节、导弹的飞行轨迹优化等。以四旋翼为例，其动力学模型可近似为线性系统，使用 LQR 控制器可实现姿态角（俯仰角、滚转角、偏航角）的快速收敛与稳定。
应用场景表格：
| 应用领域 | 系统模型 | 控制目标 | LQR优势 |
|---|---|---|---|
| 四旋翼无人机 | 线性姿态动力学 | 姿态角稳定 | 实时性强、反馈控制精度高 |
| 飞行器轨迹跟踪 | 线性化运动模型 | 轨迹跟踪误差最小化 | 支持多变量控制 |
| 工业机器人 | 关节动力学模型 | 关节角度快速响应 | 控制能量优化 |
| 自动驾驶车辆 | 纵向与横向动力学 | 路径跟踪与避障 | 稳定性与鲁棒性高 |
LQR 在这些系统中的应用不仅提高了控制精度，也增强了系统对扰动的鲁棒性。通过调整代价函数矩阵 $ Q $ 和 $ R $，还可以在控制性能与能耗之间进行灵活权衡。
通过本章内容的介绍，我们深入理解了 LQR 控制的基本建模思想、动态规划解法及其在实际系统中的应用方式。LQR 控制不仅在理论分析中具有清晰的数学结构，在工程实现中也表现出良好的性能和可扩展性。
5. 动态规划的数值求解方法与性能分析
动态规划作为优化控制和决策系统的核心方法，其数值求解过程直接影响到算法的实用性与可扩展性。随着问题规模的增大，尤其是状态空间维度的增加，传统的动态规划方法面临“维度灾难”（Curse of Dimensionality）的挑战。因此，如何高效地实现动态规划算法，并在有限计算资源下获得高质量的近似解，成为研究的重点。本章将围绕动态规划的数值计算框架、近似动态规划（ADP）技术以及算法性能评估与优化策略展开深入探讨。
5.1 动态规划的数值计算框架
5.1.1 网格化与插值技术
在动态规划的实际计算中，状态空间通常需要离散化为有限的网格点，以便在这些点上近似价值函数。这种离散化方法称为 网格化 （Grid-based Discretization）。对于一维状态空间，网格划分相对简单，但随着状态维度增加，网格点数量呈指数级增长，导致计算复杂度急剧上升。
网格化示例
以一个二维状态空间 $ x = (x_1, x_2) \in [0, 1] \times [0, 1] $ 为例，若每个维度划分 $ N = 100 $ 个点，则总网格点数为 $ N^2 = 10,000 $。若维度增加到 $ d = 5 $，则网格点数将变为 $ 10^{10} $，这显然在计算上不可行。
import numpy as np
# 二维状态空间网格化
x1 = np.linspace(0, 1, 100)
x2 = np.linspace(0, 1, 100)
X1, X2 = np.meshgrid(x1, x2)
# 价值函数在网格点上的近似值
V = np.zeros_like(X1)
# 简单赋值（模拟动态规划迭代）
for i in range(len(x1)):
for j in range(len(x2)):
V[i, j] = np.sin(X1[i, j]) + np.cos(X2[i, j])
代码逻辑分析：
- np.linspace 创建等间距网格点；
- np.meshgrid 生成二维网格坐标；
- 初始化价值函数矩阵 V ；
- 通过双重循环对每个网格点赋值，模拟动态规划迭代中价值函数的更新。
参数说明：
- x1 , x2 ：状态变量的离散化范围；
- X1 , X2 ：网格化后的状态空间；
- V[i,j] ：状态 $ (x1[i], x2[j]) $ 对应的价值函数近似值。
5.1.2 高维问题的计算挑战
高维状态空间带来的主要问题在于状态数量的指数增长（维度灾难），这使得网格化方法在实际问题中难以应用。例如在机器人控制、金融投资或大规模资源调度问题中，状态维度可能达到数十甚至上百维，传统动态规划的网格化方法无法胜任。
为应对这一挑战，研究者提出了多种替代方法，包括稀疏网格（Sparse Grid）、蒙特卡洛采样（Monte Carlo Sampling）、以及函数逼近技术等。这些方法旨在减少计算复杂度，同时保持对价值函数的良好近似能力。
5.2 近似动态规划（ADP）与函数逼近
5.2.1 线性基函数与神经网络的应用
近似动态规划（Approximate Dynamic Programming, ADP） 是处理高维状态空间问题的有效方法。其核心思想是使用函数逼近器（Function Approximator）来近似价值函数或策略函数，从而避免显式地维护状态网格。
线性基函数逼近
线性基函数（Linear Basis Function）是 ADP 中最常用的函数逼近方式之一。假设价值函数 $ V(x) $ 可以表示为一组基函数 $ \phi_i(x) $ 的线性组合：
V(x) \approx \sum_{i=1}^n w_i \phi_i(x)
其中 $ w_i $ 为权重系数，$ \phi_i(x) $ 为基函数。
例如，可以选取多项式基函数：
\phi_1(x) = 1,\quad \phi_2(x) = x,\quad \phi_3(x) = x^2
示例：线性回归逼近价值函数
from sklearn.linear_model import LinearRegression
# 生成模拟数据
X = np.random.rand(100, 1) * 10
y = np.sin(X) + np.random.normal(0, 0.1, X.shape)
# 使用线性回归拟合价值函数
model = LinearRegression()
model.fit(X, y)
# 预测
X_test = np.linspace(0, 10, 100).reshape(-1, 1)
y_pred = model.predict(X_test)
代码逻辑分析：
- 模拟状态空间数据 X 和价值函数 y ；
- 使用线性回归模型 LinearRegression 来逼近价值函数；
- 对测试集 X_test 进行预测，获得近似价值函数曲线。
参数说明：
- X : 状态变量输入；
- y : 真实价值函数输出；
- model : 训练好的线性逼近模型；
- y_pred : 近似价值函数预测值。
神经网络逼近
当状态空间高度非线性或复杂时，使用神经网络作为函数逼近器更为有效。神经网络具有强大的非线性拟合能力，可以捕捉复杂的函数关系。
import tensorflow as tf
from tensorflow.keras import layers, models
# 构建神经网络模型
model = models.Sequential([
layers.Dense(64, activation='relu', input_shape=(1,)),
layers.Dense(64, activation='relu'),
layers.Dense(1)
])
# 编译模型
model.compile(optimizer='adam', loss='mse')
# 训练模型
model.fit(X, y, epochs=200, verbose=0)
# 预测
y_pred_nn = model.predict(X_test)
代码逻辑分析：
- 构建包含两个隐藏层的神经网络；
- 使用均方误差（MSE）作为损失函数；
- 使用 Adam 优化器进行训练；
- 预测并输出神经网络逼近的价值函数。
参数说明：
- layers.Dense : 全连接层；
- activation='relu' : 激活函数；
- optimizer='adam' : 优化器；
- loss='mse' : 损失函数。
5.2.2 基于采样的近似方法
在实际系统中，状态空间可能非常大或连续，无法进行完整采样。此时， 基于采样的近似方法 （Sampling-based ADP）提供了一种有效的解决方案。
重要采样（Importance Sampling）
重要采样是一种从经验数据中估计价值函数的方法。其核心思想是根据策略的分布对状态进行采样，并计算加权平均。
示例：基于采样的价值函数估计
def importance_sampling(trajectories, target_policy, behavior_policy):
weights = []
returns = []
for traj in trajectories:
states, actions, rewards = traj
weight = 1.0
G = 0
for t in range(len(states)):
pi = target_policy(states[t], actions[t])
mu = behavior_policy(states[t], actions[t])
weight *= pi / mu
G += rewards[t]
weights.append(weight)
returns.append(G)
return np.average(returns, weights=weights)
代码逻辑分析：
- 输入 trajectories 表示采样轨迹；
- 对每条轨迹计算策略比率（target_policy / behavior_policy）；
- 计算加权回报；
- 返回加权平均值作为估计值。
参数说明：
- trajectories : 采样轨迹列表；
- target_policy : 目标策略；
- behavior_policy : 行为策略；
- G : 累计回报；
- weight : 策略比率乘积。
5.3 算法性能评估与优化策略
5.3.1 收敛速度与计算效率
动态规划算法的性能主要体现在两个方面： 收敛速度 （Convergence Rate）和 计算效率 （Computational Efficiency）。收敛速度决定了算法在给定误差范围内达到稳定解所需的迭代次数；计算效率则反映了算法在单位时间内完成的任务量。
收敛性分析
以价值迭代（Value Iteration）为例，其收敛性依赖于 Bellman 算子的压缩性质。假设折扣因子 $ \gamma < 1 $，则 Bellman 算子是压缩映射，保证价值函数最终收敛到最优解。
def value_iteration(P, R, gamma=0.9, theta=1e-8):
V = np.zeros_like(R)
while True:
delta = 0
for s in range(len(V)):
v = V[s]
V[s] = max([sum([P[s, a, s1]*(R[s, a, s1] + gamma*V[s1]) for s1 in range(len(V))]) for a in range(len(P[s]))])
delta = max(delta, abs(v - V[s]))
if delta < theta:
break
return V
代码逻辑分析：
- 初始化价值函数 V 为零；
- 迭代更新每个状态的价值；
- 若最大误差 delta 小于阈值 theta ，则停止迭代；
- 返回收敛后的价值函数。
参数说明：
- P[s, a, s1] : 状态转移概率；
- R[s, a, s1] : 奖励函数；
- gamma : 折扣因子；
- theta : 收敛阈值。
计算效率分析
计算效率通常通过算法运行时间、内存占用、迭代次数等指标衡量。高维状态空间或复杂策略会导致计算效率下降。为提升效率，可采用以下策略：
- 并行计算 ：将状态空间划分，利用多核 CPU 或 GPU 并行计算；
- 增量更新 ：仅更新变化部分，而非重新计算整个价值函数；
- 启发式剪枝 ：忽略低概率状态或动作，减少计算量。
5.3.2 并行计算与GPU加速方案
随着 GPU 在科学计算领域的广泛应用，动态规划算法的并行化成为提升性能的重要手段。利用 GPU 的大规模并行架构，可以显著加速高维状态空间的处理。
使用 GPU 加速动态规划（CUDA 示例）
__global__ void valueIterationKernel(float *V, float *V_new, float *P, float *R, float gamma, int n_states, int n_actions) {
int s = threadIdx.x + blockIdx.x * blockDim.x;
if (s >= n_states) return;
float max_val = -INFINITY;
for (int a = 0; a < n_actions; a++) {
float sum_val = 0;
for (int s1 = 0; s1 < n_states; s1++) {
sum_val += P[s * n_actions * n_states + a * n_states + s1] *
(R[s * n_actions * n_states + a * n_states + s1] + gamma * V[s1]);
}
max_val = fmaxf(max_val, sum_val);
}
V_new[s] = max_val;
}
代码逻辑分析：
- 定义 CUDA 内核函数 valueIterationKernel ；
- 每个线程处理一个状态 s ；
- 对每个动作 a ，计算期望回报；
- 更新价值函数 V_new[s] ；
- 使用 fmaxf 获取最大值。
参数说明：
- V , V_new : 当前与更新后的价值函数；
- P , R : 状态转移与奖励矩阵；
- gamma : 折扣因子；
- n_states , n_actions : 状态与动作空间维度。
性能对比表格
| 方法 | 状态维度 | 迭代次数 | GPU加速比（vs CPU） |
|---|---|---|---|
| 串行 CPU | 1000 | 100 | 1x |
| 多线程 CPU | 1000 | 100 | 4x |
| GPU加速 | 1000 | 100 | 20x |
| 并行 GPU + ADP | 10000 | 150 | 50x |
流程图：动态规划数值求解流程
graph TD
A[状态空间网格化] --> B[初始化价值函数]
B --> C[策略评估与更新]
C --> D{是否收敛?}
D -- 是 --> E[输出最优策略]
D -- 否 --> F[函数逼近更新]
F --> G[并行计算加速]
G --> C
通过本章的深入分析，我们系统地探讨了动态规划的数值求解方法，包括网格化与插值、高维状态处理、函数逼近技术、收敛性分析与并行加速方案。这些内容不仅为动态规划的实际应用提供了理论支撑，也为后续章节中处理非确定性问题与实际系统建模打下了坚实基础。
6. 非确定性动态规划与马尔科夫决策过程（MDP）
在现实世界的许多决策问题中，系统状态的转移并非总是确定性的。例如，在自动驾驶、金融投资和机器人路径规划等应用中，系统的下一状态往往受到多种随机因素的影响。为了处理这类问题，非确定性动态规划（Stochastic Dynamic Programming）应运而生，其中马尔科夫决策过程（Markov Decision Process, MDP）是最核心的建模工具之一。本章将系统讲解MDP的基本结构、Bellman期望方程的建立、以及其求解方法，并通过实际案例展示其在复杂系统中的应用。
6.1 不确定环境下的动态规划建模
6.1.1 概率状态转移与期望代价函数
在传统的确定性动态规划中，系统状态在给定控制输入后，会以唯一确定的方式转移到下一状态。然而，在非确定性环境中，状态转移具有概率性。例如，自动驾驶车辆在某个控制动作（如转向或加速）后，由于传感器误差、路面状况等因素，其下一时刻的状态可能是多个可能状态中的某一个，每个状态转移具有一定的概率。
为了刻画这种不确定性，引入 状态转移概率函数 $ P(s’ | s, a) $，表示在状态 $ s $ 下采取动作 $ a $ 后，转移到状态 $ s’ $ 的概率。此外，系统的代价函数也从单一确定值扩展为期望值形式：
J(s) = \mathbb{E} \left[ \sum_{t=0}^{\infty} \gamma^t c(s_t, a_t) \right]
其中：
- $ c(s_t, a_t) $：在状态 $ s_t $ 下采取动作 $ a_t $ 的即时代价；
- $ \gamma \in [0, 1] $：折扣因子，用于衡量未来代价的当前价值。
6.1.2 MDP的基本结构与Bellman期望方程
马尔科夫决策过程（MDP）是一个五元组 $ (S, A, P, c, \gamma) $，其中：
| 符号 | 含义 |
|---|---|
| $ S $ | 状态空间，表示所有可能的状态集合 |
| $ A $ | 动作空间，表示所有可能的控制输入集合 |
| $ P $ | 状态转移函数，即 $ P(s’ |
| $ c $ | 即时代价函数，$ c: S \times A \rightarrow \mathbb{R} $ |
| $ \gamma $ | 折扣因子，用于衡量未来代价的重要性 |
MDP 的目标是找到一个策略 $ \pi: S \rightarrow A $，使得从任意状态 $ s $ 出发的期望总代价最小化。
Bellman期望方程
对于给定策略 $ \pi $，其价值函数 $ V^\pi(s) $ 表示从状态 $ s $ 出发、按照策略 $ \pi $ 执行的期望代价：
V^\pi(s) = c(s, \pi(s)) + \gamma \sum_{s’} P(s’ | s, \pi(s)) V^\pi(s’)
这是MDP中最基础的Bellman期望方程。它反映了当前状态的代价由当前动作的即时代价加上下一状态的期望代价构成。
6.2 MDP的求解方法
6.2.1 价值迭代与策略迭代算法
在MDP中，求解最优策略的核心在于找到满足Bellman最优方程的价值函数：
V^ (s) = \min_{a \in A} \left[ c(s, a) + \gamma \sum_{s’} P(s’ | s, a) V^ (s’) \right]
该方程没有解析解，通常采用 价值迭代（Value Iteration） 和 策略迭代（Policy Iteration） 两种算法进行求解。
价值迭代算法（Value Iteration）
def value_iteration(S, A, P, c, gamma, theta=1e-6):
V = {s: 0 for s in S}
while True:
delta = 0
for s in S:
v = V[s]
# 更新V(s)
V[s] = min([c(s, a) + gamma * sum(P(s, a, s_prime) * V[s_prime] for s_prime in S) for a in A])
delta = max(delta, abs(v - V[s]))
if delta < theta:
break
# 提取最优策略
policy = {}
for s in S:
policy[s] = min(A, key=lambda a: c(s, a) + gamma * sum(P(s, a, s_prime) * V[s_prime] for s_prime in S))
return policy, V
代码逻辑分析：
-
S、A：状态空间与动作空间。 -
P(s, a, s_prime)：状态转移概率。 -
c(s, a)：即时代价函数。 -
gamma：折扣因子。 -
theta：收敛阈值，控制迭代终止条件。
在每次迭代中，算法遍历所有状态，计算每个动作下的期望代价，并更新当前状态的最优价值。当价值函数的变化小于设定的阈值时，迭代终止，并根据最终价值函数提取最优策略。
策略迭代算法（Policy Iteration）
策略迭代算法分为两个步骤： 策略评估（Policy Evaluation） 和 策略改进（Policy Improvement） 。
def policy_evaluation(policy, S, P, c, gamma, theta=1e-6):
V = {s: 0 for s in S}
while True:
delta = 0
for s in S:
v = V[s]
a = policy[s]
V[s] = c(s, a) + gamma * sum(P(s, a, s_prime) * V[s_prime] for s_prime in S)
delta = max(delta, abs(v - V[s]))
if delta < theta:
break
return V
def policy_improvement(V, S, A, P, c, gamma):
policy = {}
for s in S:
policy[s] = min(A, key=lambda a: c(s, a) + gamma * sum(P(s, a, s_prime) * V[s_prime] for s_prime in S))
return policy
def policy_iteration(S, A, P, c, gamma):
policy = {s: A[0] for s in S} # 初始化策略
while True:
V = policy_evaluation(policy, S, P, c, gamma)
new_policy = policy_improvement(V, S, A, P, c, gamma)
if new_policy == policy:
break
policy = new_policy
return policy, V
代码逻辑分析：
-
policy_evaluation：对当前策略进行价值函数评估； -
policy_improvement：基于当前价值函数改进策略； -
policy_iteration：不断迭代评估与改进，直到策略稳定。
策略迭代通常比价值迭代更快收敛，但每次评估可能需要较多的计算资源。
6.2.2 蒙特卡洛方法与时间差分学习
当状态空间和动作空间非常大甚至连续时，传统的价值迭代和策略迭代难以处理。此时可以采用 基于采样的方法 ，如：
- 蒙特卡洛方法（Monte Carlo, MC）
- 时间差分学习（Temporal Difference, TD）
蒙特卡洛策略评估（Monte Carlo Policy Evaluation）
MC方法通过完整的 episode（即从初始状态到终止状态的一次完整运行）来估计价值函数：
def monte_carlo_evaluation(policy, env, gamma=0.9, num_episodes=1000):
V = defaultdict(float)
returns = defaultdict(list)
for _ in range(num_episodes):
episode = generate_episode(policy, env)
states, actions, rewards = zip(*episode)
G = 0
for t in reversed(range(len(episode))):
G = gamma * G + rewards[t]
returns[states[t]].append(G)
V[states[t]] = np.mean(returns[states[t]])
return V
代码逻辑分析：
-
generate_episode：生成一个episode，记录状态、动作和即时奖励； -
G：累计折扣回报； - 每个状态的价值是其所有episode中G的平均值。
时间差分学习（TD Learning）
TD方法在每个时间步更新价值函数，无需等到episode结束：
def td_learning(policy, env, gamma=0.9, alpha=0.1, num_episodes=1000):
V = defaultdict(float)
for _ in range(num_episodes):
state = env.reset()
done = False
while not done:
action = policy[state]
next_state, reward, done, _ = env.step(action)
target = reward + gamma * V[next_state] * (not done)
V[state] += alpha * (target - V[state])
state = next_state
return V
代码逻辑分析：
-
alpha：学习率，控制更新步长； -
target：TD目标，即当前奖励加上下一状态的估计价值； - 每个状态的价值函数通过当前TD误差进行更新。
6.3 非确定性动态规划的应用实例
6.3.1 自动驾驶路径规划中的MDP建模
在自动驾驶路径规划中，车辆需要在复杂环境中选择最优路径，同时考虑各种不确定因素，如行人突然横穿、其他车辆变道等。
MDP建模要素：
| 要素 | 内容 |
|---|---|
| 状态 $ s $ | 车辆当前位置、速度、周围交通情况 |
| 动作 $ a $ | 加速、减速、左转、右转 |
| 转移函数 $ P(s’ | s, a) $ |
| 即时代价 $ c(s, a) $ | 与安全距离、偏离目标路径、能量消耗等相关的加权代价 |
| 折扣因子 $ \gamma $ | 通常设为0.95，以平衡短期与长期目标 |
示例流程图（mermaid）：
graph TD
A[状态s] --> B{执行动作a}
B --> C[环境反馈]
C --> D[得到下一状态s'和即时代价c]
D --> E[更新价值函数]
E --> F{是否到达目标状态？}
F -->|是| G[终止]
F -->|否| A
通过MDP建模，自动驾驶系统可以在不确定性环境中做出最优路径决策，提升安全性与效率。
6.3.2 金融投资决策中的风险控制策略
在金融投资中，资产价格的波动具有高度不确定性，因此采用MDP建模可以帮助制定风险控制策略。
MDP建模要素：
| 要素 | 内容 |
|---|---|
| 状态 $ s $ | 当前资产组合、市场趋势、风险指标 |
| 动作 $ a $ | 买入、卖出、持有 |
| 转移函数 $ P(s’ | s, a) $ |
| 即时代价 $ c(s, a) $ | 收益损失、交易成本、风险惩罚项 |
| 折扣因子 $ \gamma $ | 通常设为0.9，更关注长期收益 |
示例表格：不同动作下的期望收益与风险
| 动作 | 平均收益（%） | 波动率（%） | 风险惩罚项 | 期望净收益 |
|---|---|---|---|---|
| 买入 | 5.2 | 8.1 | 0.5 | 4.7 |
| 卖出 | 0.8 | 2.3 | 0.1 | 0.7 |
| 持有 | 3.1 | 4.9 | 0.3 | 2.8 |
基于MDP模型，投资策略可以在收益与风险之间找到最优平衡点，适用于高频交易、组合优化等多种金融场景。
本章系统介绍了非确定性动态规划的核心建模工具——马尔科夫决策过程（MDP），并深入分析了其Bellman方程、价值迭代、策略迭代等求解方法，以及蒙特卡洛和时间差分等基于采样的学习方法。最后通过自动驾驶路径规划与金融投资的实际应用案例，展示了MDP在复杂不确定性环境中的强大建模与决策能力。这些内容为后续在强化学习、智能决策系统等领域的深入研究与应用提供了坚实的理论基础。
7. 动态规划在实际系统中的应用实践
7.1 资源调度中的动态规划应用
7.1.1 多阶段调度模型的构建
在资源调度问题中，动态规划的核心在于将问题建模为一个多阶段决策过程。每个阶段代表一个决策点，状态变量表示当前资源的分配情况，而决策变量则代表当前阶段所采取的资源调度动作。
以一个简化版的电力调度问题为例，假设我们需要在T个时间周期内调度电力资源，以最小化总成本。目标函数可表示为：
$$ J = \sum_{t=0}^{T-1} C_t(x_t, u_t) + C_T(x_T) $$
其中：
- $ x_t $：表示第 $ t $ 个时间周期的状态（如剩余电量）；
- $ u_t $：表示第 $ t $ 个时间周期的控制输入（如用电量）；
- $ C_t(x_t, u_t) $：表示第 $ t $ 个周期的即时成本；
- $ C_T(x_T) $：终端成本（如电力不足的惩罚）；
状态转移函数可以表示为：
$$ x_{t+1} = f(x_t, u_t) $$
例如，在电力系统中，状态转移可能表示为：
$$ x_{t+1} = x_t - u_t + \text{发电量}_t $$
动态规划通过递归地从终端状态向前计算，求解出每一步的最优策略。我们可以使用价值迭代方法进行求解，具体步骤如下：
def value_iteration(states, actions, transition, cost, T):
V = {s: 0 for s in states} # 初始化价值函数
for t in reversed(range(T)):
V_new = {}
for s in states:
min_cost = float('inf')
for a in actions:
next_s = transition(s, a)
c = cost(s, a)
total_cost = c + V[next_s]
if total_cost < min_cost:
min_cost = total_cost
V_new[s] = min_cost
V = V_new
return V
说明：
-states：所有可能的状态集合；
-actions：所有可能的动作集合；
-transition：状态转移函数；
-cost：代价函数；
-T：总时间步数；
- 每次迭代计算每个状态的最优价值函数值；
7.1.2 动态规划在能源管理中的实践
在能源管理系统中，动态规划被广泛用于优化发电与储能策略。例如，风力发电的波动性使得电网调度变得复杂。通过动态规划建模，可以在不同时间步中动态调整储能系统的充放电策略，以最小化能源浪费和调度成本。
一个典型的优化目标函数如下：
$$ \min \sum_{t=0}^{T} \left[ \alpha (u_t^{\text{charge}})^2 + \beta (u_t^{\text{discharge}})^2 + \gamma (E_t - E_{\text{target}})^2 \right] $$
其中：
- $ u_t^{\text{charge}} $：第 $ t $ 个周期的充电功率；
- $ u_t^{\text{discharge}} $：第 $ t $ 个周期的放电功率；
- $ E_t $：当前储能容量；
- $ \alpha, \beta, \gamma $：权重系数；
通过建立状态转移函数与目标函数，结合动态规划算法，可以在每个时间点动态优化储能系统的控制策略。
7.2 机器人路径规划中的动态规划实现
7.2.1 状态空间建模与目标函数设计
机器人路径规划是动态规划的典型应用场景之一。通常，机器人在一个二维或三维空间中移动，其状态包括位置、速度、方向等。目标是找到一条从起点到终点的最优路径，使得路径的代价最小。
假设机器人的状态为 $ x_t = (x, y, \theta) $，即位置与朝向，控制输入为 $ u_t = (\nu, \omega) $，即线速度和角速度。目标函数可定义为：
$$ J = \sum_{t=0}^{T} | x_t - x_{\text{goal}} | + \lambda | u_t | $$
其中：
- $ | x_t - x_{\text{goal}} | $：当前位置与目标的距离；
- $ | u_t | $：控制输入的大小，用于平滑路径；
- $ \lambda $：控制输入的权重系数；
状态转移函数为：
$$ x_{t+1} = x_t + \nu \cos(\theta) \Delta t $$
$$ y_{t+1} = y_t + \nu \sin(\theta) \Delta t $$
$$ \theta_{t+1} = \theta_t + \omega \Delta t $$
使用动态规划求解该问题时，可以采用价值函数近似或策略迭代方法，结合网格化或采样方法进行状态空间的离散化处理。
7.2.2 实时路径优化与避障策略
在动态环境中，机器人需要根据传感器反馈实时调整路径。此时可以采用在线动态规划方法（Online DP）或模型预测控制（MPC）结合动态规划思想进行实时路径优化。
例如，使用MPC框架，每一步预测未来N步的状态与控制输入，求解有限时间范围内的最优控制序列：
$$ \min_{u_{t|t}, \dots, u_{t+N|t}} \sum_{k=0}^{N} \left( |x_{t+k|t} - x_{\text{goal}}|^2 + \lambda |u_{t+k|t}|^2 \right) $$
约束条件包括：
- 动力学模型：$ x_{t+k+1|t} = f(x_{t+k|t}, u_{t+k|t}) $
- 障碍物避开：$ g(x_{t+k|t}) \geq 0 $
通过滚动优化策略，机器人可以在动态环境中实时调整路径，实现避障与目标趋近的双重目标。
7.3 通信网络优化中的动态规划方法
7.3.1 网络流量调度与资源分配问题
在通信网络中，动态规划被用于优化流量调度和资源分配。假设我们有一个具有多个节点的网络，每个节点在每个时间周期内可以发送、接收或转发数据包。
目标是设计一种调度策略，使得整个网络的吞吐量最大化，同时最小化延迟和丢包率。
建模为动态规划问题时，状态变量可以是每个节点的当前缓冲区大小，动作变量为路由决策（选择下一跳节点），目标函数可定义为：
$$ J = \sum_{t=0}^{T} \left( -\alpha \cdot \text{delay}_t + \beta \cdot \text{throughput}_t \right) $$
其中：
- $ \text{delay}_t $：第 $ t $ 个周期的平均数据包延迟；
- $ \text{throughput}_t $：第 $ t $ 个周期的总吞吐量；
- $ \alpha, \beta $：权衡延迟与吞吐量的参数；
状态转移函数则由网络的拓扑结构和数据包的路由策略决定。
7.3.2 动态路由优化与QoS保障机制
服务质量（QoS）保障是通信网络优化中的关键问题。通过动态规划方法，可以在多阶段决策中动态调整路由策略，以满足带宽、时延、抖动等指标。
例如，可以使用强化学习结合动态规划的思想（如Q-learning）来实现动态路由优化：
def q_learning_update(Q, state, action, reward, next_state, alpha=0.1, gamma=0.99):
best_next_action = max(Q[next_state], key=Q[next_state].get)
td_target = reward + gamma * Q[next_state][best_next_action]
td_error = td_target - Q[state][action]
Q[state][action] += alpha * td_error
return Q
说明：
-Q：Q值表，表示每个状态-动作对的预期回报；
-state：当前状态；
-action：当前动作；
-reward：执行动作后的即时奖励；
-next_state：执行动作后的新状态；
-alpha：学习率；
-gamma：折扣因子；
通过不断更新Q值表，网络可以学习到在不同状态下选择最优的路由路径，从而实现QoS保障与资源动态调度。
（章节内容未完，继续阅读下一章以了解深度强化学习与动态规划的融合应用）
简介：《动态规划与最优控制》是Bertsekas教授的经典教材，系统介绍了动态规划的基本理论与最优控制方法。作为优化理论与决策过程的重要工具，动态规划被广泛应用于工程、经济、计算机科学等多个领域。本书深入浅出地讲解了状态空间、Bellman方程、离散与连续时间动态规划、线性二次型控制（LQR）、数值求解方法、最优控制原理以及马尔科夫决策过程等内容。通过大量应用实例，如资源调度、机器人路径规划和通信网络优化，帮助读者掌握动态规划在实际问题中的建模与求解技巧，是控制理论与优化领域的重要参考书。
北京人形旗下天工造物具身智能开源社区，聚焦具身天工与慧思开物两大平台
更多推荐
17
0- 0
已为社区贡献7条内容
所有评论(0)