前两天终于把基于ε-纳什均衡的航天器末端追逃博弈仿真完整跑通了。整个过程比我预想的曲折：理论上看，微分对策和纳什均衡的框架很清晰，但一落到Matlab里，打靶法初值敏感、哈密顿函数数值奇异、ε值的收敛判据这些问题轮番上场。今天把这套“理论+代码+调参经验”一起整理出来，希望能给正在啃追逃博弈仿真的朋友一点参考。

这套仿真的核心是：在航天器末端接近过程中，追逃双方都具备轨道机动能力，各自按照自己的目标函数选择控制策略，最终达到一个双方都不愿意单方面改变策略的“均衡状态”。由于严格纳什均衡在连续时间动态博弈里往往难以数值实现，所以引入了ε-纳什均衡——允许策略偏离带来的收益损失不超过一个很小的ε，从而让数值求解变得可行，同时又保留了博弈论意义上的稳定性。

适合谁看？如果你是刚接触微分对策、想用Matlab实现追逃博弈仿真，或者已经实现了单纯的最优控制但觉得“不够对抗性”，这篇文章能帮你把理论模型、代码结构和调试思路串起来。

1. 为什么末端追逃要引入ε-纳什均衡：先把博弈模型说清楚

1.1 一个“追得上”和“逃得掉”的博弈

末端追逃问题本质上是一个双人零和动态博弈。追击者的目标是让末端相对距离尽可能小，最好在限定时间内实现“足够接近”；逃逸者的目标恰好相反，尽可能拉开距离，避免被“捕获”。这里的“追逃”不是单纯的一方做最优控制，另一方被动接受，而是双方都在线决策，互相影响。

举一个生活化的类比：两辆车在路口博弈，一辆想并线超车，另一辆想守住车道。超车车辆踩油门的同时，守道车辆也在调整位置；任何一方的策略效果都依赖于对方的即时反应。航天器末端追逃也是这个道理，只不过运动方程换成了轨道力学，控制输入换成了推力加速度，目标函数换成了终端相对距离或其他收益指标。

这个博弈的数学表达式可以写成：

• 追击者：最小化终端距离 J_p = ||r(t_f)||，控制约束 ||u_p|| ≤ u_{p,max}
• 逃逸者：最大化终端距离 J_e = ||r(t_f)||，控制约束 ||u_e|| ≤ u_{e,max}

注意，追击者想最小化，逃逸者想最大化，所以 J_p = -J_e，是一个典型的零和博弈。如果把终端距离写作对双方统一的代价函数，那么追击者是极小化者，逃逸者是极大化者。

1.2 从严格纳什到ε-纳什：数值求解的现实妥协

在有限时间、连续状态、连续控制的微分对策中，严格纳什均衡往往对应一组哈密顿-雅可比-艾萨克斯方程（HJIE）的解。理论很漂亮，但实际求解析解极其困难，绝大多数情况下只能数值求解。

即使采用间接法，把问题转化为两点边值问题（TPBVP），我们面对的也是一个强非线性、初值敏感的方程组。用打靶法迭代时，稍微偏离正确的协态初值，终端状态就会“飞出银河系”。严格纳什均衡要求双方策略对各自目标函数实现精确最优，这在数值计算中几乎是不可达的。

ε-纳什均衡适时地松动了这个条件。它的定义很直接：

如果一方单方面改变策略，其收益最多只能额外增加ε，那么当前策略组合就是一个ε-纳什均衡。

换句话说，双方当前策略已经足够好，即使对方或自己“反悔”，所能获得的收益提升也不会超过一个很小的阈值ε。这在工程上是可以接受的，因为任何数值求解本身就有离散误差、迭代误差，强行追求精确均衡反而会让算法不收敛。

1.3 这个约定对仿真结果意味着什么

引入ε-纳什均衡之后，我们可以把仿真输出解释为“近似最优策略对”。在仿真复盘时，你可以做单方面策略扰动实验：固定追逸者的策略，只改变追击者的控制序列，看看终端距离能比当前结果提升多少；如果提升量小于ε，说明当前追击者策略是ε-最优的。对逃逸者也做同样的检验。

这一步在理论分析中特别重要。因为数值解本身不告诉你“这是均衡”，只有通过扰动检验，才能确认这个策略组合确实具备博弈意义上的稳定性。

2. 末端追逃的建模：相对轨道动力学与博弈描述

2.1 坐标系与相对运动方程

