航天器追逃博弈的ε-纳什均衡与Matlab数值复现

1. 项目整体拆解与复现思路

1.1 这个项目到底在做什么

我刚开始接触“航天器末端追逃博弈”这个题目时，第一反应是：这不就是经典的微分博弈问题换个马甲吗？但真正复现之后才发现，把博弈论放进航天器的末端交会场景里，事情远没有想象中那么简单——尤其在引入“ε-纳什均衡”这个概念之后，整个问题的求解思路跟教科书上那些线性二次型微分博弈完全是两码事。

先说清楚这个项目要解决的实际问题：一个追踪航天器（追方）试图接近并拦截一个目标航天器（逃方），双方在末端阶段都只剩很短的机动窗口。追方希望最小化末端相对距离和自身燃料消耗，逃方希望最大化末端相对距离同时节省自己的燃料。两边的控制目标直接对立，这就是一个典型的两方零和微分博弈。

这个项目跟你平时看到的“导弹拦截弹道仿真”有本质区别。导弹拦截更多是单边最优控制问题——目标机动是预设的或已知的；而追逃博弈里，逃方是一个具备自主决策能力的智能体，他会根据追方的行为实时调整策略。这就引出了“均衡”的概念：双方都找不到更好的单方面偏离策略时，才能说这个博弈有了一个解。

1.2 ε-纳什均衡和严格纳什均衡的区别在哪

教科书里的纳什均衡要求“没有任何玩家能通过单方面改变策略来获得严格更好的收益”。这是一个精确到小数点后无数位的理想状态。但在工程实践里，你不可能通过数值方法找到一个精度无限高的均衡解——就算找到了，控制执行机构也实现不了那么精细的推力矢量。

ε-纳什均衡放松了条件：只要每个玩家单方面偏离后，收益改善不超过某个给定的小正数ε，这个策略组合就被接受为近似均衡。这个思想放到工程里非常实用，因为卫星的推力器有最小脉冲限制、姿态控制有响应延迟、数值求解有截断误差，这些因素本身就决定了你不可能也不需要找到精确均衡。

这个项目标题里的“ε-纳什均衡”，本质上就是承认工程实现有误差容限，然后用一个可以量化的指标（ε）来评估你所求策略到底离理想均衡有多远。后续仿真里我做的第一件事就是计算实际策略组合的ε值，再根据ε大小判断这个解是否可用——这个环节是理论分析和工程实现之间的关键桥梁。

1.3 复现前的准备：工具链和知识储备

复现这个项目不需要什么特殊的硬件环境，一台普通笔记本就够了。我的工作环境是：Matlab R2022a、Windows 11、CPU为Intel i5-1240P。整个仿真中计算量最大的部分是求解一个带状态约束和控制约束的两点边值问题，在末段追逃时间只有30秒、离散节点200个的情况下，单次求解耗时不到3秒。这是可以接受的。

知识储备方面，建议先弄清楚以下内容再动手：

• 微分博弈的基本定义和Hamilton-Jacobi-Isaacs（HJI）方程的推导过程
• 两点边值问题（BVP）的数值求解方法，特别是打靶法和有限差分法的适用条件
• Matlab自带的bvp4c/bvp5c求解器的用法和局限性
• 非线性规划求解器fmincon的基本用法，特别是如何处理非线性约束
• 航天器相对运动动力学的基本假设，比如Clohessy-Wiltshire（CW）方程的适用条件

说实话，如果是零基础直接上手这个项目，比较吃力。但反过来讲，如果你对上面这些工具都有基本了解，这个项目就是一个把它们串起来的绝佳案例。

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-Isaacs（HJI）方程：

$$
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）的框架：

1. 先猜一个初始的逃方策略$u_E^{(0)}(t)$，比如从当前状态出发，沿远离追方的方向全力机动。
2. 固定$u_E^{(0)}$，求解追方的最优控制问题，得到$u_P^{(1)}$。
3. 固定$u_P^{(1)}$，求解逃方的最优控制问题，得到$u_E^{(1)}$。
4. 重复步骤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 收敛判据和ε值的计算

迭代求解完成后，需要验证所得策略组合确实是ε-纳什均衡。这里我的做法是：对追方和逃方分别做“单方面偏离测试”。以追方为例：

1. 固定逃方使用均衡策略$u_E^*$。
2. 求解追方的最优控制问题（允许追方任意偏离均衡策略），得到偏离后的最优收益$J_P^{dev}$。
3. 计算追方偏离收益与均衡收益的差值：$\Delta_P = J_P^{dev} - J_P^*$。
4. 如果$\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轮才收敛——因为侧向几何下双方的策略空间更复杂，早期迭代中策略更新幅度较大，直到后期才逐渐稳定。

ε值的偏离测试结果如下：

场景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$不低于150，200到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慢，而且对初始猜测更敏感。我的经验是：先用一个粗离散（$N=50$）求解一个粗略解，然后将这个解通过插值扩展为细离散（$N=200$）的初始猜测，再次求解精细解。这个过程能显著提升fmincon的收敛速度和稳定性。

