动态规划逆向求解实战从LQR控制到HJB方程的Python实现在控制理论领域动态规划提供了一种系统性的方法来解决最优控制问题。与传统的变分法不同动态规划通过逆向递推的方式从终点回溯求解最优控制策略这种方法特别适合处理具有多阶段决策特性的控制问题。本文将带您从离散时间的LQR控制入手逐步深入到连续时间的HJB方程求解通过Python代码实现完整的求解链路。1. 动态规划基础与离散LQR控制动态规划的核心思想是贝尔曼最优性原理一个最优策略的子策略也必须是最优的。在离散时间系统中这个原理可以转化为逆向递推的数值计算方法。考虑一个经典的线性二次调节器(LQR)问题其状态方程和代价函数分别为import numpy as np from scipy.linalg import solve_discrete_are # 系统参数 A np.array([[1.0, 0.1], [0.0, 1.0]]) # 状态转移矩阵 B np.array([[0.0], [0.1]]) # 控制输入矩阵 Q np.eye(2) # 状态代价矩阵 R np.eye(1) # 控制代价矩阵 N 50 # 时间步长LQR问题的动态规划解法可以通过以下逆向递推实现def lqr_backward_recursion(A, B, Q, R, N): P [None] * (N 1) # 存储每个时间步的代价矩阵 K [None] * N # 存储反馈增益矩阵 # 终端条件 P[N] Q # 逆向递推 for k in range(N-1, -1, -1): P_k P[k1] K[k] -np.linalg.solve(R B.T P_k B, B.T P_k A) P[k] Q A.T P_k A A.T P_k B K[k] return K, P这个实现展示了动态规划在离散时间系统中的典型应用从终端条件出发逆向计算每个时间步的最优控制策略。2. 连续系统的HJB方程与数值解法当我们将问题扩展到连续时间系统时动态规划导出了著名的Hamilton-Jacobi-Bellman(HJB)方程。HJB方程描述了最优代价函数必须满足的偏微分方程$$ -\frac{\partial J^}{\partial t} \min_u \left[ L(x,u,t) \left(\frac{\partial J^}{\partial x}\right)^T f(x,u,t) \right] $$对于一般的非线性系统HJB方程难以解析求解我们需要借助数值方法。下面介绍一种基于网格的数值解法import numpy as np from scipy.interpolate import interpn from scipy.optimize import minimize def solve_hjb(f, L, phi, x_grid, t_grid, u_space): 数值求解HJB方程 参数: f: 系统动力学方程 dx/dt f(x,u,t) L: 运行时代价函数 phi: 终端代价函数 x_grid: 状态空间网格 t_grid: 时间网格 u_space: 控制输入空间 nx len(x_grid) nt len(t_grid) J np.zeros((nx, nt)) u_opt np.zeros((nx, nt)) # 设置终端条件 J[:,-1] phi(x_grid) # 逆向时间递推 for k in range(nt-2, -1, -1): t t_grid[k] dt t_grid[k1] - t for i in range(nx): x x_grid[i] # 定义优化目标函数 def objective(u): x_next x f(x, u, t) * dt J_next interpn(x_grid, J[:,k1], x_next) return L(x, u, t)*dt J_next # 求解最优控制 res minimize(objective, x0u_space[0], bounds[(u_space[0], u_space[-1])]) u_opt[i,k] res.x J[i,k] res.fun return J, u_opt这种方法虽然计算量较大但可以处理一般的非线性系统。对于高维系统还需要考虑维度灾难问题可能需要采用更高效的数值方法。3. 符号计算与自动微分在HJB求解中的应用对于复杂的系统手动推导HJB方程中的各项可能非常繁琐。我们可以利用SymPy进行符号计算自动生成必要的数学表达式from sympy import symbols, Function, diff, Eq, solve def symbolic_hjb(): # 定义符号变量 x1, x2 symbols(x1 x2) u symbols(u) t symbols(t) J Function(J)(x1, x2, t) # 定义系统动力学 f1 -x1 x2 f2 x1 - x2 u # 定义代价函数 L x1**2 x2**2 u**2 # 计算哈密顿函数 H L diff(J, x1)*f1 diff(J, x2)*f2 # 求极小值条件 u_opt solve(diff(H, u), u)[0] # 代入得到HJB方程 HJB_eq Eq(diff(J, t) H.subs(u, u_opt), 0) return HJB_eq, u_opt这种方法特别适合理论研究阶段可以自动生成复杂的数学表达式减少手工推导的错误。4. 实际案例倒立摆系统的HJB控制让我们考虑一个经典的控制问题倒立摆的平衡控制。系统动力学可以表示为$$ \dot{x}_1 x_2 \ \dot{x}_2 \frac{mgl\sin x_1 - b x_2 u}{ml^2} $$其中$x_1$是摆角$x_2$是角速度$u$是控制力矩。我们可以将前面的数值HJB求解方法应用于这个问题def inverted_pendulum(): # 系统参数 m 0.1 # 质量(kg) l 0.5 # 长度(m) b 0.1 # 摩擦系数 g 9.8 # 重力加速度 # 定义动力学 def f(x, u, t): x1, x2 x dx1 x2 dx2 (m*g*l*np.sin(x1) - b*x2 u) / (m*l**2) return np.array([dx1, dx2]) # 定义代价函数 def L(x, u, t): return x[0]**2 0.1*x[1]**2 0.01*u**2 # 终端代价 def phi(x): return 0 # 定义网格 x1_grid np.linspace(-np.pi, np.pi, 51) x2_grid np.linspace(-5, 5, 51) t_grid np.linspace(0, 10, 101) u_space np.linspace(-5, 5, 21) # 求解HJB方程 J, u_opt solve_hjb_2d(f, L, phi, x1_grid, x2_grid, t_grid, u_space) return J, u_opt这个案例展示了如何将理论方法应用于实际控制问题。虽然倒立摆系统相对简单但同样的方法可以推广到更复杂的系统。5. 性能优化与实用技巧在实际应用中HJB方程的数值求解可能面临计算效率问题。以下是一些提高性能的技巧并行计算HJB求解中每个网格点的计算是独立的适合并行化稀疏网格在高维情况下采用稀疏网格方法缓解维度灾难自适应网格在状态变化剧烈区域使用更密的网格值函数逼近使用神经网络等函数逼近器代替网格存储from joblib import Parallel, delayed def parallel_hjb_solve(f, L, phi, x_grid, t_grid, u_space): # ...其他代码... # 并行化处理每个时间步 def process_time_step(k): # 处理单个时间步的计算 pass results Parallel(n_jobs4)(delayed(process_time_step)(k) for k in range(nt-2, -1, -1)) # 合并结果 # ...这些优化技术可以显著提高求解效率使得HJB方法能够应用于更实际的控制问题。
动态规划逆向求解指南：用Python实现LQR控制与HJB方程数值解法
动态规划逆向求解实战从LQR控制到HJB方程的Python实现在控制理论领域动态规划提供了一种系统性的方法来解决最优控制问题。与传统的变分法不同动态规划通过逆向递推的方式从终点回溯求解最优控制策略这种方法特别适合处理具有多阶段决策特性的控制问题。本文将带您从离散时间的LQR控制入手逐步深入到连续时间的HJB方程求解通过Python代码实现完整的求解链路。1. 动态规划基础与离散LQR控制动态规划的核心思想是贝尔曼最优性原理一个最优策略的子策略也必须是最优的。在离散时间系统中这个原理可以转化为逆向递推的数值计算方法。考虑一个经典的线性二次调节器(LQR)问题其状态方程和代价函数分别为import numpy as np from scipy.linalg import solve_discrete_are # 系统参数 A np.array([[1.0, 0.1], [0.0, 1.0]]) # 状态转移矩阵 B np.array([[0.0], [0.1]]) # 控制输入矩阵 Q np.eye(2) # 状态代价矩阵 R np.eye(1) # 控制代价矩阵 N 50 # 时间步长LQR问题的动态规划解法可以通过以下逆向递推实现def lqr_backward_recursion(A, B, Q, R, N): P [None] * (N 1) # 存储每个时间步的代价矩阵 K [None] * N # 存储反馈增益矩阵 # 终端条件 P[N] Q # 逆向递推 for k in range(N-1, -1, -1): P_k P[k1] K[k] -np.linalg.solve(R B.T P_k B, B.T P_k A) P[k] Q A.T P_k A A.T P_k B K[k] return K, P这个实现展示了动态规划在离散时间系统中的典型应用从终端条件出发逆向计算每个时间步的最优控制策略。2. 连续系统的HJB方程与数值解法当我们将问题扩展到连续时间系统时动态规划导出了著名的Hamilton-Jacobi-Bellman(HJB)方程。HJB方程描述了最优代价函数必须满足的偏微分方程$$ -\frac{\partial J^}{\partial t} \min_u \left[ L(x,u,t) \left(\frac{\partial J^}{\partial x}\right)^T f(x,u,t) \right] $$对于一般的非线性系统HJB方程难以解析求解我们需要借助数值方法。下面介绍一种基于网格的数值解法import numpy as np from scipy.interpolate import interpn from scipy.optimize import minimize def solve_hjb(f, L, phi, x_grid, t_grid, u_space): 数值求解HJB方程 参数: f: 系统动力学方程 dx/dt f(x,u,t) L: 运行时代价函数 phi: 终端代价函数 x_grid: 状态空间网格 t_grid: 时间网格 u_space: 控制输入空间 nx len(x_grid) nt len(t_grid) J np.zeros((nx, nt)) u_opt np.zeros((nx, nt)) # 设置终端条件 J[:,-1] phi(x_grid) # 逆向时间递推 for k in range(nt-2, -1, -1): t t_grid[k] dt t_grid[k1] - t for i in range(nx): x x_grid[i] # 定义优化目标函数 def objective(u): x_next x f(x, u, t) * dt J_next interpn(x_grid, J[:,k1], x_next) return L(x, u, t)*dt J_next # 求解最优控制 res minimize(objective, x0u_space[0], bounds[(u_space[0], u_space[-1])]) u_opt[i,k] res.x J[i,k] res.fun return J, u_opt这种方法虽然计算量较大但可以处理一般的非线性系统。对于高维系统还需要考虑维度灾难问题可能需要采用更高效的数值方法。3. 符号计算与自动微分在HJB求解中的应用对于复杂的系统手动推导HJB方程中的各项可能非常繁琐。我们可以利用SymPy进行符号计算自动生成必要的数学表达式from sympy import symbols, Function, diff, Eq, solve def symbolic_hjb(): # 定义符号变量 x1, x2 symbols(x1 x2) u symbols(u) t symbols(t) J Function(J)(x1, x2, t) # 定义系统动力学 f1 -x1 x2 f2 x1 - x2 u # 定义代价函数 L x1**2 x2**2 u**2 # 计算哈密顿函数 H L diff(J, x1)*f1 diff(J, x2)*f2 # 求极小值条件 u_opt solve(diff(H, u), u)[0] # 代入得到HJB方程 HJB_eq Eq(diff(J, t) H.subs(u, u_opt), 0) return HJB_eq, u_opt这种方法特别适合理论研究阶段可以自动生成复杂的数学表达式减少手工推导的错误。4. 实际案例倒立摆系统的HJB控制让我们考虑一个经典的控制问题倒立摆的平衡控制。系统动力学可以表示为$$ \dot{x}_1 x_2 \ \dot{x}_2 \frac{mgl\sin x_1 - b x_2 u}{ml^2} $$其中$x_1$是摆角$x_2$是角速度$u$是控制力矩。我们可以将前面的数值HJB求解方法应用于这个问题def inverted_pendulum(): # 系统参数 m 0.1 # 质量(kg) l 0.5 # 长度(m) b 0.1 # 摩擦系数 g 9.8 # 重力加速度 # 定义动力学 def f(x, u, t): x1, x2 x dx1 x2 dx2 (m*g*l*np.sin(x1) - b*x2 u) / (m*l**2) return np.array([dx1, dx2]) # 定义代价函数 def L(x, u, t): return x[0]**2 0.1*x[1]**2 0.01*u**2 # 终端代价 def phi(x): return 0 # 定义网格 x1_grid np.linspace(-np.pi, np.pi, 51) x2_grid np.linspace(-5, 5, 51) t_grid np.linspace(0, 10, 101) u_space np.linspace(-5, 5, 21) # 求解HJB方程 J, u_opt solve_hjb_2d(f, L, phi, x1_grid, x2_grid, t_grid, u_space) return J, u_opt这个案例展示了如何将理论方法应用于实际控制问题。虽然倒立摆系统相对简单但同样的方法可以推广到更复杂的系统。5. 性能优化与实用技巧在实际应用中HJB方程的数值求解可能面临计算效率问题。以下是一些提高性能的技巧并行计算HJB求解中每个网格点的计算是独立的适合并行化稀疏网格在高维情况下采用稀疏网格方法缓解维度灾难自适应网格在状态变化剧烈区域使用更密的网格值函数逼近使用神经网络等函数逼近器代替网格存储from joblib import Parallel, delayed def parallel_hjb_solve(f, L, phi, x_grid, t_grid, u_space): # ...其他代码... # 并行化处理每个时间步 def process_time_step(k): # 处理单个时间步的计算 pass results Parallel(n_jobs4)(delayed(process_time_step)(k) for k in range(nt-2, -1, -1)) # 合并结果 # ...这些优化技术可以显著提高求解效率使得HJB方法能够应用于更实际的控制问题。