末端追逃仿真最常用的是目标轨道坐标系下的相对运动方程。当追击者与逃逸者之间的相对距离远小于轨道半径时，可以用线性化的CW方程（Clohessy-Wiltshire）描述。

设状态向量为：

x = [rx, ry, rz, vx, vy, vz]^T

其中 r = [rx, ry, rz] 是逃逸者相对追击者的位置向量，v = [vx, vy, vz] 是相对速度。轨道角速度 n 由参考轨道的半长轴决定，n = sqrt(μ / a^3)。

CW方程的标准形式为：

• 径向: x'' = 3n²x + 2n y' + u_x
• 迹向: y'' = -2n x' + u_y
• 法向: z'' = -n² z + u_z

这里的 u = u_p - u_e，即追击者控制加速度与逃逸者控制加速度的差。注意，如果分别对追击者和逃逸者施加控制，追逃双方的绝对轨道机动会同时作用于相对运动。

在Matlab里，我把这个方程封装成一个函数，输入当前状态和控制差，输出状态导数。这一步是所有后续求解的地基。

2.2 博弈的代价函数和终端条件

末端追逃的典型终端条件是：

• 固定终端时间 t_f
• 终端目标：追击者试图让 ||r(t_f)|| ≤ r_capture，逃逸者试图让 ||r(t_f)|| > r_capture

这里的 r_capture 可以理解为“捕获半径”，在实际仿真中依据任务场景设定，比如几十米或者几百米。终端距离小于捕获半径，就认为追击成功；大于捕获半径且时间耗尽，则逃逸成功。

但纯终端距离指标有个问题：如果追击者只需在最后时刻“冲”到很近的地方，那么轨迹可能非常激进，不符合实际轨道机动约束。所以通常会在代价函数中加入控制能量项：

J = ||r(t_f)||² + ∫[ (u_p^T R_p u_p) - (u_e^T R_e u_e) ] dt

追击者希望 J 小，逃逸者希望 J 大。这里的积分项可以理解为控制代价，避免策略无限消耗燃料。

2.3 末端博弈的“有限时间”特性

末日追逃是有限时域博弈，终端时间 t_f 固定。这个设定与现实任务一致：追逃窗口由轨道预报决定，超过窗口后的事件不再纳入博弈模型。

有限时间特性给数值求解带来两个影响：

1. 协态变量的终端条件可以直接由终端代价函数的梯度给出，不需要横截条件的复杂推导。
2. 哈密顿函数在终端附近可能出现奇异性，尤其是控制约束起作用时（bang-bang控制），需要做光滑化处理。

我把这些特征全部纳入后续的求解框架。模型越贴近实际，数值求解的“脾气”越大，但只要每一步都有明确物理含义，调试起来还是能找到方向的。

3. Matlab仿真架构设计与核心代码实现

3.1 模块划分：从初始化到后处理

写仿真代码最忌讳把整个博弈过程塞进一个几千行的脚本里。我的做法是把工程拆成四个模块：

1. 初始化模块：设置轨道参数、控制约束、终端时间、玩家策略初值。
2. 博弈求解模块：核心策略迭代与TPBVP求解。
3. 仿真推进模块：用离散化的相对运动方程做开环或闭环轨迹仿真。
4. 后处理模块：计算终端距离、绘制双方轨迹、做单方面策略扰动检验。

这种划分方式让我在调试时可以单独测试每个模块。比如先不跑策略迭代，直接给定简单的恒值控制，验证轨迹推演是否满足CW方程；然后再接入协态方程，跑完整的均衡求解。

3.2 相对运动离散化与仿真主循环

下面这段代码是仿真主循环的核心，用四阶Runge-Kutta对CW方程做时间推进。为了兼顾速度和精度，我通常取 dt = 0.01 个轨道周期单位，总步数取决于 t_f。

function [t, X, U] = simulate_trajectory(x0, uP_func, uE_func, n, tF, dt) % 基于四阶Runge-Kutta的末端追逃轨迹仿真 t = 0:dt:tF; N = length(t); X = zeros(6, N); U = zeros(6, N); X(:,1) = x0; for k = 1:N-1 tk = t(k); uP = uP_func(X(:,k), tk); uE = uE_func(X(:,k), tk); U(:,k) = [uP; uE]; k1 = relative_dynamics(X(:,k), uP, uE, n); k2 = relative_dynamics(X(:,k) + 0.5*dt*k1, uP, uE, n); k3 = relative_dynamics(X(:,k) + 0.5*dt*k2, uP, uE, n); k4 = relative_dynamics(X(:,k) + dt*k3, uP, uE, n); X(:,k+1) = X(:,k) + (dt/6)*(k1 + 2*k2 + 2*k3 + k4); end end

