ε-纳什均衡下的航天器末端追逃博弈：Matlab策略迭代仿真与实现

最近组里有朋友在复现“基于ε-纳什均衡的航天器末端追逃博弈”这类论文卡在收敛那块好几天。我在Matlab里也把整套仿真完整跑通过前前后后改了七版模型踩了各种坑最后才把基于ε-纳什均衡的追逃博弈策略迭代稳定下来。说实话这个方向看起来高大上真正动手做时难点不在博弈论公式而在模型简化、收益函数设计、以及离散化后迭代过程的数值稳定性。这篇就当一份复现笔记把末端追逃的动力学建模、ε-纳什均衡为什么会被引入、策略迭代的主流程、以及Matlab代码实现里哪些地方最容易出bug都梳理一遍给正在复现类似论文的读者做个参考。1. 追逃博弈模型先把末端运动方程写对1.1 视线坐标系下的相对运动建模航天器末端追逃博弈本质上是一个双方控制问题一方称为追击者记作P要尽量接近目标另一方称为逃避者记作E要尽量远离。末端段的特征很明显相对距离通常只有几公里到几十公里交战时间在几十秒量级此时轨道摄动、地球扁率这些因素对相对运动的影响很小所以工程上习惯把相对动力学简化到视线坐标系里去处理。我用的模型是这样的。在二维追逃平面内取追击者P与逃避者E之间的相对距离为r视线角为q那么相对运动可以写成[ \ddot{r} - r\dot{q}^2 a_{E,r} - a_{P,r} ][ r\ddot{q} 2\dot{r}\dot{q} a_{E,q} - a_{P,q} ]其中(a_{P,r})、(a_{P,q})分别是追击者在视线径向和法向的推力加速度分量(a_{E,r})、(a_{E,q})同理。这个方程组就是极坐标下的相对运动方程也是后续所有博弈计算的底层动态。需要注意这里默认双方推力连续可调且没有考虑执行器延迟属于理想控制量假设。如果要引入延迟状态量还得再加一拍迭代复杂度会明显上升所以复现论文时一般先不加。有些论文里会用惯性坐标系下的相对位置和相对速度做状态我试过之后觉得还是视线系的表达更贴合末端交会的物理直觉。因为视线角速率(\dot{q})本身就是追踪制导里经常关注的一个量用来刻画接近过程中的旋转程度。状态量我取为([r, q, \dot{r}, \dot{q}])这样收益函数可以很方便地对视线角速率做惩罚也方便后续和比例导引这类经典制导律做对比。末端追逃的整个数值博弈就是这个四维状态空间里的策略搜索问题。1.2 博弈三要素到追逃问题的映射一个博弈问题要有明确的参与者、策略集、收益函数。追逃问题里这三件套非常清晰实际建模时只需要把它们一一对应起来。参与者就是追击者P和逃避者E二者目标冲突这是双人零和微分博弈的标准场景。策略集是指每个时刻双方根据自己的观测状态选择加速度指令。为了数值求解我事先把每个控制分量离散成有限集合例如[−2, −1, 0, 1, 2] m/s²径向和法向各取一组于是每个参与者每一步大约有25种可选指令组合。控制集大小直接影响策略迭代的计算量后面我会专门讲怎么折中。收益函数我采用了终端收益加过程收益的组合方式。终端收益以末端脱靶量为主过程收益考虑燃料消耗和视线角速率的积分。具体来说追击者的收益函数设计为[ J_P \Phi(x(t_f)) \int_{t_0}^{t_f} L(x(t), u_P(t), u_E(t)) dt ]其中(\Phi(x(t_f)) \alpha r^2(t_f) \beta \dot{q}^2(t_f))(\alpha)和(\beta)是权重系数过程项(L \lambda_P |u_P|^2 \lambda_E |u_E|^2 \gamma \dot{q}^2)。数值上我让追击方最小化这个收益逃避方最大化。实际写代码时也可以定义双方各自的收益函数让追击方取负号关键是符号约定要全程一致。这里最需要注意的点收益函数设计直接决定了均衡策略长什么样。如果把视线角速率权重(\gamma)调得很大追击方的策略会明显偏向平滑接近而不是快速趋近从博弈结果来看就是末端过载小但接近时间变长。一开始我直接把某一篇论文里的权重照搬过来结果跑出来的轨迹特征完全对不上后来才发现权重必须跟状态量纲、控制量范围联动调整不能盲目套。建议的做法是先固定控制量范围把终端距离项归一化到某个量级比如让初始距离对应的终端收益为10左右再从这里出发去调其他项的权重。2. 为什么末端追逃要用ε-纳什均衡2.1 精确纳什均衡在连续微分博弈中的求解困境如果能在时间连续、状态连续的意义上求得双方的最优策略那就得到了经典微分博弈的鞍点解也就是需要求解HJI方程全称Hamilton-Jacobi-Isaacs方程。对于双人零和微分博弈我们需要解一个偏微分方程[ \frac{\partial V}{\partial t} \min_{u_P}\max_{u_E} \left( \nabla_x V \cdot f(x, u_P, u_E) L \right) 0 ]边值条件由终端收益给出。这个方程有解析解的场景非常少末端追逃问题里因为视线动力学耦合、控制约束基本没有闭式解。数值解HJI方程也确实可行但需要在状态空间的每个格点上同时求解偏微分方程维度一高计算量立刻爆炸。我一开始也直接尝试HJI数值解法结果四维状态只要把网格稍微细化一点Matlab内存就吃紧等一次迭代要十几分钟还伴随着数值耗散带来的策略抖动。这里面的矛盾在于精确纳什均衡要求每一步都找到全局最优响应而连续状态空间的全局最优响应本身就是一个高维优化问题。在末端追逃这种短时间强耦合的场景里这个代价很不划算。更现实的问题是我们最终的输出是离散化的策略表控制指令本身也是从有限集合里选的。既然控制已经是离散的状态也是格点化的那追求一个理论上的精确均衡就没有太多工程意义。能够找到一个策略组合让任何一方单方面改变策略都得不到显著收益提升就已经足够支撑后续的制导逻辑设计和效能分析。2.2 ε-纳什均衡的数学定义与工程价值正式一点对双人博弈如果策略组合((u_P^, u_E^))满足对追击方的任意可选策略(u_P)有[ J_P(u_P, u_E^) \geq J_P(u_P^, u_E^*) - \varepsilon ]对逃避方的任意可选策略(u_E)有[ J_E(u_P^, u_E) \geq J_E(u_P^, u_E^*) - \varepsilon ]就称这个策略组合是ε-纳什均衡。注意符号约定取决于收益是最大化还是最小化我这里写成两边都“不会因为单方偏离而损失超过ε”实际编程时按这个逻辑实现即可。用大白话说就算某一方偷偷换成别的策略最多也只能额外占一点便宜而且这个便宜的量不超过ε。这在工程场景里非常实用因为我们离散化了控制量、截断了状态空间这些近似操作本身就会引入误差追求ε等于零没有任何实际意义。只要ε小到某个阈值比如收益量级的1%这个策略组合就可以放心拿去用。用一个生活类比来理解两个人轮流报价抢一个标的如果无论对方怎么换报价策略你的收益都不会变差超过一个很小的数值那这个报价策略组合就算稳定了。虽然理论上可能还存在着略好一点的策略但找它的成本太高不如接受这个接近最优的结果。我在程序里做的判定很简单每次策略迭代完成后固定其中一方策略遍历另外一方所有离散控制组合计算最大收益偏差如果双方的偏差都小于预设的ε就判定收敛。这个过程中涉及的收益计算需要靠前向推演加后向回溯完成也就是下一章要展开的核心内容。3. Matlab仿真核心流程从策略迭代到收敛判定3.1 状态空间离散与控制集设计我在Matlab里把整个求解过程做成了一套策略迭代流程。第一步是对状态空间做离散化。我实际使用的网格参数如下相对距离r取0到12 km等分为50格视线角q取−60度到60度等分为40格相对速度ṙ取−80到20 m/s等分为30格视线角速率q̇取−0.15到0.15 rad/s等分为20格。这样状态总格点数约为50乘以40乘以30乘以20也就是120万个。四维状态直接存储策略矩阵的话每个策略矩阵是120万个double数值双方各存一份内存占用大概20MB左右在Matlab里还能接受。再往细了加网格计算量就明显不可控了所以这个规模算是一个比较合适的折中。控制集我设为径向和法向各5档范围是−2到2 m/s²。两个通道合起来每个玩家有25个可选控制组合。终端时间我取为60秒如果状态提前越界就按越界时刻结算。积分步长取0.01秒用固定步长RK4。这些参数的选取直接影响策略推演的精度和耗时如果积分步长过大末端收益会出现明显误差步长过小则迭代时间成倍增加。3.2 策略迭代与收益回溯过程整个迭代的大循环大概长这样初始化追击方和逃避方的策略矩阵可以全部置零也可以用比例导引启发式给一个初值固定逃避方策略对每个状态格点遍历追击方所有控制组合选择目标函数最小的控制更新追击方策略矩阵固定追击方策略对每个状态格点遍历逃避方所有控制组合选择目标函数最大的控制更新逃避方策略矩阵计算本轮与上轮策略矩阵的变化量同时做一次ε-纳什均衡判定若满足收敛条件则终止否则回到第2步。收益推演不是只走一步而是要走完整个末端交战过程直到终端时间或状态超出边界。也就是说对每个状态格点在选定双方控制后要用数值积分把相对运动推演几十秒累加过程收益并计算终端收益才能得到这一步控制对应的收益值。这个步骤是整个程序最耗时的部分也是优化空间最大的地方。一个非常关键的经验收益计算最好做成向量化形式不要在每个状态格点上单独循环。我早期写的是逐点循环四维状态下跑三百轮迭代大概要两三个小时。后来把状态格点全部展平成矩阵一次性做状态转移和收益计算同样的迭代轮数压缩到二十分钟以内差距非常明显。下面的Matlab伪代码展示了核心迭代骨架省略了具体收益函数计算细节但结构足够说明问题% 初始化策略矩阵 strategy_P zeros(nr, nq, nrdot, nqdot); strategy_E zeros(nr, nq, nrdot, nqdot); % 初始化收益矩阵 J_table zeros(nr, nq, nrdot, nqdot); for iter 1:max_iter % 固定E策略更新P策略 for idx 1:numel(strategy_P) [r, q, rd, qd] ind2sub(size(strategy_P), idx); best_u []; best_J inf; for uP control_set_P J_val compute_payoff(r, q, rd, qd, uP, strategy_E(idx)); if J_val best_J best_J J_val; best_u uP; end end new_strategy_P(idx) best_u; end % 固定P策略更新E策略类似上面但取最大化 % 判断策略变化量和epsilon均衡条件 % 满足则break end实际写代码时compute_payoff函数必须写成支持矩阵输入的形式否则性能会很差。如果读者刚开始接触这个方向建议先按最清晰的逐点循环写通再逐步改成向量化不要一上来就追求性能否则bug会很难排查。3.3 收敛参数设置与效果对比收敛判据方面我用了两个条件同时判断。第一个是策略变化量本轮与上一轮策略矩阵中取值不同的格点占比小于某个阈值比如0.1%。第二个是ε-纳什均衡检验固定对手策略分别做单方偏离测试最大收益增量小于预先给定的ε。ε的取值推荐结合收益量级来定。比如终端收益项大约几百到几千的量级那ε取0.5到1就足够了如果收益函数整体量级很小ε过大会导致策略粗糙过小则迭代很难收敛。我尝试过不同参数组合结果如下表网格规模控制档位数迭代轮数单轮耗时收敛ε30×20×15×10587约4秒1.050×40×30×205156约30秒0.870×55×40×257230以上约2分钟0.5未完全收敛表中的数据来自我某次实际运行记录可以清楚看到网格规模和收敛速度的矛盾。第一次运行用最粗的网格收敛很快但画出来的策略云图有很多锯齿细化后策略变得平滑但迭代时间明显拉长。实际上对于系统分析50×40×30×20配合5档控制已经够用。如果想拿去做实时制导应用这个求解规模仍需进一步压缩比如用查表加降维的方法那是另一个话题。4. 复现结果分析均衡策略到底长什么样4.1 不同初始条件下的均衡路径对比跑通迭代之后最有价值的工作是分析均衡状态下的相对运动轨迹。我选了三种典型初始状态做对比初始状态相对距离视线角径向速度视线角速率工况A8 km20度−50 m/s0.08 rad/s工况B5 km−15度−30 m/s−0.05 rad/s工况C10 km5度−70 m/s0.02 rad/s在ε-纳什均衡策略下三条轨迹有一个共同特征前期双方都在积极试探相对距离下降不快当相对距离缩小到某一阈值后追击方的径向推力开始加大而逃避方则尽量把视线角速率往大里拉。因为视线角速率越大追击方的法向通道就越难把视线稳住末端过载需求越大。有一个现象值得注意在工况B里初始视线角是负的逃避方最优策略是先往正视线角方向甩头而不是直接远离。一开始我看这个结果总觉得不对后来画了控制量时间曲线才明白逃避方这么做是通过改变视线角来增加追击方法向通道的负担牺牲一部分径向逃逸速度来换取更大的整体收益。这个思路在纯工程直觉里很容易被忽略但博弈解会自然地展现出这类策略。我把这个结果和用比例导引跑出来的轨迹放一起对比差异非常明显比例导引下追击方会一路把视线角速率压到接近零轨迹呈现经典的直线接近特征而博弈均衡策略下追击方的轨迹有明显的横向扰动因为逃避方在持续改变角度追击方必须不断修正视线方向末端法向过载曲线波动也更大。4.2 ε-纳什均衡策略的解释粒度另一个值得说的点是ε-纳什均衡策略对初始状态敏感但又不是完全敏感。具体表现为如果两个初始状态在状态空间里的距离很小它们的均衡策略在高维空间里依然会比较接近一旦初始状态跨越了某条策略分界线比如径向速度由负变正策略结构会发生整体切换。这其实就是博弈策略的bang-bang特性而离散控制集的引入让这种现象更明显。我的解释是在离散控制集且游戏时长有限的前提下双方策略从本质上更接近多阶段最优响应而ε-纳什均衡刻画的就是这种数值意义上的稳态。实际复现的时候不需要对结果做过度理论解读重点是把轨迹、控制曲线和收益值三者对上账。如果轨迹合理但收益数值很奇怪通常问题出在收益函数某一项的权重没有归一化而不是求解过程有错。5. 踩坑记录与常见问题排查5.1 策略不收敛的几个典型原因我在复现过程中遇到最典型的坑有三个。第一个是初始策略设置。如果不加区分地全部置零某些状态下追击方和逃避方谁也不动收益完全由终端状态决定迭代很容易陷入局部近似解。后来我改用启发式初值追击方初始策略先按视线角速率零点为目标逃避方按使视线角速率绝对值增大的方向设收敛速度和最终收益值都明显改善。第二个是控制离散集过于稀疏。5档控制如果覆盖范围太大比如±5 m/s²容易让收益曲面在离散点上出现多峰导致策略迭代在两个状态之间震荡表现在结果里就是相邻轮次策略矩阵变化率始终降不下来。解决办法是缩小控制量范围或者先粗后细地做两轮求解。先用大步长探索大概方向再缩小范围细化这样能有效避免震荡。第三个是终端时间一致性。有一个版本我把双方的最大仿真时长设成一样的但追击方如果提前达到极小距离就会提前终止导致两个玩家的游戏视野不一致迭代后期出现收益突变。修正方案是统一按先触界先结束的方式处理双方同时结算收益保证评估口径一致。5.2 网格选择与Matlab计算性能的平衡Matlab跑这个模型的瓶颈在于收益函数的批量推演因为每个状态格点都要完成RK4积分。如果整成双层循环运算速度会非常感人。我推荐的做法是把状态格点按列组织成矩阵控制指令也用同样的矩阵形式广播这样一次状态更新就变成矩阵乘法收益函数的累加也能用sum这类向量化函数完成。另外要注意内存预分配。策略矩阵、收益矩阵在迭代前就初始化成固定大小不要在循环内动态扩维。我最初没做预分配Matlab在循环里频繁重新分配内存一百轮迭代之后速度明显下降改成预分配后立刻恢复。还有个细节是尽量把不会变的参数比如控制集、网格边界、积分步长放到循环外面避免重复计算。5.3 代码实现中的细节提醒最后列几个容易忽略的细节每一个都是我实际踩过的坑角度的单位统一。状态里如果用角度制来表示视线角那么速度方程里的角速率项必须转回弧度制否则数值上会出现一个量级差10倍以上的隐性bug。我在第一版就吃了这个亏排查了很久才发现是单位混用。边界状态的处理。相对距离r接近零时相对动力学方程里会有除以r的项要设置一个最小距离限制或者直接在这个状态下跳出推演并按终端收益结算。我把最小距离设成0.1 km低于这个值就认为已满足追及条件。收益函数归一化。因为终端收益和过程收益的量纲不同直接把几个项加起来会让某个项主导建议把每个项除以一个参考值再乘以权重系数。这样调试权重时更直观也更容易复现论文中的结果。下面整理了一份常见问题速查表方便排查现象可能原因解决思路策略矩阵长期不更新初始策略全零双方都不动改用启发式初值迭代曲线震荡不收敛控制离散集过疏或多峰缩小控制范围分两轮细化轨迹末端收益突变双方终端时间不一致统一先触界先结束规则法向控制量异常大视线角速率权重过小增大γ系数或做归一化收敛很快但轨迹很奇怪网格太粗细化状态空间网格程序越跑越慢循环内动态扩维预先分配矩阵内存我个人在实际操作中的体会是这套仿真真正值钱的部分在于把均衡策略可视化之后那些很反直觉的机动轨迹它逼着你重新理解追逃对抗的本质。代码本身不难难的是一遍遍调整收益函数和网格参数去逼近稳定的ε-纳什均衡过程。跑通之后如果有余力可以把这个框架往下游扩展比如加上终端角约束、饱和控制约束甚至是多追击者协同追逃底层的策略迭代骨架基本不用大改这也是我觉得这个方向值得花时间复现的原因。