6.4 终端距离计算结果对离散节点数不敏感性的检验

在做数值实验时，一定要做收敛性检验——否则你不知道你的结果到底是物理规律还是数值假象。我的做法是：分别用$N=100$、$200$、$400$跑同一组参数，比较终端距离和ε值的变化。

结果如下：$N=100$时终端距离为105米，$N=200$时为112米，$N=400$时为113.5米。从100到200的变化幅度约7%，从200到400的变化幅度约1.3%，说明$N=200$已经足够收敛，$N=400$的结果更精确但计算时间增加了近5倍。对于工程实践，$N=200$是一个合理的平衡点。

6.5 一个容易被忽略的数值陷阱：参考轨道的线性化误差

前面的CW方程是基于圆轨道假设的线性化模型。当博弈终端距离较大或者相对速度较高时，线性化假设可能不再成立。我在场景C中做过一个验证：用非线性相对运动方程重新仿真均衡策略对应的开环轨迹，结果发现终端相对距离与线性化模型计算值有约8米的偏差。

对航天器交会来说，这个偏差不可忽略。如果你希望仿真结果更接近真实情况，有两种路径：

• 在规划层就采用非线性模型，但这会显著增加求解难度；
• 先在线性化模型上求出均衡策略，再在非线性模型上进行闭环修正（比如引入反馈控制项）。

第二种路径工程上更可行。我建议在博弈规划完成后，添加一个基于线性二次型调节器（LQR）的跟踪控制器，把实际轨迹锁定在规划轨迹附近，这样线性化误差的影响会显著减小。

7. 从复现到扩展：这个项目的应用价值

7.1 学术研究层面的延展

这个项目在学术上有明确的延展方向。目前我们求解的是开环均衡策略——双方的策略只是时间的函数，不依赖于当前状态反馈。如果要走向更接近实战的应用，需要考虑闭环反馈策略：即策略是当前状态的函数$u_P(X, t)$而非$u_P(t)$。

实现闭环有两种思路。第一种是“滚动时域近似”：在每个采样时刻重新求解开环均衡，只执行第一个采样周期内的控制量，然后更新状态并重新计算。这种方法计算量较大，但实现相对直接。第二种是“策略函数逼近”：利用神经网络近似均衡策略，让网络输入当前状态、输出控制指令。这个方向是目前研究的前沿热点，但训练数据的生成和收敛性保证都需要仔细设计。

7.2 相关领域的应用启发

虽然项目在航天器场景下做演示，但微分博弈和ε-纳什均衡的数学框架可以迁移到很多其他领域：

• 无人机对抗中的追逃博弈和攻防对抗策略设计
• 自动驾驶中自车与对手车辆的交互决策（对方不是简单避让，而是有博弈性质的决策者）
• 多智能体系统中的资源竞争和对抗性策略设计
• 机器人足球比赛中的攻防策略生成

在这些场景中，目标函数、动力学模型和控制约束都不一样，但“通过迭代最优响应来寻找ε-纳什均衡”的核心算法框架是通用的。这也解释了为什么微分博弈在近年来自动驾驶和机器人领域应用越来越广泛——车辆之间的交互本质上就是多重博弈的过程。

8. 复现过程中的几点心得

前前后后花了两周时间复现这个项目，踩了不少坑，也积累了一些值得分享的经验。

第一点心得是：理论文献中那些干净的公式推导和工程实现之间有相当大的距离。比如HJI方程在理论上非常漂亮，但真要去数值求解就撞上维数灾难；再比如纳什均衡的定义很简洁，但数值求解中如何验证均衡性质、如何设置ε阈值，都是文献中很少详细讲但实践中必须解决的问题。

第二点心得：参数选取往往比算法本身更影响最终结果。我花了整整一天调整权重系数、离散节点数和迭代容差，才让三组场景都得到满意的均衡解。刚开始调参时，一组参数在场景A表现良好，换到场景B就出现迭代不收敛；后来我发现这是权重系数设置不当导致的，调整后三组场景就都能一致地收敛了。

第三点心得：视觉化是发现问题的关键手段。在代码完成后，我认真画了追逃轨迹图和控制量曲线图。轨迹图上如果追方的轨迹出现明显的不合理弯曲（比如本可以直飞却绕路），大概率是目标函数权重设置有误。控制量曲线上如果出现高频抖动，多半是离散点数不够或约束处理有问题。这些直观的图形反馈，比盯着终端指标的数字更高效。

第四点心得：这个项目的代码框架本身就很有复用价值。只要你把dynamics.m中的相对运动方程换成无人机质点模型、把目标函数换成自动驾驶场景下的安全距离指标，整个迭代求解引擎可以不加修改地运行。我在复现后顺手把框架用到了另一个模拟任务中，验证了这个扩展性。

如果后面有机会，我会在这个框架上尝试引入闭环策略和神经网络逼近。但就目前的项目版本而言，它在“理论分析和工程复现之间搭桥”这个目标上，已经是一个相当完整的学习案例了。