基于ε-纳什均衡的航天器末端追逃博弈Matlab仿真实现

前两天终于把基于ε-纳什均衡的航天器末端追逃博弈仿真完整跑通了。整个过程比我预想的曲折理论上看微分对策和纳什均衡的框架很清晰但一落到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 固定。这个设定与现实任务一致追逃窗口由轨道预报决定超过窗口后的事件不再纳入博弈模型。有限时间特性给数值求解带来两个影响协态变量的终端条件可以直接由终端代价函数的梯度给出不需要横截条件的复杂推导。哈密顿函数在终端附近可能出现奇异性尤其是控制约束起作用时bang-bang控制需要做光滑化处理。我把这些特征全部纳入后续的求解框架。模型越贴近实际数值求解的“脾气”越大但只要每一步都有明确物理含义调试起来还是能找到方向的。3. Matlab仿真架构设计与核心代码实现3.1 模块划分从初始化到后处理写仿真代码最忌讳把整个博弈过程塞进一个几千行的脚本里。我的做法是把工程拆成四个模块初始化模块设置轨道参数、控制约束、终端时间、玩家策略初值。博弈求解模块核心策略迭代与TPBVP求解。仿真推进模块用离散化的相对运动方程做开环或闭环轨迹仿真。后处理模块计算终端距离、绘制双方轨迹、做单方面策略扰动检验。这种划分方式让我在调试时可以单独测试每个模块。比如先不跑策略迭代直接给定简单的恒值控制验证轨迹推演是否满足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(:,k1) X(:,k) (dt/6)*(k1 2*k2 2*k3 k4); end endrelative_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 策略迭代算法严格意义上同时求解哈密顿函数已经给出了均衡候选但在数值上打靶法可能不收敛。我改用了一种更鲁棒的策略迭代算法初始化为追击者和逃逸者分别设置一组初始控制策略比如恒值控制或者比例导引类策略。固定逃逸者策略追击者求解最优响应单方最优控制问题。固定追击者策略逃逸者求解最优响应。重复步骤2-3直到双方策略变化引起的终端距离变化小于给定阈值。对收敛后的策略组合做单方面扰动检验计算ε值。这个过程很像博弈论里的“虚拟博弈”或“迭代最优响应”。理论上如果博弈是严格竞争的迭代会收敛到一个近似鞍点解如果出现振荡说明可能存在多个均衡或者需要降低步长。伪代码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; end4.3 ε值怎么选在仿真中ε不是一个固定常数而是一个后验计算的量级。计算方式固定逃逸策略让追击者单方面偏离当前策略尝试若干扰动方向记录终端距离的最小改善量 δ_p。固定追击策略让逃逸者单方面偏离当前策略记录终端距离的最大恶化量 δ_e。ε max(δ_p, δ_e)。如果 ε 远小于终端距离的特征尺度比如终端距离是500米ε只有0.3米那么我们就可以说当前策略是“工程意义上的均衡”。我在仿真中发现ε值与控制幅值、终端时间、策略迭代收敛阈值直接相关参考下表控制幅值m/s²终端时间s迭代收敛阈值实测εm0.5301e-30.120.5601e-30.851.0301e-30.211.0601e-40.42这个表告诉我们终端时间越长博弈窗口越大微小策略扰动对终端距离的影响越容易被放大所以ε会增大。如果你希望ε保持在很小范围要么缩短终端时间要么把迭代收敛阈值设得更严。5. 仿真结果分析与典型场景解读5.1 同一轨道高度下的对称初始条件我首先跑了一个“对称”场景追击者与逃逸者初始相对位置和速度都对称双方控制能力相同。这在理论上应该得到一个对称的均衡轨迹——双方在相对运动空间中划出镜像轨迹终端距离处于中间水平。仿真结果确实是这样的。追击者从初始点向逃逸者方向机动逃逸者则沿近似镜像的方向“躲”。终端相对速度不大双方在最后时刻都接近控制饱和。这个结果的意义在于验证模型连对称场景都给出不对称轨迹那大概率是代码有bug而不是博弈本身有趣。5.2 初始相位差对博弈结果的影响轨道追逃和平面追逃最大的不同在于相对运动方程里的耦合项。初始相位差——也就是逃逸者在轨道的不同位置——会显著影响博弈结果。我设计了四组初始相位差观察终端距离的变化初始相位差°终端距离m终端时刻双方相对速度m/sε值m0214.32.80.3145182.53.60.4590158.24.50.52135236.72.20.38有趣的是初始相位差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实现里的每一次调参、每一个数值技巧都是在跟这种平衡反复磨合。希望这篇文章能帮你少走一些弯路把精力集中到真正有意思的博弈现象分析上。