用MATLAB搞定最优控制：梯度法实战教程（附完整代码）

MATLAB梯度法实战：最优控制问题的高效数值解法

引言：最优控制问题的工程挑战

在工程实践中，我们经常遇到需要动态系统在满足特定约束条件下达到最优性能的问题。这类问题在航空航天、机器人控制、工业过程优化等领域尤为常见。传统解析解法在面对复杂系统时往往束手无策，而数值解法凭借其强大的适应性和计算效率成为工程师的首选工具。

MATLAB作为科学计算领域的标准工具，提供了完善的数值计算环境和丰富的算法库，使其成为实现最优控制数值解法的理想平台。梯度法（Gradient Method）作为最优控制问题中最基础、最直观的数值解法之一，不仅原理易于理解，实现也相对简单，特别适合作为学习最优控制数值解法的切入点。

本文将聚焦于无约束最优控制问题的梯度法实现，通过一个典型实例，详细讲解如何将数学理论转化为可执行的MATLAB代码。不同于纯理论推导，我们更关注算法实现中的实际问题：如何设置迭代终止条件？如何处理数值积分误差？如何选择适当的步长加速收敛？这些都是在教科书上很少涉及，但在实际编程中至关重要的问题。

1. 问题建模与算法原理

1.1 最优控制问题的数学表述

我们考虑如下一类无约束最优控制问题：给定系统状态方程

dx/dt = f(x,u,t), x(t₀) = x₀

性能指标（目标函数）为

J(u) = Φ(x(T),T) + ∫[t₀,T] L(x,u,t) dt

目标是找到控制量u*(t)使得J(u)最小化。

实例分析：我们以《最优化与最优控制》中的经典问题为例：

• 状态方程：dx/dt = -x² + u
• 初始条件：x(0) = 10
• 性能指标：J(u) = 0.5∫₀¹ (x² + u²) dt

1.2 梯度法的核心思想

梯度法基于一个直观的优化原理：沿着目标函数下降最快的方向（负梯度方向）逐步调整控制量，直至达到最优解。具体到最优控制问题，算法流程可分为以下几个关键步骤：

1. 构造哈密顿函数：

   H = L + λᵀf

   其中λ(t)称为协态变量，满足伴随方程：

   dλ/dt = -∂H/∂x, λ(T) = ∂Φ/∂x(T)

2. 控制更新规则： 梯度方向由∂H/∂u给出，控制量更新公式为：

   u⁽ᵏ⁺¹⁾ = u⁽ᵏ⁾ - K⁽ᵏ⁾(∂H/∂u)⁽ᵏ⁾

   其中K⁽ᵏ⁾为步长参数。

3. 迭代终止条件： 当梯度范数‖∂H/∂u‖小于预设阈值ε时，算法终止。

1.3 数值实现的关键挑战

将上述理论转化为实际代码时，我们需要解决几个关键问题：

• 时间离散化：连续时间问题需要离散化处理
• 正/反向积分：状态方程正向积分，伴随方程反向积分
• 插值处理：不同方程在不同时间网格上求解，需要数据对齐
• 步长选择：固定步长与自适应步长的取舍

以下表格对比了理论方程与数值实现的关键差异：

2. MATLAB实现详解

2.1 基础环境设置

我们首先初始化算法参数和计算环境：

%% 基础参数设置 clear; clc; close all; eps = 1e-3; % 收敛阈值 t0 = 0; % 初始时间 tf = 1; % 终止时间 t_segment = 50; % 时间分段数 Tu = linspace(t0, tf, t_segment); % 控制量时间网格 % 初始猜测控制量（全零初始化） u = zeros(1, t_segment); % 初始状态和协态 initx = 10; % 初始状态x(0) initp = 0; % 终端协态λ(T)=0 max_iteration = 50; % 最大迭代次数 options = odeset('RelTol',1e-4, 'AbsTol',1e-4); % 积分器选项

2.2 核心算法实现

梯度法的主循环包含三个关键部分：正向积分状态方程、反向积分伴随方程、控制量更新。

%% 梯度法主循环 for iter = 1:max_iteration fprintf('迭代次数: %d\n', iter); % 1. 正向积分状态方程 [Tx, X] = ode45(@(t,x) stateEq(t,x,u,Tu), [t0 tf], initx, options); x = X'; % 2. 反向积分伴随方程 [Tp, P] = ode45(@(t,p) costateEq(t,p,x,Tx), [tf t0], initp, options); p = P'; p = interp1(Tp, p, Tx); % 对齐时间网格 % 3. 计算梯度 dH = gradient(p, Tx, u, Tu); grad_norm = norm(dH); % 4. 检查收敛条件 if grad_norm < eps fprintf('收敛于迭代次数: %d\n', iter); break; end % 5. 更新控制量（固定步长） K = 0.55; % 步长参数 u = u - K*dH; end

