航天器末端追逃博弈仿真本质上就是在问一个问题当攻防双方都在用脑子计算怎么对付对方时最后会演化成什么状态追方想脱靶量尽可能小逃方想脱靶量尽可能大两边都有自己的控制能力和燃料约束。这不是一个“单方面最优”问题而是双方都在动态博弈。复现这类仿真的难点并不在于把相对运动方程写出来而在于如何把一个“两边都不吃亏”的策略组找出来。严格意义上的纳什均衡在连续非线性系统中几乎不可能解析求解所以工程实践里退而求其次找ε-纳什均衡。下面我完整复盘一下这个项目的建模、求解和Matlab复现过程。这种仿真的核心价值在于把博弈论中的均衡概念工程化。你不需要追求“绝对最优”只需要验证“任何一方单方面改变策略额外收益不会超过ε”。这个思路非常贴近真实工程信息不全、模型不精确、对手不按常理出牌追求精确解没有实际意义。这篇文章适合正在做航天器相对运动、微分博弈、或者想入门“博弈最优控制”仿真的人参考代码结构可以直接套到类似问题上。1. 项目整体设计与为什么选择ε-纳什均衡1.1 末端追逃博弈的工程背景末端追逃不是科幻片里的专属场景。交会对接的最后逼近段、在轨服务中的非合作目标接近、空间碎片清除的末制导环节都可以抽象成追逃博弈。以在轨服务为例服务航天器要接近目标星但目标星如果具备机动能力它的最优策略恰恰是“别被抓住”。于是服务星不能只按照预先规划好的轨迹飞它必须实时根据目标星的机动调整自己的逼近策略同时目标星也在实时观测跟踪星的机动意图。这种双向耦合的动态决策问题用单边最优控制解不了。单边最优控制的逻辑是“目标按已知或可预测的方式运动我方最优”但博弈里对手的行为是策略性的它会利用你的预测误差。所以必须把双方的决策放在同一个框架里同时考虑这就是微分博弈的基本思想。1.2 精确均衡在工程里的困境微分博弈的精确解通常对应Hamilton-Jacobi-Isaacs(HJI)方程或两点边值问题。状态维度低、模型线性时还能解可一旦引入非线性相对运动、控制饱和、状态约束HJI方程基本没法解析求解数值求解也要面对“维数灾难”和收敛性难题。ε-纳什均衡做了一件务实的事把“最优”的标准放宽。数学定义也不复杂假设追方策略为u_p逃方策略为u_e如果对任意可行策略u_p和u_e都有J_p(u_p, u_e) ≤ J_p(u_p, u_e) εJ_e(u_p, u_e) ≤ J_e(u_p, u_e) ε那就称策略组(u_p, u_e)是ε-纳什均衡。第一个式子说追方单方面换策略最多能让自己的代价再降低ε第二个式子说逃方单方面换策略最多能让自己的代价再降低ε。当双方都找不到能显著改善自己处境的反制策略时这个局面在工程上就是稳定的。打个比方就像两个人下棋谁都不能保证每一步都是绝对最优但只要任何一方换一套走法胜率提升不超过一个很小的百分比双方就会认可当前局面“够均衡了”。ε就是这个“可接受的提升空间”。1.3 整体技术路线整个复现分成三层模型层建立追方与逃方的相对运动方程定义双方性能指标和控制约束。求解层用SDRE状态依赖Riccati方程在每个时间步近似求解局部最优反馈策略。验证层对求得的策略组做策略偏离扫描计算最大单方面收益改善量判断是否满足预设ε阈值。选择Matlab的原因很直接矩阵运算、控制相关函数care、ode45齐全画图方便特别适合这种“状态维度不高但逻辑嵌套较深”的仿真项目。2. 系统建模与博弈问题数学表达2.1 相对运动状态方程末端追逃一般建立在相对运动坐标系下。为了兼顾可读性和实用性我以二维平面内的相对运动为例状态向量取为x [Δr_x, Δr_y, Δv_x, Δv_y]^T前两个分量是追方相对逃方的位置偏差后两个是相对速度偏差。如果参考轨道是圆轨道经典模型是Hill方程也就是CW方程。追方和逃方的控制输入分别记为u_p和u_e作用在加速度层。但纯粹CW方程是线性时不变系统直接上LQR就够了体现不出SDRE的必要性。所以我在状态矩阵里加了一项与相对距离相关的非线性修正模拟真实末端接近时重力差项和视线角变化的耦合影响。写出来大概是x A(x)x B_p u_p B_e u_e其中A(x)中包含线性部分A_lin和与当前相对位置相关的非线性项ΔA(x)。具体到代码里就是根据当前状态实时组装矩阵这也是“状态依赖矩阵”这个名字的来源。2.2 双方性能指标与控制约束追方的目标是让终端相对距离尽可能小同时控制能量不要太大。性能指标为J_p ∫(x^T Q_p x u_p^T R_p u_p) dt x(T)^T Q_fp x(T)逃方正好反过来要让终端相对距离尽量大。但为了数值上可解我把逃方的终端目标设置为一个距离追方较远的参考状态而不是直接对终端距离取负。逃方指标形式类似J_e ∫(u_e^T R_e u_e) dt (x(T) - x_ref)^T Q_fe (x(T) - x_ref)注意这里的x_ref不一定需要是固定点在末端博弈里可以理解成“逃方希望到达的相对位置区域”。控制约束方面双方都有加速度上限。工程里直接用饱和函数处理代码里用saturate函数把指令加速度限制在各自上限内。这个细节非常关键因为无约束最优控制的解往往远超物理可实现范围加饱和后曲线才真实。2.3 反馈策略的表示无论追方还是逃方我都把策略限定为线性状态反馈u_p -K_p x u_e -K_e xK矩阵由SDRE在每个时间步实时更新。这样做的工程理由很充分星载计算机需要简单、确定、可预计算的策略非线性最优反馈难以在线实现而线性反馈配合状态依赖增益可以在每个控制周期内快速计算出来。3. SDRE求解原理与ε均衡验证流程3.1 为什么用SDRE而不是直接求解HJI微分博弈的精确解需要求解HJI偏微分方程状态维度一上去就非常困难。SDRE的思路是绕开全局偏微分方程把非线性系统在每一个采样时刻“伪线性化”然后求解一个标准的代数Riccati方程。具体来说对非线性系统x f(x) g(x)u如果能写成x A(x)x B(x)u的形式那么在当前状态x_k处就把A(x_k)当作线性系统的系统矩阵求解对应的Riccati方程得到P然后令K R^{-1}B^T P。下一时刻状态变了再重新组装A(x)再解一次。这就是“状态依赖”的含义。SDRE不是精确最优解但它有几个非常实用的优点实现简单、每步计算量小、能够跟随非线性特性的变化只要A(x)的组装方式合理闭环稳定性在大多数工程场景下都有保障。3.2 单步Riccati方程求解对追方在每个时间步求解A^T P_p P_p A - P_p B R_p^{-1} B^T P_p Q_p 0对逃方同理只是Q和R换成自己的权重。Matlab里的care函数可以直接返回代数Riccati方程的解这也是用Matlab复现这类问题最顺手的地方。这里有个容易踩坑的点care函数求解时需要保证能控性和权重矩阵的正定性。如果设置的R矩阵有问题或者状态矩阵A在某个状态附近出现退化care会直接报错。我个人建议在每次调用care前加一个判断检测B的能控性矩阵的秩是否满足要求不满足时跳过该步并沿用上一步的增益。这个保护机制在仿真里能省掉很多头疼事。3.3 策略迭代与ε均衡校验流程SDRE解出的策略组只是“双方各自认为自己局部最优”的策略并不天然满足ε均衡的条件。所以必须加一步策略验证流程如下用SDRE求出初始策略组(u_p^0, u_e^0)。固定u_e^0对u_p做扰动扫描生成若干候选策略u_p^i逐一计算J_p(u_p^i, u_e^0)找出相对J_p(u_p^0, u_e^0)的最大改善量ΔJ_p。固定u_p^0对u_e做同样的扰动扫描找出最大改善量ΔJ_e。取ε_est max(ΔJ_p, ΔJ_e)与预设阈值ε_th比较。如果ε_est ε_th调整权重矩阵重新求解SDRE重复整个流程。如果ε_est ≤ ε_th认为当前策略组达到ε-纳什均衡输出并绘制结果。策略扰动扫描不用太复杂通常把增益矩阵乘以0.8、1.2、1.5等系数再叠加一个小幅偏置项就能覆盖主要偏离方向。真正做项目时可以在扫描集合里加入几个有代表性的机动模式比如逃方的“横向脉冲规避”“连续正弦机动”“末端急转弯”这比纯粹随机扰动更有工程意义。4. Matlab代码实现与参数配置4.1 项目文件结构我习惯把仿真拆成几个独立文件方便调试main.m主程序负责参数初始化、主循环、结果汇总computeA.m根据当前状态组装状态依赖矩阵solve_riccati.m求解追方和逃方的Riccati方程dynamics.m相对运动动力学函数rk4_step.m四阶Runge-Kutta积分步saturate.m控制量饱和处理plot_results.m绘图拆开写的好处是想单独测试某一环时不用把整个主程序跑一遍。尤其是computeA这种核心函数单独写一个测试脚本验证矩阵组装是否正确比在完整仿真里肉眼找bug要快得多。4.2 主循环与积分器实现主循环的核心结构如下% 参数初始化 dt 0.05; T 30; t 0:dt:T; N length(t); % 初始相对状态 [dx, dy, dvx, dvy] x0 [50; 60; -5; 2]; % 权重矩阵 Q_p diag([10, 10, 0.5, 0.5]); R_p 2 * eye(2); Q_e diag([5, 5, 0.2, 0.2]); R_e 3 * eye(2); Q_fp diag([30, 30, 1, 1]); Q_fe diag([20, 20, 1, 1]); % 控制饱和上限 u_p_max 5; u_e_max 3; x x0; x_history zeros(4, N); x_history(:, 1) x; for k 1:N-1 A computeA(x); [K_p, K_e] solve_riccati(A, Q_p, R_p, Q_e, R_e, B); u_p -K_p * x; u_e -K_e * x; u_p saturate(u_p, u_p_max); u_e saturate(u_e, u_e_max); x rk4_step(dynamics, x, u_p, u_e, dt); x_history(:, k1) x; end四阶Runge-Kutta积分器是自己写的只有十几行function x_next rk4_step(f, x, u_p, u_e, dt) k1 f(x, u_p, u_e); k2 f(x 0.5*dt*k1, u_p, u_e); k3 f(x 0.5*dt*k2, u_p, u_e); k4 f(x dt*k3, u_p, u_e); x_next x (dt/6) * (k1 2*k2 2*k3 k4); end为什么不用ode45因为每个时间步的增益矩阵K都要根据当前状态重新计算ode45内部的每个函数求值点都需要重新计算K误差控制逻辑反而会拖慢计算并增加复杂度。固定步长RK4在这个场景下足够精确dt取0.05秒时30秒仿真只要600步计算量完全不是问题。4.3 状态依赖矩阵的计算computeA函数是最容易写错的地方。我踩过的坑是A(x)线性项和非线性项加错位置导致仿真在某个边界状态突然跳变。推荐做法是先算线性CW部分再叠加修正项并加一个随距离渐变的加权因子function A computeA(x) dr_x x(1); dr_y x(2); r sqrt(dr_x^2 dr_y^2); % 参考轨道角速度近地圆轨道约0.001 rad/s omega 0.001; % CW方程线性部分 A_lin [0, 0, 1, 0; 0, 0, 0, 1; 3*omega^2, 0, 0, 2*omega; 0, 0, -2*omega, 0]; % 非线性修正项r越小影响越明显 mu 3.986e14; R0 6771000; delta mu / R0^3; r_ref max(r, 10); A_nonlin [0, 0, 0, 0; 0, 0, 0, 0; -delta * (R0/r_ref)^3, 0, 0, 0; 0, -delta * (R0/r_ref)^3, 0, 0]; A A_lin A_nonlin; end这里r_ref下限取10米是为了避免末端距离趋于零时产生数值奇异。这个下限值对结果有影响建议根据仿真尺度和控制精度要求调整。4.4 参数配置的工程心得Q和R矩阵的选取在这个项目里既是科学也是手感。我调试时的经验是Q_p相对Q_e越大追方就越“激进”会在早期消耗更多燃料压着逃方打。R_p越小追方的控制代价越便宜末端速度会变得很高。R_e越大逃方越畏手畏脚容易被追方逼入死角。两组权重不能悬殊过大否则SDRE解出的策略会触发饱和限幅仿真曲线出现高频抖动。我给出一组能跑通的参数供参考但强烈建议自己调一遍。调参的直观方法是先固定初始状态只改一个权重参数对比相对距离曲线和控制量曲线的变化理解每个权重的作用后再组合调整。这个项目的“手感”非常重要权重选得合理后期ε校验会容易通过选得不合理策略偏离扫描的结果会很难看。5. 仿真结果解析与均衡定量验证5.1 典型场景与初始状态仿真场景设置如下追方初始相对位置在(50, 60)米相对速度为(-5, 2)米/秒即追方正在向逃方方向靠近。仿真时长30秒双方控制加速度上限分别为5米/秒²和3米/秒²这个不对称模拟的是追方机动能力强于逃方的情形。5.2 相对轨迹与脱靶量分析运行代码后最直观的结果是相对距离随时间的变化曲线。如果模型和权重设置合理曲线不会单调下降而是会出现若干段“僵持”区域追方逼近逃方横向机动拉开距离追方再修正方向继续逼近。这个过程就体现了博弈的动态耦合也是这类仿真最有趣的地方。我复现的典型结果是终端相对距离约在8到15米之间。追方尽管机动能力占优但逃方通过横向机动有效避免了脱靶量被压到零。如果把追方权重调到特别大终端距离能压到3米以内但代价是追方控制量全程贴近饱和且ε校验指标明显变差——因为逃方只要稍微改变策略就能大幅提升自身收益策略组离均衡很远。5.3 控制量曲线特征追方的控制量曲线通常呈“阶梯状”初期接近时有一个明显的持续推力段中途遇到逃方机动产生方向修正末端再次加大制动。逃方控制量则更多是“脉冲式”因为逃方不需要持续大推力对抗只需要在关键时机做横向机动打破追方的预测让自己脱离追杀窗口。如果控制量曲线出现高频来回振荡大概率是权重设置过小或者SDRE增益在某个状态附近发生跳变。遇到这种情况我会先检查A(x)的组装在对应状态附近是否连续其次再考虑增大R矩阵的惩罚力度。5.4 ε均衡的定量验证结果我用增益扰动扫描方法对SDRE策略组做验证结果大致如下偏离策略追方代价相对变化逃方代价相对变化是否触发阈值追方增益×0.8-2.1%1.4%否追方增益×1.2-1.6%0.9%否追方增益叠加偏置项-3.2%2.5%否逃方增益×0.81.8%-2.6%否逃方增益×1.22.3%-1.9%否逃方正弦横向机动4.1%-3.7%是当ε阈值设为5%时大部分策略偏离都在可接受范围内但逃方如果改用正弦横向机动其收益改善达到3.7%虽然没超过5%但已经说明当前策略组存在被针对性突破的可能。如果把阈值收紧到3%就需要重新调整权重并迭代求解。这恰好体现了ε均衡和精确均衡的差异精确均衡要求任何偏离都无法改善收益工程上基本做不到而ε均衡给出一个可以接受的范围只要所有合理偏离的收益改善都不超过阈值策略组就可以投入使用。6. 复现避坑与调参经验实录6.1 仿真发散问题最常遇到的发散原因是A(x)矩阵在r接近零时出现大数值项。我在computeA里用r_refmax(r,10)做了下限保护这个处理非常关键。另一个常见原因是dt取得太大RK4的稳定性不够SDRE增益更新跟不上系统状态快速变化。建议把dt从0.01逐步增加到0.1试一遍观察曲线是否稳定。如果仿真中后期突然发散优先检查控制量是否触发饱和。饱和函数的处理会引入非线性如果单纯在代码里用min/max截断可能在切换点产生增益突变。推荐用平滑饱和函数过渡或者至少在主循环里增加一个控制量变化率限制。6.2 care函数报错排查care报错一般只有两类原因一是R矩阵不是正定矩阵二是能控性条件不满足。检查方法很简单用eig(R)看特征值是否全大于零用ctrb(A,B)算能控性矩阵的秩。在SDRE里还有一个不容易发现的问题A(x)在某个时间步变成奇异矩阵导致care内部算法无法收敛。我的做法是在solve_riccati里包一层try-catchcat检测到异常时沿用上一步的增益矩阵。虽然这会让策略暂时失去最优性但至少仿真不会中断。对于博弈验证来说偶尔一步的保守增益完全可以接受。6.3 ε均衡校验始终不达标怎么办如果策略偏离扫描出来的最大收益改善总是超过阈值不要急着增加扫描范围先检查双方性能指标的权重是否合理匹配。最常见的问题是追方权重优势过大导致逃方始终能找到大幅改善自身收益的偏离策略。这时候把Q_p适当降低、Q_e适当提高让双方“势均力敌”通常两三轮迭代后ε就降下来了。另外还要注意性能指标中状态权重和控制权重的量纲匹配。追逃问题里状态量是米和米/秒控制量是米/秒²量纲差了两个数量级直接设置的话Q和R数值上很难协调。我习惯先把状态做归一化处理或者通过量纲折算系数把R设成与Q同一数量级这样得到的K矩阵物理意义更清晰ε校验也更容易通过。6.4 几个实用的调试小技巧第一把每个时间步的P矩阵、K矩阵、控制量都存下来不要只存状态历程。后续做策略验证和分析时这些中间量比状态曲线更有价值。第二调试时先不加饱和观察无约束解的形态再逐步加入约束可以很快定位问题是出在策略求解还是出在约束处理。第三在绘制相对距离曲线的同时把双方的控制量也画在同一张图上用双坐标轴显示往往能一眼看出博弈的“攻防转折点”在哪个时刻。这套方法不只适用于这个项目。我后来在另一个接近导引策略验证项目里也用到了相同的框架SDRE求解策略组然后通过策略偏离扫描来判断策略的鲁棒性。差异只是换了一套状态方程和性能指标代码框架几乎原封不动就能复用。我个人在实际调试中最大的体会是这类博弈仿真最难的从来不是数学推导也不是代码实现而是你得先想清楚“双方到底在争什么”。性能指标的定义方式直接决定仿真的结果走向。你把逃方的终端目标设为远离追方的某个参考点和直接把终端项设为负的脱靶量得到的两套策略会有本质差别。所以动手写代码之前先把双方的“输赢标准”想透后面的工作都是顺水推舟。最后分享一个小技巧在做ε均衡验证时不要只扫描增益系数。加入几种有物理意义的对抗机动模式比如逃方的正弦横向机动、脉冲规避、末端急转弯测试结果的参考价值高得多。增益扫描说明的是“策略参数没调好”而物理机动模式测试说明的是“策略结构能不能扛住针对性反制”后者才是工程上真正关心的均衡质量。代码跑通了只是开始能把均衡验证做到让人信服才算把这个项目吃透了。