relative_dynamics里就是CW方程。注意uP和uE是绝对控制加速度在相对运动方程中的体现，所以差分项是uP - uE。

function xdot = relative_dynamics(x, uP, uE, n) % x = [rx; ry; rz; vx; vy; vz] xdot = zeros(6,1); xdot(1:3) = x(4:6); xdot(4) = 3*n^2*x(1) + 2*n*x(5) + uP(1) - uE(1); xdot(5) = -2*n*x(4) + uP(2) - uE(2); xdot(6) = -n^2*x(3) + uP(3) - uE(3); end

如果你用的是绝对轨道六根数或球形坐标，代码会复杂很多。但对末端追逃这个尺度，CW方程已经能抓住主要动态特征。仿真前记得做无量纲化处理，否则数值量级差异会让ODE45的步长控制非常痛苦。

3.3 基于间接法的两点边值问题求解

真正让博弈求解变得“硬核”的是哈密顿函数和协态方程。对于零和微分对策，哈密顿函数定义为：

H = λ^T * f(x, u_p, u_e) + L(x, u_p, u_e)

其中 f 是相对运动方程，L 是代价函数中的积分项。追击者最小化末值 + 积分，逃逸者最大化。根据庞特里亚金极小值原理，最优控制满足：

u_p* = argmin H, u_e* = argmax H

对于控制约束 ||u|| ≤ u_max，最优控制就是“方向由协态决定、幅值饱和”的结构。比如追击者：

u_p* = - u_{p,max} * λ_{u,p} / ||λ_{u,p}||

这里的 λ_{u,p} 是协态向量中与速度分量对应的部分。逃逸者方向相反。

将最优控制带回哈密顿函数，得到关于状态和协态的哈密顿系统：

dx/dt = ∂H/∂λ, dλ/dt = -∂H/∂x

再配上终端条件：

r(t_f) = r_f（自由或固定），λ(t_f) = ∂φ/∂x(t_f)

就构成了一个标准的两点边值问题。我用打靶法求解，代码框架如下：

function lambda0 = solve_tpbvp(x0, n, tF, uP_max, uE_max) % 猜测协态初值 lambda0_guess = zeros(6,1); opts = optimoptions('fsolve', 'Display','iter', 'TolFun',1e-10); lambda0 = fsolve(@(lam) terminal_residual(lam, x0, n, tF, uP_max, uE_max), ... lambda0_guess, opts); end function res = terminal_residual(lambda0, x0, n, tF, uP_max, uE_max) % 积分哈密顿系统到终端 % 这里需要把最优控制表达式内联到动力学里 % 返回终端状态与目标值的差 [~, X] = ode45(@(t,x) hamiltonian_dynamics(x, n, uP_max, uE_max), [0 tF], [x0; lambda0]); res = [X(end,1); X(end,2); X(end,3)] - 目标终端位置; end

实际调试时，fsolve经常给出“无解”或陷入局部极小。我的经验是：先用小规模随机搜索生成多个候选初值，再逐个精细打靶，最后挑选终端残差最小的那组协态作为均衡策略。

4. ε-纳什均衡的数值求解与策略迭代

4.1 并行哈密顿求解

上一节提到的最优控制表达式其实依赖于一个前提：追击者最小化的同时逃逸者最大化，双方同时求解哈密顿函数。这不同于单方最优控制——你不能先固定对方向后求自己的最优，因为对方也会同步改变策略。

在数值实现上，我把哈密顿系统写成一个函数，在每一时间步内先根据当前协态计算双方最优控制，再将控制代入状态方程和协态方程。这样能保证求解的是“同时决策”的博弈问题，而不是“序列决策”的主从问题。

function dydt = hamiltonian_dynamics(y, n, uP_max, uE_max) % y = [r; v; lambda_r; lambda_v] r = y(1:3); v = y(4:6); lambda_r = y(7:9); lambda_v = y(10:12); % 最优控制：方向取协态速度子向量方向，幅值饱和 uP = -uP_max * lambda_v / (norm(lambda_v) + 1e-8); uE = uE_max * lambda_v / (norm(lambda_v) + 1e-8); % 状态方程：CW方程 dydt(1:3) = v; dydt(4) = 3*n^2*y(1) + 2*n*y(5) + uP(1) - uE(1); dydt(5) = -2*n*y(4) + uP(2) - uE(2); dydt(6) = -n^2*y(3) + uP(3) - uE(3); % 协态方程：-∂H/∂x dydt(7:9) = -[3*n^2*y(10); 0; -n^2*y(12)]; dydt(10) = -2*n*y(11) - 3*n^2*y(9); % 这里需要仔细推导 dydt(11) = 2*n*y(10); dydt(12) = 0; end

