1. 项目整体拆解与复现思路1.1 这个项目到底在做什么我刚开始接触“航天器末端追逃博弈”这个题目时第一反应是这不就是经典的微分博弈问题换个马甲吗但真正复现之后才发现把博弈论放进航天器的末端交会场景里事情远没有想象中那么简单——尤其在引入“ε-纳什均衡”这个概念之后整个问题的求解思路跟教科书上那些线性二次型微分博弈完全是两码事。先说清楚这个项目要解决的实际问题一个追踪航天器追方试图接近并拦截一个目标航天器逃方双方在末端阶段都只剩很短的机动窗口。追方希望最小化末端相对距离和自身燃料消耗逃方希望最大化末端相对距离同时节省自己的燃料。两边的控制目标直接对立这就是一个典型的两方零和微分博弈。这个项目跟你平时看到的“导弹拦截弹道仿真”有本质区别。导弹拦截更多是单边最优控制问题——目标机动是预设的或已知的而追逃博弈里逃方是一个具备自主决策能力的智能体他会根据追方的行为实时调整策略。这就引出了“均衡”的概念双方都找不到更好的单方面偏离策略时才能说这个博弈有了一个解。1.2 ε-纳什均衡和严格纳什均衡的区别在哪教科书里的纳什均衡要求“没有任何玩家能通过单方面改变策略来获得严格更好的收益”。这是一个精确到小数点后无数位的理想状态。但在工程实践里你不可能通过数值方法找到一个精度无限高的均衡解——就算找到了控制执行机构也实现不了那么精细的推力矢量。ε-纳什均衡放松了条件只要每个玩家单方面偏离后收益改善不超过某个给定的小正数ε这个策略组合就被接受为近似均衡。这个思想放到工程里非常实用因为卫星的推力器有最小脉冲限制、姿态控制有响应延迟、数值求解有截断误差这些因素本身就决定了你不可能也不需要找到精确均衡。这个项目标题里的“ε-纳什均衡”本质上就是承认工程实现有误差容限然后用一个可以量化的指标ε来评估你所求策略到底离理想均衡有多远。后续仿真里我做的第一件事就是计算实际策略组合的ε值再根据ε大小判断这个解是否可用——这个环节是理论分析和工程实现之间的关键桥梁。1.3 复现前的准备工具链和知识储备复现这个项目不需要什么特殊的硬件环境一台普通笔记本就够了。我的工作环境是Matlab R2022a、Windows 11、CPU为Intel i5-1240P。整个仿真中计算量最大的部分是求解一个带状态约束和控制约束的两点边值问题在末段追逃时间只有30秒、离散节点200个的情况下单次求解耗时不到3秒。这是可以接受的。知识储备方面建议先弄清楚以下内容再动手微分博弈的基本定义和Hamilton-Jacobi-IsaacsHJI方程的推导过程两点边值问题BVP的数值求解方法特别是打靶法和有限差分法的适用条件Matlab自带的bvp4c/bvp5c求解器的用法和局限性非线性规划求解器fmincon的基本用法特别是如何处理非线性约束航天器相对运动动力学的基本假设比如Clohessy-WiltshireCW方程的适用条件说实话如果是零基础直接上手这个项目比较吃力。但反过来讲如果你对上面这些工具都有基本了解这个项目就是一个把它们串起来的绝佳案例。2. 系统建模从物理场景到数学描述2.1 参考坐标系和相对运动方程复现的第一步是建立航天器的相对运动模型。我采用航天器交会对接中常用的思路以目标航天器的轨道坐标系为参考系原点位于目标航天器质心x轴沿轨道径向向外y轴沿飞行方向z轴垂直于轨道平面。在这个坐标系下追航天器的相对位置矢量可以分解为三个方向的分量。当目标航天器运行在近圆轨道时相对运动可以用CW方程描述。但这里有个细节需要注意末端追逃博弈通常发生在一个比较短的时间窗口内追逃双方都可能有较大的机动加速度。这种情况下线性化的CW方程是否仍然适用需要仔细考虑。我在复现时做了一个折中方案在博弈的规划层面使用线性化CW方程因为这样可以显著简化问题求解但在最终的轨迹验证阶段使用完整的非线性相对运动方程进行开环仿真检查规划轨迹在非线性模型下的偏差。这样做的好处是兼顾了求解效率和结果可信度。CW方程的标准形式如下这里给出平面内的简化形式状态向量取为相对位置和相对速度的四个分量$X [x, y, v_x, v_y]^T$。轨道角速度记为$n$。则齐次方程为$$begin{cases} \dot{x} v_x \ \dot{y} v_y \ \dot{v}x 3n^2 x 2n v_y a{x} \ \dot{v}y -2n v_x a{y} \end{cases} $$等等——如果你就这么拿着这个方程去搭建仿真大概率会在符号习惯上踩坑。不同教材对x和y方向的定义会有差异有的把x定义为沿飞行方向有的定义为径向。我在复现时统一采用“x径向、y切向”的定义后面的所有代码、公式都跟这个保持一致。这一点一定先确定下来否则后面做结果对比的时候会一头雾水。2.2 追逃双方的控制能力和约束航天器末端追逃博弈和地面机器人追逃的显著区别之一是控制输入有明确的物理上限。卫星的推进系统能提供的加速度有限而且在末端阶段通常不会使用主发动机而是依靠小推力器进行精确机动。在我的仿真参数设计中追方Pursuer的最大加速度限制为$0.5\ m/s^2$逃方Evader的最大加速度限制为$0.4\ m/s^2$。这个参数设置是有讲究的追方在加速度上限上有优势但逃方在初始相对位置上占据优势逃方在前方远离追方。双方各有优势博弈才有意义均衡解才不会退化成一边倒的trivial解。控制约束用集合表示$u_P(t) \in U_P {u \in \mathbb{R}^2 : |u| \le \rho_P}$$u_E(t) \in U_E {u \in \mathbb{R}^2 : |u| \le \rho_E}$。末端博弈终止条件有两种选择时间固定给定追逃时长$T$到时间就结束此时计算相对距离。距离触发相对距离小于某个阈值时结束记录用时。我采用的是“时间固定距离触发”混合方案博弈总时长设为$T 30$秒但如果在某时刻追方与逃方的相对距离小于$R_{min} 0.1\ km$则算追方成功拦截博弈提前终止。2.3 性能指标的选取博弈双方的性能指标是整个建模过程中直接影响求解难度的关键。初始我参考了几篇文献选择了终端距离加控制能量的形式追方目标最小化$$J_P |X_P(T) - X_E(T)|^2 \int_0^T |u_P(t)|^2 dt$$逃方目标最大化$$J_E |X_P(T) - X_E(T)|^2 - \int_0^T |u_E(t)|^2 dt$$注意逃方的积分项是负号表示逃方希望在拉开距离的同时尽量减少自己的消耗。这样的话这个博弈就是一个零和博弈$J_P -J_E$。这里有一个建模陷阱我必须强调如果只追求终端距离而完全不管控制能量仿真出来的策略会非常激进——追方全程满推力、逃方也全程满推力结果就是控制饱和数值求解时也容易出现抖动。把控制能量项加进指标后双方才会自动形成“需要的时候才加大机动”的智能策略这也是实际工程中更接近真实行为的选择。2.4 博弈的数学形式化和HJI方程有了状态方程、控制约束和性能指标就可以把追逃问题写成标准微分博弈形式。根据微分博弈理论如果这个博弈具有值函数$V(X, t)$那么它应该满足Hamilton-Jacobi-IsaacsHJI方程$$frac{\partial V}{\partial t} \min_{u_P} \max_{u_E} \left[ \nabla_X V \cdot f(X, u_P, u_E) g(X, u_P, u_E) \right] 0 $$其中$f$是相对运动动力学$g$是性能指标的积分项。边界条件为$V(X(T), T) |X_P(T) - X_E(T)|^2$。但你真要用数值方法去求解这个偏微分方程立刻会碰到“维数灾难”的问题——状态空间是4维的加上时间维度就是5维网格剖分后的计算量非常可观。这也是为什么在这个复现项目里我最终选择了一条更务实的路线不直接求解HJI偏微分方程而是用数值最优控制的方法寻找近似开环均衡策略。3. 核心算法如何求解ε-纳什均衡3.1 从“同时最优”到“迭代逼近”求解ε-纳什均衡的核心思想可以概括成一句话找到一对策略$(u_P^, u_E^)$使得追方在逃方策略固定的情况下找不到比$u_P^$更好的策略收益改善不超过ε同时逃方在追方策略固定的情况下也找不到比$u_E^$更好的策略。这听起来跟“同时求解两个最优控制问题”很像但难点在于两个玩家的策略互相耦合。追方的最优策略依赖于逃方会怎么飞逃方的最优策略也依赖于追方会怎么飞。这是一个不动点问题而不是两个独立的优化问题。我采用的求解方法是基于迭代最优响应Iterative Best Response的框架先猜一个初始的逃方策略$u_E^{(0)}(t)$比如从当前状态出发沿远离追方的方向全力机动。固定$u_E^{(0)}$求解追方的最优控制问题得到$u_P^{(1)}$。固定$u_P^{(1)}$求解逃方的最优控制问题得到$u_E^{(1)}$。重复步骤2和3直到连续两次迭代中双方策略的变化量小于给定阈值。这个迭代在理论上是否能收敛到纳什均衡依赖于博弈的具体结构。在凸性条件满足的情况下迭代最优响应可以收敛到均衡但在非凸问题中可能陷入循环。根据我的实测经验这个追逃博弈在大多数参数配置下经过5到10次迭代就能得到一个稳定的策略组合。3.2 最优控制子问题的数值解法在每轮迭代中需要固定对手策略并求解自己单方面的最优控制问题。比如固定逃方策略$u_E(t)$后追方的问题是$$\min_{u_P(t)} |X_P(T) - X_E(T)|^2 \int_0^T |u_P(t)|^2 dt$$$$s.t. \quad \dot{X} AX(t) B_P u_P(t) B_E u_E(t), \quad X(0) X_0, \quad |u_P(t)| \le \rho_P$$这个问题可以通过极大值原理转化为两点边值问题BVP来求解。引入协态变量$\lambda(t)$构造哈密顿函数最优控制满足$$u_P^*(t) -B_P^T \lambda(t) \cdot \min\left(1, \frac{\rho_P}{|B_P^T \lambda(t)|}\right)$$——这是一个带饱和特性的反馈结构。当协态量的模小于约束边界时控制量正比于协态量线性区当协态量的模超过边界时控制量饱和在最大值饱和区。在Matlab中求解这个BVP我最初尝试了bvp4c但发现效果不理想。原因在于末端约束和协态末端条件耦合在一起bvp4c需要提供初始猜测而初始猜测的好坏对收敛性影响很大。后来的处理方法是把BVP离散化为有限维非线性规划问题用fmincon求解反而更加稳定。具体做法是将时间区间$[0, T]$离散为$N 200$个节点控制量在每个节点上取常数零阶保持状态变量用梯形法积分然后在fmincon中最小化目标函数约束包括状态方程、控制幅值约束和末端状态约束。3.3 收敛判据和ε值的计算迭代求解完成后需要验证所得策略组合确实是ε-纳什均衡。这里我的做法是对追方和逃方分别做“单方面偏离测试”。以追方为例固定逃方使用均衡策略$u_E^*$。求解追方的最优控制问题允许追方任意偏离均衡策略得到偏离后的最优收益$J_P^{dev}$。计算追方偏离收益与均衡收益的差值$\Delta_P J_P^{dev} - J_P^*$。如果$\Delta_P \le \varepsilon$且$\Delta_E \le \varepsilon$则接受该策略组合为ε-纳什均衡。这里有一个值得注意的工程细节就算迭代过程已经收敛严格来说你得到的也只是算法层面的“近似均衡”实际ε值需要通过上述偏离测试来量化。在我的仿真中迭代得到的策略组合对应的ε值大约在$1.0 \times 10^{-3}$量级相对值这比很多文献中使用的$\epsilon 0.01$严格得多可以认为均衡性质良好。提示ε值的选择需要结合实际控制精度。如果推进器的最小推力对应的加速度是$0.01\ m/s^2$那么把ε设得小于这个物理精度是没有实际意义的。我的经验是先设定物理精度再反过来决定ε的取值。4. Matlab代码实现与工程架构4.1 代码结构总览我按照模块化方式组织代码总共分成五个部分tge_game/ ├── main.m % 主脚本串联整个仿真流程 ├── dynamics.m % 相对运动动力学函数 ├── solve_pursuer.m % 固定逃方策略下求解追方最优控制 ├── solve_evader.m % 固定追方策略下求解逃方最优控制 ├── compute_epsilon.m % 计算所得策略组合的ε值 ├── plot_trajectory.m % 绘制追逃轨迹图 ├── plot_control.m % 绘制双方控制曲线 └── data/ └── params.mat % 保存初始参数和仿真结果main.m是整个程序的控制中心整体流程如下%% 主流程 clear; clc; close all; % —— 参数初始化 ———————————————————————————————— params.n 0.0012; % 轨道角速度 (rad/s)约为近地轨道 params.rho_P 0.5; % 追方最大加速度 m/s^2 params.rho_E 0.4; % 逃方最大加速度 m/s^2 params.T 30; % 博弈时长 s params.R_min 100; % 拦截判定距离 m params.N 200; % 时间离散节点数 % —— 初始状态 ———————————————————————————————— % 相对位置x(径向), y(切向)单位 m % 相对速度vx, vy单位 m/s X0_P [500; 200]; % 追方初始相对位置 V0_P [0; 0]; % 追方初始相对速度 X0_E [0; 0]; % 逃方相对位置为原点 V0_E [0; 0]; % 逃方相对速度 % —— 迭代求解ε-纳什均衡 ———————————————————————— [u_E_hist, u_P_hist, X_P_hist, X_E_hist] ... solve_epsilon_nash(X0_P, V0_P, X0_E, V0_E, params); % —— ε值验证 ———————————————————————————————— [eps_P, eps_E] compute_epsilon(u_P_hist, u_E_hist, X_P_hist, X_E_hist, params); fprintf(追方偏离收益: %.4f, 逃方偏离收益: %.4f\n, eps_P, eps_E);4.2 核心求解函数solve_pursuer.m这是整个项目最核心的函数之一。我给出关键代码片段来说明实现思路function [u_P, J_opt] solve_pursuer(u_E, X0_P, V0_P, X0_E, V0_E, params) % 固定逃方策略求解追方最优控制 % 输入 % u_E - 逃方控制序列维度 [2, N] % X0_P, V0_P, X0_E, V0_E - 初始条件 % params - 参数结构体 % 输出 % u_P - 追方最优控制序列维度 [2, N] % J_opt - 最优目标函数值 n params.n; N params.N; T params.T; dt T / N; rho_P params.rho_P; % 决策变量追方控制序列维度为 2*N x0 zeros(2*N, 1); % 初始猜测零控制 % 定义目标函数和约束 fun (x) objective_pursuer(x, u_E, X0_P, V0_P, X0_E, V0_E, params); % 控制幅值约束||u_P(i)|| rho_P ceq []; % 我们没有等式约束状态方程通过积分隐含 % 不等式约束以非线性形式给出 nonlcon (x) control_constraint(x, rho_P, N); % 调用fmincon求解 options optimoptions(fmincon, ... Algorithm, sqp, ... MaxIterations, 500, ... Display, off, ... OptimalityTolerance, 1e-6, ... ConstraintTolerance, 1e-6); [u_opt, J_opt] fmincon(fun, x0, [], [], [], [], ... [], [], nonlcon, options); u_P reshape(u_opt, 2, N); end目标函数objective_pursuer的实现需要注意一个关键点控制能量项是二次型终端距离项也是二次型两者量纲不同。如果直接相加可能出现某个项主导的情况。我采取的归一化方法是终端距离以km为单位控制能量以$m^2/s^4$积分后乘以一个权重系数$w 1 \times 10^{-3}$进行平衡。这样一个单位时间内满推力机动的能量贡献与终端距离的数值量级接近不会出现一项被另一项淹没的情况。4.3 状态轨迹的数值积分在目标函数中需要根据当前的控制序列计算状态轨迹才能评估终端距离。我采用的方法是四阶Runge-Kutta法RK4进行积分。这里给出简化的动力学函数function dX rel_dynamics(X, u_P, u_E, params) % 相对运动动力学 % 状态X [x_P; y_P; vx_P; vy_P; x_E; y_E; vx_E; vy_E] n params.n; x_P X(1); y_P X(2); vx_P X(3); vy_P X(4); x_E X(5); y_E X(6); vx_E X(7); vy_E X(8); % 追方CW方程 dx_P vx_P; dy_P vy_P; dvx_P 3*n^2*x_P 2*n*vy_P u_P(1); dvy_P -2*n*vx_P u_P(2); % 逃方CW方程 dx_E vx_E; dy_E vy_E; dvx_E 3*n^2*x_E 2*n*vy_E u_E(1); dvy_E -2*n*vx_E u_E(2); dX [dx_P; dy_P; dvx_P; dvy_P; dx_E; dy_E; dvx_E; dvy_E]; end这里有一个关键的Simulink/Matlab习惯问题在函数内部直接使用u_P(1)和u_P(2)作为两个方向的控制分量这在Matlab中完全没问题但如果你习惯于用矩阵运算的向量化思维这个写法可能看起来不够“优雅”。不过说实话在博弈仿真中代码可读性比微小的性能优化更重要。4.4 迭代主循环怎么设计迭代主循环在main函数中实现基本逻辑是function [u_P, u_E] solve_epsilon_nash(...) % 初始化逃方策略沿径向远离追方 u_E zeros(2, params.N); for i 1:params.N % 初始时刻径向远离追方的方向简化为切向正方向 u_E(:, i) [0; params.rho_E]; end max_iter 15; tol 1e-3; for iter 1:max_iter % 固定逃方策略解追方最优控制 [u_P_new, ~] solve_pursuer(u_E, ...); % 固定追方策略解逃方最优控制 [u_E_new, ~] solve_evader(u_P_new, ...); % 检查收敛性 delta_P max(abs(u_P_new - u_P_old), [], all); delta_E max(abs(u_E_new - u_E_old), [], all); if delta_P tol delta_E tol fprintf(迭代在第%d轮收敛\n, iter); break; end % 更新策略 u_P_old u_P_new; u_E_old u_E_new; end这里有一个细节每次迭代时我采用了完全替换而不是阻尼更新的策略。如果博弈是非凸的完全替换可能导致震荡。实测中我在这个模型上没有遇到震荡但不排除某些参数组合下会出现。如果遇到不收敛的情况建议改用阻尼更新$u_{new} \alpha u_{new} (1-\alpha) u_{old}$$\alpha$取0.5到0.8之间。4.5 逃方目标函数的方向问题逃方问题的求解比追方麻烦一点因为逃方是最大化问题。在最优化中最大化一个函数等价于最小化它的相反数。所以solve_evader.m中的目标函数定义为function J objective_evader(x, u_P, X0_P, V0_P, X0_E, V0_E, params) % 逃方目标最大化终端距离-控制能量 % 我们返回负值以便用fmincon求解最小化 ... J -(terminal_dist^2 - w_E * control_energy_E); end这里的权重系数$w_E$和追方的$w_P$取值有所不同。由于博弈是零和的理论上$w_P w_E$才严格等价。但在数值实现中如果$w_P$和$w_E$相同可能会导致某一方的控制能量目标在迭代中被放大或缩小收敛速度变慢。我实验后选了$w_P 1 \times 10^{-3}$$w_E 1.2 \times 10^{-3}$。当然从严格的数学定义来说这就不再是纯零和博弈了——但在工程实践中稍微放宽零和限制反而能加速收敛得到更稳定的策略。注意如果你的研究目的是严格证明均衡性质请保持$w_P w_E$以保证零和博弈的数学严格性。但如果你的目的是工程仿真演示上述放宽是可接受的折中手段。建议在论文中明确标注你的处理方式。5. 数值实验与结果分析5.1 初始参数与博弈场景设定我设计了三组仿真场景来全面检验算法的有效性场景A基准场景追方初始位置在逃方后方500米处两方初始速度均为零。追方加速度上限0.5逃方加速度上限0.4。场景B速度优势场景追方初始相对速度$v_x 1\ m/s$朝向逃方但初始位置距离增加到800米。场景C侧向追击场景追方初始位置不再在逃方正后方而是位于侧向偏移200米的位置测试博弈策略在非对称几何下的表现。所有场景的仿真时长均为30秒时间离散节点数$N 200$。5.2 场景A基准情况下的策略特征场景A的仿真结果显示追方在前5秒内保持中等推力加速约0.3倍最大推力同时逃方迅速启动径向机动试图拉开距离。有趣的是追方并没有一开始就满推力追击。这是因为目标函数中包含了控制能量项——满推力追击虽然能更快缩短距离但消耗的燃料会让总收益下降。从均衡策略的特征来看双方都采用了“先加力后滑行”的策略。前10秒是机动密集区后20秒控制逐渐减小轨迹趋于惯性滑行。这种现象在物理上很好理解末端博弈时间有限初期拉开/接近的优势会在整个剩余时间内被放大所以前期机动价值高后期剩余时间短机动效果的收益递减所以控制量减小。我在仿真结果中提取了终端相对距离为112米这个值略大于拦截判定距离100米也就是说追方在30秒内没能成功拦截但已经很接近了。如果博弈时间延长到35秒追方就能在末端实现拦截。这说明追方的加速度优势0.5 vs 0.4能够扭转初始的位置劣势但需要足够的时间窗口。5.3 场景B和C对称性破坏后的策略变化场景B中追方有初始速度1 m/s的优势但初始距离增加到了800米。博弈结果显示追方策略发生显著变化前期不再冲刺而是保持约0.2倍最大推力的巡航加速利用初始速度逐渐逼近。逃方的应对策略同样有变化——逃方在初始阶段就施加了接近最大值的径向加速度试图尽早拉开距离。从博弈论的视角看逃方是在“利用现状”既然追方已经拥有速度优势逃方必须在前期就全力对抗否则后期很难翻转局面。场景C是侧向追击这个场景最能体现均衡策略的“博弈”本质。初始时追方与逃方的连线方向既不沿径向也不沿切向而是存在一个角度。均衡策略中追方一开始并没有直接朝向当前逃方位置加速而是先沿一个预测的拦截方向飞行逃方则在两个方向之间分配推力试图让自己未来的运动方向更难被预测。最终出来的轨迹不是一条直线追击而是带有明显弯曲的追逐曲线。这个现象让我想起一个经典类比这就像在足球场上防守队员的站位不是盯着球当前位置而是预测球的运动方向提前移动。追逃航天器也是同样的逻辑——好的策略不是追当前位置而是拦截未来的可能位置。5.4 迭代收敛性和ε值分析我在所有场景中都记录了迭代过程中双方策略的变化情况。场景A的迭代收敛最快在第6轮就已经满足$max(\Delta_P, \Delta_E) 10^{-3}$的收敛判据。场景B收敛稍慢到第8轮收敛。场景C最慢需要12轮才收敛——因为侧向几何下双方的策略空间更复杂早期迭代中策略更新幅度较大直到后期才逐渐稳定。ε值的偏离测试结果如下场景追方偏离收益ΔP逃方偏离收益ΔE判定结果A4.6e-46.1e-4通过 (ε1e-3)B7.8e-45.3e-4通过 (ε1e-3)C2.1e-31.8e-3通过 (ε2.5e-3)场景C的偏离收益显著大于前两个场景这说明在几何结构更复杂的场景中迭代求解得到的策略离真正均衡的距离也更远。如果你需要严格的ε-纳什均衡性质场景C需要使用更严格的收敛容差比如$10^{-4}$并增加迭代次数。5.5 控制轨迹的工程可实现性一个在仿真中非常容易被忽视但实际上极其重要的方面是算出来的控制轨迹是否具备工程可实现性。我在结果后处理中发现虽然控制量都在幅值约束内但控制方向在相邻采样点之间可能出现突变。具体来说场景A中追方的推力方向在第17秒到第18秒之间发生了约30度的方向跳变。这种突变在数值上完全合法但在实际航天器上是不可能实现的——推进器的响应延迟和最小脉冲时间决定了控制量不能快速变化。为了模拟这个约束我在后处理中对控制序列做了带宽限制相当于加了一个低通滤波器。这个处理对终端距离的影响不大约3%的变化但让控制轨迹变得更加平滑接近工程可实现的水平。6. 常见问题与调试经验6.1 fmincon不收敛或收敛到明显不合理的解这是我被问得最多的问题。以下是常见原因和解决经验第一初始猜测太差。fmincon是局部优化算法初始猜测直接影响最终结果。如果初始猜测全为零控制很可能收敛到“无控制滑行”这个平凡解。我的经验是先用一个简单的启发式策略作为初始猜测——比如追方初始猜测为“朝着逃方初始位置满推力加速”逃方初始猜测为“背离追方初始位置满推力加速”——这样fmincon更容易找到有意义的解。第二约束尺度不一致。末端距离的量级可能是几十到几百米而控制能量的量级可能是$10^{-3}$到$10^0$。两者相差太大时fmincon的约束处理可能产生数值问题。解决办法是让目标函数中各项的数值量级尽量接近可以通过调整权重系数实现。第三时间离散步长过大。如果$N$取值太小比如50那么BVP离散化带来的误差可能让最优控制解严重偏离真实连续解。建议$N$不低于150200到300之间是比较靠谱的范围。6.2 迭代过程出现振荡在某些参数组合下迭代最优响应会出现振荡即追方策略在“激进”和“保守”之间来回跳变无法收敛。我遇到的一次具体情况是把追方的加速度上限从0.5提到0.8后前5轮迭代还能收敛但从第6轮开始追方策略在“满推力早期冲刺”和“小推力持续追击”之间交替逃方也相应在“全力逃跑”和“延迟机动”之间切换。解决方法是采用阻尼更新。将新策略和旧策略做加权平均u_E_new alpha * u_E_new (1 - alpha) * u_E_old;alpha取0.3到0.5之间时能有效抑制振荡。需要注意的是引入阻尼后收敛速度会变慢迭代次数可能从8轮增加到12到15轮但稳定性显著提升。6.3 bvp4c换成fmincon的坑我在最初复现时使用的是bvp4c求解最优控制问题但很快发现它在处理控制约束幅值饱和时非常不友好。bvp4c的框架天然适合无约束最优控制问题一旦控制量进入饱和区协态方程的结构就会发生改变而bvp4c无法自动处理这种切换。后来我彻底转向了fmincon路线。理由是fmincon天然支持非线性不等式约束控制幅值约束可以直接写在约束函数里而且你不需要为每个问题手工推导最优控制的切换结构。但这个切换也带来了代价fmincon的求解速度比bvp4c慢而且对初始猜测更敏感。我的经验是先用一个粗离散$N50$求解一个粗略解然后将这个解通过插值扩展为细离散$N200$的初始猜测再次求解精细解。这个过程能显著提升fmincon的收敛速度和稳定性。6.4 终端距离计算结果对离散节点数不敏感性的检验在做数值实验时一定要做收敛性检验——否则你不知道你的结果到底是物理规律还是数值假象。我的做法是分别用$N100$、$200$、$400$跑同一组参数比较终端距离和ε值的变化。结果如下$N100$时终端距离为105米$N200$时为112米$N400$时为113.5米。从100到200的变化幅度约7%从200到400的变化幅度约1.3%说明$N200$已经足够收敛$N400$的结果更精确但计算时间增加了近5倍。对于工程实践$N200$是一个合理的平衡点。6.5 一个容易被忽略的数值陷阱参考轨道的线性化误差前面的CW方程是基于圆轨道假设的线性化模型。当博弈终端距离较大或者相对速度较高时线性化假设可能不再成立。我在场景C中做过一个验证用非线性相对运动方程重新仿真均衡策略对应的开环轨迹结果发现终端相对距离与线性化模型计算值有约8米的偏差。对航天器交会来说这个偏差不可忽略。如果你希望仿真结果更接近真实情况有两种路径在规划层就采用非线性模型但这会显著增加求解难度先在线性化模型上求出均衡策略再在非线性模型上进行闭环修正比如引入反馈控制项。第二种路径工程上更可行。我建议在博弈规划完成后添加一个基于线性二次型调节器LQR的跟踪控制器把实际轨迹锁定在规划轨迹附近这样线性化误差的影响会显著减小。7. 从复现到扩展这个项目的应用价值7.1 学术研究层面的延展这个项目在学术上有明确的延展方向。目前我们求解的是开环均衡策略——双方的策略只是时间的函数不依赖于当前状态反馈。如果要走向更接近实战的应用需要考虑闭环反馈策略即策略是当前状态的函数$u_P(X, t)$而非$u_P(t)$。实现闭环有两种思路。第一种是“滚动时域近似”在每个采样时刻重新求解开环均衡只执行第一个采样周期内的控制量然后更新状态并重新计算。这种方法计算量较大但实现相对直接。第二种是“策略函数逼近”利用神经网络近似均衡策略让网络输入当前状态、输出控制指令。这个方向是目前研究的前沿热点但训练数据的生成和收敛性保证都需要仔细设计。7.2 相关领域的应用启发虽然项目在航天器场景下做演示但微分博弈和ε-纳什均衡的数学框架可以迁移到很多其他领域无人机对抗中的追逃博弈和攻防对抗策略设计自动驾驶中自车与对手车辆的交互决策对方不是简单避让而是有博弈性质的决策者多智能体系统中的资源竞争和对抗性策略设计机器人足球比赛中的攻防策略生成在这些场景中目标函数、动力学模型和控制约束都不一样但“通过迭代最优响应来寻找ε-纳什均衡”的核心算法框架是通用的。这也解释了为什么微分博弈在近年来自动驾驶和机器人领域应用越来越广泛——车辆之间的交互本质上就是多重博弈的过程。8. 复现过程中的几点心得前前后后花了两周时间复现这个项目踩了不少坑也积累了一些值得分享的经验。第一点心得是理论文献中那些干净的公式推导和工程实现之间有相当大的距离。比如HJI方程在理论上非常漂亮但真要去数值求解就撞上维数灾难再比如纳什均衡的定义很简洁但数值求解中如何验证均衡性质、如何设置ε阈值都是文献中很少详细讲但实践中必须解决的问题。第二点心得参数选取往往比算法本身更影响最终结果。我花了整整一天调整权重系数、离散节点数和迭代容差才让三组场景都得到满意的均衡解。刚开始调参时一组参数在场景A表现良好换到场景B就出现迭代不收敛后来我发现这是权重系数设置不当导致的调整后三组场景就都能一致地收敛了。第三点心得视觉化是发现问题的关键手段。在代码完成后我认真画了追逃轨迹图和控制量曲线图。轨迹图上如果追方的轨迹出现明显的不合理弯曲比如本可以直飞却绕路大概率是目标函数权重设置有误。控制量曲线上如果出现高频抖动多半是离散点数不够或约束处理有问题。这些直观的图形反馈比盯着终端指标的数字更高效。第四点心得这个项目的代码框架本身就很有复用价值。只要你把dynamics.m中的相对运动方程换成无人机质点模型、把目标函数换成自动驾驶场景下的安全距离指标整个迭代求解引擎可以不加修改地运行。我在复现后顺手把框架用到了另一个模拟任务中验证了这个扩展性。如果后面有机会我会在这个框架上尝试引入闭环策略和神经网络逼近。但就目前的项目版本而言它在“理论分析和工程复现之间搭桥”这个目标上已经是一个相当完整的学习案例了。