2.3 辅助函数实现

三个核心辅助函数分别对应状态方程、伴随方程和梯度计算：

% 状态方程 function dx = stateEq(t, x, u, Tu) u_interp = interp1(Tu, u, t); % 控制量插值 dx = -x^2 + u_interp; end % 伴随方程 function dp = costateEq(t, p, x, Tx) x_interp = interp1(Tx, x, t); % 状态量插值 dp = -x_interp + 2*p.*x_interp; end % 梯度计算 function dH = gradient(p, Tx, u, Tu) u_interp = interp1(Tu, u, Tx); % 控制量插值 dH = u_interp + p; % ∂H/∂u = u + λ end

关键细节说明：

1. interp1函数用于在不同时间网格间插值，确保变量对齐
2. 伴随方程从tf到t0反向积分，需注意时间方向
3. 步长K的选择直接影响收敛速度，需要权衡稳定性和效率

3. 算法优化与进阶技巧

3.1 步长选择策略

固定步长虽然简单，但可能影响收敛效率。我们可以引入Armijo准则实现自适应步长选择：

function step_size = armijo(u, dH, J, x, Tx, Tu) beta = 0.55; % 收缩系数 sigma = 0.5; % 可接受条件参数 m = 0; max_m = 20; while m <= max_m u_new = u - beta^m * dH; J_new = compute_J(u_new, x, Tx, Tu); if J_new <= J + sigma*beta^m*dH'*dH break; end m = m + 1; end step_size = m; end function J = compute_J(u, x, Tx, Tu) % 计算性能指标J的具体实现 % ... end

3.2 共轭梯度法改进

共轭梯度法通过引入历史梯度信息，可以加速收敛：

% 在梯度法基础上修改控制更新部分 if iter == 1 d = -dH; else beta = (dH'*dH)/(dH_prev'*dH_prev); d = -dH + beta*d_prev; end % 更新控制量 u = u + K*d; dH_prev = dH; d_prev = d;

3.3 结果可视化与分析

完整的可视化代码可以帮助我们直观理解算法性能：

%% 结果可视化 figure('Position', [100, 100, 900, 700]); % 状态和控制量轨迹 subplot(2,2,1); plot(Tx, x, 'LineWidth', 2); hold on; stairs(Tu, u, 'r--', 'LineWidth', 1.5); xlabel('时间'); ylabel('值'); legend('状态x(t)', '控制u(t)'); title('状态与控制量'); % 性能指标迭代过程 subplot(2,2,2); semilogy(1:iter, J_history(1:iter), 'o-', 'LineWidth', 2); xlabel('迭代次数'); ylabel('J'); title('性能指标收敛过程'); % 梯度范数变化 subplot(2,2,3); semilogy(1:iter, grad_norm_history(1:iter), 's-', 'LineWidth', 2); xlabel('迭代次数'); ylabel('‖∇J‖'); title('梯度范数变化'); % 哈密顿函数沿时间分布 subplot(2,2,4); plot(Tx, H, 'LineWidth', 2); xlabel('时间'); ylabel('H'); title('哈密顿函数');

4. 工程实践中的经验分享

在实际项目应用中，我们发现以下几个经验特别值得分享：

1. 时间网格的选择：

   • 对于平滑问题，50-100个时间点通常足够
   • 对于剧烈变化的控制，可能需要局部加密网格
   • 可以尝试非均匀网格，在变化剧烈区域增加节点

2. 初始猜测的影响：

   • 好的初始猜测能显著减少迭代次数
   • 对于周期性问题，可以考虑用前一周期的解作为初始猜测
   • 物理直觉常常能提供有价值的初始控制形状

3. 调试技巧：

   • 先验证单个迭代步的正确性
   • 检查哈密顿函数是否单调递减
   • 可视化中间结果有助于发现问题

4. 性能瓶颈分析：

   • 使用MATLAB Profiler识别耗时环节
   • 数值积分通常是计算瓶颈，可适当降低精度要求
   • 考虑将关键部分用MEX文件实现

以下是一个典型工程问题的参数选择参考表：

通过本教程的详细讲解和完整代码实现，读者应该能够掌握梯度法求解最优控制问题的核心要点，并能够根据具体工程问题进行调整和优化。数值解法虽然不像解析解那样完美，但其强大的适用性使其成为解决复杂工程问题的有力工具。