上面的协态方程是我随手写的形式，真正使用前一定要从哈密顿函数重新推导一遍。数值上哪怕多了一个符号，轨迹都会完全走样。

4.2 策略迭代算法

严格意义上，同时求解哈密顿函数已经给出了均衡候选，但在数值上打靶法可能不收敛。我改用了一种更鲁棒的策略迭代算法：

1. 初始化：为追击者和逃逸者分别设置一组初始控制策略（比如恒值控制或者比例导引类策略）。
2. 固定逃逸者策略，追击者求解最优响应（单方最优控制问题）。
3. 固定追击者策略，逃逸者求解最优响应。
4. 重复步骤2-3，直到双方策略变化引起的终端距离变化小于给定阈值。
5. 对收敛后的策略组合做单方面扰动检验，计算ε值。

这个过程很像博弈论里的“虚拟博弈”或“迭代最优响应”。理论上，如果博弈是严格竞争的，迭代会收敛到一个近似鞍点解；如果出现振荡，说明可能存在多个均衡或者需要降低步长。

伪代码：

uP_old = @(x,t) zeros(3,1); uE_old = @(x,t) zeros(3,1); for iter = 1:20 % 固定逃逸者，求解追击者的最优控制 uP_new = solve_optimal_pursuer(uE_old); % 固定追击者，求解逃逸者的最优控制 uE_new = solve_optimal_evader(uP_new); % 检查策略变化是否小于阈值 delta = evaluate_strategy_change(uP_new, uP_old, uE_new, uE_old); if delta < 1e-4 break; end uP_old = uP_new; uE_old = uE_new; end

4.3 ε值怎么选

在仿真中，ε不是一个固定常数，而是一个后验计算的量级。计算方式：

1. 固定逃逸策略，让追击者单方面偏离当前策略，尝试若干扰动方向，记录终端距离的最小改善量 δ_p。
2. 固定追击策略，让逃逸者单方面偏离当前策略，记录终端距离的最大恶化量 δ_e。
3. ε = max(δ_p, δ_e)。

如果 ε 远小于终端距离的特征尺度，比如终端距离是500米，ε只有0.3米，那么我们就可以说当前策略是“工程意义上的均衡”。

我在仿真中发现，ε值与控制幅值、终端时间、策略迭代收敛阈值直接相关，参考下表：

这个表告诉我们：终端时间越长，博弈窗口越大，微小策略扰动对终端距离的影响越容易被放大，所以ε会增大。如果你希望ε保持在很小范围，要么缩短终端时间，要么把迭代收敛阈值设得更严。

5. 仿真结果分析与典型场景解读

5.1 同一轨道高度下的对称初始条件

我首先跑了一个“对称”场景：追击者与逃逸者初始相对位置和速度都对称，双方控制能力相同。这在理论上应该得到一个对称的均衡轨迹——双方在相对运动空间中划出镜像轨迹，终端距离处于中间水平。

仿真结果确实是这样的。追击者从初始点向逃逸者方向机动，逃逸者则沿近似镜像的方向“躲”。终端相对速度不大，双方在最后时刻都接近控制饱和。这个结果的意义在于验证模型：连对称场景都给出不对称轨迹，那大概率是代码有bug，而不是博弈本身有趣。

5.2 初始相位差对博弈结果的影响

轨道追逃和平面追逃最大的不同在于相对运动方程里的耦合项。初始相位差——也就是逃逸者在轨道的不同位置——会显著影响博弈结果。我设计了四组初始相位差，观察终端距离的变化：

有趣的是，初始相位差90°时终端距离反而最小。原因是此时的相对位置方向使得追击者可以利用CW方程中的科里奥利耦合进行“抄近道”，而逃逸者的最优规避方向受限于轨道几何，无法完全抵消追击者的优势。

这个结果也说明，末端追逃博弈不能只盯着追击者的最大加速度或逃逸者的初始速度看，初始相位差是一个容易被忽略却非常关键的因素。

5.3 控制能力不对称时的均衡移动

