1. 项目背景与研究价值追逃博弈到底在解决什么问题先把这个标题拆开看。基于ε-纳什均衡的航天器末端追逃博弈仿真本质上是在做这么一件事在航天器交会对接、逼近拦截这类末端场景里用微分博弈的语言描述“追的一方”和“逃的一方”各自该怎么做然后通过数值方法求出双方都能接受的一组策略再用Matlab把整个过程跑出来看效果。这类问题我在制导控制领域见得不少但它远不只是航天方向的专属话题。追逃博弈本身是一个经典的双人零和微分博弈模型逃方想让两者距离维持在安全范围之外追方则想尽量缩小甚至消除这个距离。难点在于这不是单边优化问题而是双方都在实时决策、相互对抗的闭环博弈。末端场景的特殊性在于时间窗口很短、相对运动状态变化剧烈而且双方的控制能力都有限比如追方有视线角约束逃方有可用过载限制所以简单套用LQR或者比例导引这类传统制导律往往不够用必须把“对手会怎么做”显式地建模进去。既然要复现这个项目就得回答几个核心问题什么是ε-纳什均衡跟理想纳什均衡有什么区别在航天器追逃场景下双方的数学模型和支付函数怎么设计怎么用Matlab数值求解出一组满足ε-均衡条件的策略得到仿真结果之后怎么验证这组策略真的是均衡而非随意挑的一组“看起来合理”的策略适合看这篇文章的人我认为有三类。第一类是正在做航天制导、飞行器博弈方向的研究生这个复现过程能帮你们把微分博弈从公式推导落到可运行的代码上。第二类是学习现代控制理论、想了解非线性最优控制和博弈论如何结合的人这里面的方程推导和数值解法值得推敲。第三类是单纯对Matlab数值仿真感兴趣、想看看如何用bvp4c、fsolve这类工具求解大规模边值问题的工程师这个项目的求解框架很有参考价值。2. 数学模型拆解从物理场景到可计算的博弈2.1 追逃双方的运动学建模与坐标系选择建模的第一步是确定坐标系。末端追逃场景通常有两种选择惯性坐标系和相对运动坐标系视线坐标系。我做这个项目复现时选的是视线坐标系加线性化相对运动方程原因有三一是末端场景里追逃双方的距离相对较小在视线方向上做线性化误差可控二是后续的博弈论推导需要状态方程尽量简洁线性形式可以避免非线性带来的均衡求解困难三是工程上视线坐标系直观脱靶量、视线角速率这些量可以直接作为制导信息使用。设相对位置矢量为(r)相对速度矢量为(v)追逃双方的控制加速度分别为(u_p)和(u_e)运动方程写成[ \dot{r} v ][ \dot{v} u_e - u_p ]这里必须注意符号约定。追方控制量(u_p)作用于相对加速度的负方向意思是追方通过施加自身加速度来减小相对位置变化逃方控制量(u_e)刚好反过来是往增大相对距离的方向施加加速度。实际实现中我把双方的最大加速度限制分别写成(A_p^{max})和(A_e^{max})这直接对应了导弹拦截里“可用过载”的概念。如果完全不做约束求解出来的控制量在实际飞行器上是飞不出来的这点在后面的数值求解中也要处理。状态向量我取为([x, y, v_x, v_y]^T)分别是视线坐标系下的两个方向的相对位置和相对速度。之所以没有直接把脱靶量作为终端条件约束进微分方程而是放在支付函数里是为了让博弈双方在最小化自身代价的同时自然产生终端趋势这个设计思路跟一般的轨迹优化略有不同。2.2 支付函数设计与ε-纳什均衡的定义博弈论里的支付函数决定了双方“在乎什么”。这个项目里追方和逃方的支付是零和的即一方的收益就是另一方的损失。我把支付函数写成终端脱靶量的平方加控制能量的积分惩罚[ J \frac{1}{2} r(t_f)^T Q_f r(t_f) \frac{1}{2} \int_{t_0}^{t_f} \left( u_p^T R_p u_p - u_e^T R_e u_e \right) dt ]解释一下为什么要这样设计。终端项(\frac{1}{2} r(t_f)^T Q_f r(t_f))要求追方在末端时刻把相对位置压到最小对应拦截或交会的核心目标。控制能量项则是对双方策略的“经济性”约束——如果控制完全不受惩罚数值求出来的控制量会非常大这在工程上没有意义。逃方控制项前面的负号体现零和博弈的本质逃方消耗同样的控制量对支付函数的贡献方向相反它想最大化(J)而追方想最小化(J)。这里顺带说一个细节很多刚接触这个项目的人会问既然逃方要最大化J那为什么积分项里逃方的控制还是以惩罚形式出现答案在于这是一个博弈问题不是单纯的最优控制问题。逃方在最大化可以由它控制的那部分支付的同时也要考虑追方的最优反应。把逃方控制项写为负号惩罚是对称地构造博弈结构的自然做法跟单方最优化中“惩罚就是代价”的直觉并不冲突。理想纳什均衡要求双方中的任何一方单独偏离策略都不能获得更优的支付。但这个定义在数值求解中有个麻烦如果没有严格等于那就完全不是均衡而在连续状态空间中由于数值精度和离散误差精确的纳什均衡往往很难求出。ε-纳什均衡把这个条件放宽为任何单方面偏离的策略最多只能让自己的支付改善不超过ε。这个ε在实际仿真中就充当了一个松弛量让我们能在数值可解的范围内找到“足够接近均衡”的策略组合。我给项目设定的ε值范围在(1 \times 10^{-3})到(1 \times 10^{-1})之间具体取值取决于初值猜测质量和求解精度。当ε取得过大比如大于0.5得到的策略实际上跟随机策略没什么区别均衡验证也形同虚设取得过小比如低于(1 \times 10^{-4})边值问题的残差又很难收敛到这个容差以下。这个尺度需要在实际仿真中反复试出来。3. 求解路径与算法选型不是所有均衡都能直接“算出来”3.1 最优性条件Hamiltonian函数与共态方程求解博弈均衡和求解单个最优控制问题的思路一脉相承。我对系统的Hamiltonian函数做一阶变分再结合极值条件就能推导出双方的最优策略形式。定义Hamiltonian[ H \frac{1}{2} \left( u_p^T R_p u_p - u_e^T R_e u_e \right) \lambda^T \left( \begin{bmatrix} 0 I \ 0 0 \end{bmatrix} x B_e u_e - B_p u_p \right) ]其中(\lambda)是共态向量。极值条件要求(\frac{\partial H}{\partial u_p} 0)和(\frac{\partial H}{\partial u_e} 0)直接解出双方控制量与共态向量的关系[ u_p^* R_p^{-1} B_p^T \lambda ][ u_e^* -R_e^{-1} B_e^T \lambda ]注意这里(u_p)和(u_e)的表达式只差一个符号和权重矩阵这正是零和博弈在最优性条件上的镜像对称性。把控制量代回状态方程和共态方程就得到了一个关于状态(x)和共态(\lambda)的两点边值问题。状态从(t_0)出发初始条件(x(t_0))给定共态在终端满足横截条件(\lambda(t_f) Q_f r(t_f))。这个边值问题没有解析通解必须用数值方法求解。3.2 数值解法间接法与直接法的取舍求解这个边值问题的数值方案有几种间接法打靶、多重打靶或者直接做配点法。我在复现时最开始用的是单重打靶法思路很直白猜测一个初始共态(\lambda(t_0))正向积分状态和共态方程然后检查末端条件(\phi(\lambda(t_0)) \lambda(t_f) - Q_f r(t_f))是否为零。非零就迭代修正(\lambda(t_0))直到残差足够小。单重打靶实现简单但缺点也很明显对初值猜测极度敏感尤其在末端时刻较长、系统矩阵有正特征值的情况下微小的初值误差会被指数放大出现“面目全非”的终端状态。这是我第一次跑通之前踩得最深的坑。后来换成多重打靶才稳定下来——把整个时间弧段切成多个子区间每个子区间内分别积分然后强制相邻区间的状态和共态连续。如果不想自己写打靶逻辑也可以直接用Matlab自带的bvp4c或bvp5c求解器。bvp4c本身就是一个配点法实现配合bpset函数设置误差容差能省掉很多底层麻烦。但我发现bvp4c在强非线性或较大时间尺度下偶尔会收敛到物理上不合理的解这时候还是得回到打靶法自己控制迭代过程。对于ε-纳什均衡的验证还要再做一步额外的计算把所有可选的扰动态策略离散化逐一计算偏离支付然后统计最大偏离量是否超过ε。这一步理论上是必须的否则你只能说“我解出来的控制量满足极值条件”不能声称“这是一个ε-纳什均衡”。4. Matlab仿真实现与关键代码逻辑4.1 代码整体结构与数据流整个Matlab工程我按功能拆成五个文件主脚本、动力学函数、代价函数、求解器包装、绘图脚本。因为这一行代码解决的问题比较固定所以我尽量不搞成堆在一个文件里的面条代码下面是我在复现时用的文件组织方案% main_simulation.m - 主脚本负责参数初始化与流程调度 % dynamics.m - 相对运动微分方程与共态方程右端项 % cost_function.m - 支付函数计算 % solve_sweep.m - 打靶法/边值问题求解入口 % plot_results.m - 结果可视化主脚本的流程可以用几句话概括先定义双方最大加速度、初始相对位置、初始相对速度、终端时间范围和ε阈值然后调用求解器解出均衡策略再分别固定一方策略、让另一方做有界扰动态寻优计算偏离支付最终把状态轨迹和控制量剖面画出来并把脱靶量和策略变化轨迹放进结果报告中。4.2 动力学函数与边值条件实现动力学函数是整个仿真的核心。状态方程是线性相对运动但共态方程是耦合的必须同时积分。我这里给出一个典型的右端项代码片段关键部分是确保输入输出维度和索引顺序一致这一行出了问题会浪费很多调试时间function dydt dynamics_rhs(t, y, params) % 状态变量索引: y(1:4) - 相对状态, y(5:8) - 共态变量 x y(1:4); lambda y(5:8); % 解出双方最优控制 u_p params.inv_Rp * params.Bp * lambda; u_e -params.inv_Re * params.Be * lambda; % 状态方程 dx/dt A*x Be*u_e - Bp*u_p dxdt params.A * x params.Be * u_e - params.Bp * u_p; % 共态方程 dlambda/dt -dH/dx dlambda_dt -params.A * lambda; dydt [dxdt; dlambda_dt]; end下面是主脚本里的仿真参数初始化。这些数值不能随便拍脑袋定最好先用量纲分析确认一下单位一致性。我做的是无量纲化处理把距离归一化到初始相对距离时间归一化到总博弈时长。好处是矩阵元素数值都在同一量级避免出现跨越十几个数量级的数值运算边值问题求解收敛性也会好很多。% 仿真参数设定 params.A [0 0 1 0; 0 0 0 1; 0 0 0 0; 0 0 0 0]; params.Bp [0 0; 0 0; 1 0; 0 1]; params.Be -params.Bp; % 相对运动方程中逃方加速度方向与追方相反 params.Qf 8 * eye(2); % 终端位置权重 params.Rp 0.1 * eye(2); % 追方控制能量权重 params.Re 0.1 * eye(2); % 逃方控制能量权重 % 初始相对状态 x0 [5; 2; -1; 0.8]; % 位置5km,2km; 速度-1km/s,0.8km/s t0 0; tf 5; % 末端时间, 单位s这里还包含一个容易被忽视的陷阱相对运动方程中(A)矩阵的上三角块是单位矩阵对应速度和位置的关系但如果仿真之前在状态量顺序上不一致比如把共态量插到了中间某两行最后画出来的轨迹就会完全乱掉。4.3 打靶残差函数与循环迭代逻辑打靶法的核心代码在一个残差函数里。你需要给定初始共态猜测积分到终端时刻然后返回横截条件残差。我用的迭代方法是fsolve选项中的trust-region-dogleg算法但初值猜测质量不行时这个算法也经常翻车。后来我加了一个最外层的小技巧先把终端时间缩短一半求出解再把上一组结果插值后作为长终端时间的初值猜则相当于做了一个路径延拓稳定性大幅提升。function F shooting_residual(lambda0_guess, params) y0 [params.x0; lambda0_guess]; opts odeset(RelTol, 1e-10, AbsTol, 1e-12); [~, Y] ode45((t, y) dynamics_rhs(t, y, params), ... [params.t0, params.tf], y0, opts); yT Y(end, :); F1 yT(5:8) - params.Qf * yT(1:2); % 把终端位置残差也叠加进来加速收敛 F2 yT(1:2) - [0; 0]; F [F1; F2]; end加位置残差项是我自己试验出来的。如果不加初始猜测很差的时候迭代初期只靠共态横截条件驱动收敛很慢甚至震荡加了位置残差之后相当于给迭代过程加了一个“把终端拉向原点”的引导力能让求解过程更快进入物理上合理的区域。不过要注意位置残差会改变方程组形式严格来说这时已经不是原始的横截条件问题了所以最终判断收敛时应该以横截条件残差为主位置残差只作为辅助参考。5. 仿真结果分析怎么判断均衡真的被找到了5.1 状态轨迹与控制量剖面解读仿真跑完之后第一件事不是看脱靶量而是看双方的相对运动轨迹。我画的是以追方质心为原点的相对位置曲线横坐标为视线方向位置纵坐标为法向位置。正常情况下追方和逃方的轨迹在初始阶段会有一段明显的“对峙”过程双方都在试探性地调整方向随后在末端时刻快速收敛到一个临界区域。如果画出来的轨迹在中间出现异常的S形扭结或者振荡多半不是真实博弈行为而是打靶法收敛到了非物理解的征兆。真正的均衡策略剖面应该相对平滑因为二次型支付函数和控制能量惩罚决定了控制量不可能在短时间内剧烈跳变除非末端时间太短导致边界层效应。控制量剖面图更能说明问题。追方控制量(u_p)和逃方控制量(u_e)在时间轴上是镜像关系末端时刻双方的控制量一般都会放大因为双方需要在最后一刻尽可能实现各自目标。但放大过量就要注意了如果超过了你设定的过载上限(A_p^{max})和(A_e^{max})说明这个解不满足实际约束要么放宽支付函数里的控制权重要么缩短博弈终端时间要么把最大加速度约束直接作为路径约束加进求解过程。5.2 ε-纳什均衡验证偏离策略测试我在文章开头提到解出极值条件还不够必须做偏离策略测试。这个测试的做法是固定追方的最优策略不变允许逃方在一组预设的有界扰动策略集合里选择任意一个然后计算逃方实际得到的支付相比均衡支付的最大提升量。如果这个最大提升量不超过设定的ε值那么该策略组合就是ε-纳什均衡。反过来固定逃方策略、扰动追方策略验证追方的支付损失也不超过ε。实际仿真中我让每个策略的扰动系数在[-0.5, 0.5]区间里等间隔采样每个维度取21个点组合出来的扰动策略集有几千个。逐一遍历计算支付确实比较耗时但这一步的价值在于它能给出一个统计意义上的“均衡质量”报告而不仅仅是一根优美的轨迹曲线。典型的结果是在ε0.01时随机扰动产生的最大支付改善基本在0.003到0.008之间满足均衡条件随着扰动幅度增大到接近甚至超过均衡策略本身的量级偏离带来的支付变化才可能超过ε。这个结果跟理论预期一致——ε-纳什均衡在局部范围内是稳定的但并不意味着全局唯一。5.3 参数敏感性分析与鲁棒性观察参数敏感性分析是我临时加进去的一步因为光看一条仿真曲线很难说服自己这个方法有实用价值。我做了三组扫描改变初始相对速度、改变控制能量权重比值、改变终端时刻每组都重新求解均衡并统计脱靶量变化。结果发现这个模型在终端时刻从3秒变化到8秒时脱靶量整体趋势稳定只在某个中等时间值附近出现了一个局部极小值点控制权重比值从0.05提高到0.2脱靶量会缓慢上升符合“控制越贵、精度越差”的直觉。这里有个有意思的发现在末端追逃博弈里终端时刻对均衡的影响不是单调的。如果时间太短双方来不及做机动追方即使使用了均衡策略也无法在终端时刻收敛到位脱靶量自然很大如果时间太长逃方有更多机会利用机动耗尽追方的控制资源追方后期可能出现控制饱和脱靶量反而回升。这个非单调特性说明末端时间的选择本身就是博弈设计的一部分不应该拍脑袋定。6. 常见问题、调试经验与复现避坑指南6.1 初值猜测太差导致求解发散这是我在这类边值问题上遇到过最频繁的问题。打靶法本质上是一个非线性方程求根问题初值猜测直接决定了迭代是否会收敛到某个解。如果(\lambda(t_0))的初值猜测偏离真实解太远几乎必然发散。我的经验是“多管齐下”解决从简单情况开始先把终端时间设成很小的值比如1秒这样最优控制问题和它的边值问题都比较温和容易收敛。拿到解之后把共态曲线和时间网格延拓到更大时间范围。善用线性二次型近似先把问题当作不带终端约束的LQR来解用LQR得到的Riccati解作为共态的初始轮廓再喂给打靶迭代初值质量大幅改善。已经发散时不要反复改变fsolve的迭代次数而是回退初值猜测在同一个初值上多试几种算法比如Levenberg-Marquardt和trust-region换着用。6.2 数值结果对相对运动方程比例尺度的敏感我强烈建议做一个无量纲化处理。如果状态量直接使用国际单位制比如相对位置量级在千级别以km为单位而控制量量级在几到几十之间状态和共态的数值尺度差异能差出成百上千倍fsolve的Jacobian会变成近似病态矩阵迭代一步就飞掉。做无量纲化其实不难把位置、速度、时间分别除以参考值和终端时间得到的新状态量都在1附近波动求解稳定性立刻上去。6.3 控制量超出边界与支付函数权重调节的配合如果仿真结果里控制量剖面明显超出最大加速度限制光靠调小(R_p)和(R_e)不一定能解决问题因为缩小权重意味着控制量整体上移。我最终的做法是先把求解出来的控制量裁剪到边界内然后重新用裁剪后的策略做一次偏离统计分析看看裁剪造成的支付损失是否仍然保持在ε范围内。如果裁剪导致的均衡破坏不大这个策略依然有工程参考价值如果破坏很大就需要把控制约束以充分必要条件的形式加入最优性条件比如用Pontryagin极值原理中的控制约束变体来重推边界条件。7. 复现过程中的一些额外心得这个项目的完整技术链路从我搭建模型、推导最优性条件、迭代求解到最终验证均衡整个过程不需要太多“灵光一闪”的天才成分更多是工程上的耐心和数值上的细致。几个最值得记住的经验第一不要把任何一个数值结果当作“真实解”。边值问题的解在数学上严格成立但工程上是否合理还要依靠物理直觉来判断比如控制量剖面异常、终端位置残差恰好为零但轨迹明显偏离物理规律这些都是求解器在告诉你某处设置出了问题。第二验证比求解重要。在很多仿真报告里大家都乐意展示好看的控制轨迹和收敛曲线但“这组策略为什么是均衡策略”这个问题反而被一笔带过。亲手写一遍偏离策略测试之后你会对ε-纳什均衡这个概念的工程含义有完全不同的理解它不只是理论上锦上添花的界定更是判断仿真结果可信度的关键依据。第三无量纲化和参数延拓是边值问题求解的两个“救命稻草”尤其是后者。我做过的几乎每一次边值问题复现都会先用延拓法把参数从简单情况推演到复杂情况。如果直接跳进去挑战完整难度配置大概率会把大量时间耗在初值尝试和发散排查上这跟实际的项目进度管理规律完全一致——先把路探通再走大路。最后建议想做这个方向的朋友在拿到Matlab代码后不要急着跑图。先把仿真参数表整理出来明确每个参数对应的物理含义和单位再把动力学函数单独拿出来用一组开环控制输入验证状态更新逻辑是否正确最后才进入均衡求解环节。这就像写任何一段复杂程序之前先写单元测试一样能帮你把后面出问题的范围缩小一半以上。