现实中追逃双方的控制能力很少相同。我调整了追击者的最大控制加速度从0.5到1.5 m/s²，保持逃逸者0.5不变。终端距离的变化曲线大致是单调递减的，但没有想象中那么快：

• 追击者0.5，终端距离约230米
• 追击者1.0，终端距离约160米
• 追击者1.5，终端距离约120米

控制能力翻倍，终端距离并没有减半。原因在于逃逸者在博弈均衡中也会调整策略，把有限的控制集中在威胁最大的方向上。换句话说，你更强，但对方也不是吃素的。

这个现象在单方最优控制里看不到。如果没有博弈均衡分析，很容易给出“追击者推力翻倍，捕获时间减半”这种过度乐观的结论。

6. 复现过程中踩过的坑与可复用建议

6.1 打靶法初值敏感与归一化处理

第一次跑打靶法时，我直接用了物理单位下的状态和协态。结果数值量级差距大到离谱：位置几百米，速度米每秒，协态可能到1e-5甚至更小。ODE45的绝对和相对误差容限很难同时满足，终端残差也像过山车一样忽大忽小。

解决办法是先做无量纲化。把长度除以轨道半径特征尺度，时间除以轨道周期特征尺度，控制加速度除以特征加速度。这样所有状态量都落在0.1到10的范围内，打靶法的收敛性立刻得到改善。

具体来说，我建议至少做三件事：

• 状态与协态同时归一化，保证数量级一致。
• 对控制饱和函数做光滑化，比如用双曲正切替代硬饱和。
• 打靶变量用对数缩放或符号变换，避免非负约束。

6.2 离散化步长与终端约束的匹配

仿真主循环我一开始用固定步长dt = 0.1，结果终端距离对步长的敏感性极高。有一次把步长从0.1改到0.02，终端距离直接从140米跳到190米，差点以为博弈重新出现了新均衡。

后来才反应过来：CW方程里有2n的科里奥利项，时间尺度和轨道周期挂钩。固定步长如果大于轨道周期的1/50，数值积分误差就会明显吃掉动态特性。我的建议是：要么用ode45这类自适应步长积分器，要么在仿真前做一次步长收敛性分析，画出终端距离随步长变化的曲线，选择平台区对应的步长。

6.3 从二维平面到三维空间的扩展思路

很多文献示例先从二维平面问题讲起，我也这么做了，但真要写三维仿真时发现事情没那么简单。三维CW方程的径向、迹向和法向动态是解耦的？表面上看，z方向确实独立，但实际博弈中控制方向在三维空间里会带来更丰富的轨迹形态，尤其是逃逸者可以从垂直轨道面方向“逃跑”，追击者必须权衡是沿迹向追赶还是先改变轨道面。

扩展时最省力的方法是把2D的哈密顿系统直接改成6维状态 + 6维协态，控制向量变成三维。但注意，协态与控制方向的关系在高维下有可能出现“方向抖动”，尤其是在最优控制接近奇异弧时。我的做法是在控制更新时加入低通滤波，避免相邻步之间控制方向突变。

6.4 基于个人体验的仿真调试建议

最后分享几点调试心得，可能比上面的公式更有用。首先，不要一上来就跑完整博弈，先做一个“单边最优”对照组。固定逃逸者策略为零，让追击者做纯最优控制，看终端距离是多少；再固定追击者为零，让逃逸者做纯最优规避。这两个单边结果可以作为合理区间参考，均衡结果应该介于两者之间，否则说明策略迭代出了问题。

其次，每次修改完模型，都要重新跑一次对称场景。对称场景的结果如果不对称，首先要怀疑fsolve求出的协态初值是局部最优而非全局最优。遇到这种情况，我的土办法是用多组随机初值竞争，选终端残差最小且残差随迭代次数稳定下降的那组。

最后，ε-纳什均衡的检验不能只看数值，还要结合仿真曲线观察。我会画出多个时间剖面上的控制方向与相对位置关系，看看是否违背基本的物理直觉。比如追击者明明可以直飞接近，却绕了一大圈，这通常是哈密顿函数推导或协态方程符号错误导致的。

在跑完这轮仿真之后，我最大的体会是：基于ε-纳什均衡的追逃仿真，本质上是在“严格博弈”和“数值可实现”之间寻找平衡点。理论模型保证了问题的博弈论意义，ε松弛保证了数值求解的鲁棒性，而Matlab实现里的每一次调参、每一个数值技巧，都是在跟这种平衡反复磨合。希望这篇文章能帮你少走一些弯路，把精力集中到真正有意思的博弈现象分析上。