不完全信息ε-Nash均衡下的航天器追逃博弈Matlab复现

复现期刊论文这件事最让人头疼的从来不是公式有多长而是从公式到代码之间那一段没人替你踩的坑。这次复现的项目是不完全信息ε-Nash均衡下的航天器末端追逃博弈策略Matlab代码实现。它解决的具体问题是这样的追踪航天器在末端阶段接近逃逸航天器时双方对彼此的推力性能并不是完全清楚的各自需要在“不知道对方真实底牌”的前提下选择相对最优的控制策略而最终允许找到一个近似范围内的均衡解。这个动态过程本质上是微分博弈加入了不完全信息后又变成贝叶斯博弈论文里给出的均衡概念从精确纳什均衡松绑成了ε-Nash均衡。项目适合三类人参考做航天器轨道博弈仿真的研究生、想学微分博弈数值解法的工程师、以及需要把不完全信息博弈方法落到仿真代码里的开发者。下面从建模、求解、代码实现到调参排查完整过一遍。1. 项目复现前的三个关键认知这类题目看着很唬人但拆开看无非是三件事追逃博弈怎么建模型、不完全信息怎么处理、ε-Nash均衡怎么落地。我在动手写代码之前先花了很长时间把这三个概念在脑子里理顺了因为后面的每个函数和数据流都建立在这三件事的理解上。1.1 追逃博弈到底在求解什么追逃博弈属于微分博弈的一种。你有一个追踪者P和一个逃逸者E两者的状态都由一组微分方程描述各自通过控制输出去影响状态的演化。追踪者的目标是让终端相对距离尽可能小逃逸者的目标正好相反这构成了一个零和博弈追方极小化终端距离逃方极大化终端距离。在航天器末端场景里状态通常是追踪器和逃逸器之间的相对位置、相对速度。控制输入是推进器产生的加速度带有明显的饱和约束。整个博弈在固定末端时刻内完成追踪器要在这段时间内尽量逼近逃逸器要尽量拉开距离。两端目标相反中间靠控制博弈这就是微分博弈最典型的鞍点求解问题。我第一次理解这个概念的时候直接用了一个生活化的类比两个人站在一个对称的迷宫里对跑一个人要把两人距离缩到最小一个人要把距离拉大到最大双方都能看到对方当前的动作同时不断调整自己的路线。每一时刻的选择都同时影响对方下一步的选择这不是“我先走你后走”的静态博弈而是连续时间里的动态对抗。如果双方信息完全透明、控制能力完全已知那存在一个精确的鞍点策略也就是极小极大意义下的均衡解。但论文里这个模型还有一个前置条件末端。也就是说这个博弈不是从很远的距离就开始而是双方已经进入一个相对较近的空域博弈窗口有限。末端博弈的特点是时间短、动态变化快、初始状态敏感这种场景下拿到一个理论上严格最优的控制率不如拿到一个快速可算的近似均衡策略来得实在。这也是后面引入ε-Nash均衡的实际动机之一。1.2 不完全信息为什么要引入贝叶斯框架所谓不完全信息在这个题目里指的是双方对对方的某些关键参数不是完全掌握的。最常见的情况是不知道对方的推力上限。你可能知道对方是某个型号级别但具体到它在末端还能给出多大的加速度存在不确定性。这就像打扑克的时候你知道对方手上有牌但不知道他具体的牌面是什么。如果你假设对手总能打出最大推力那你的策略会偏向过于保守如果你假设对手很弱又可能在关键时刻被突破。实际工程中推力器性能、燃料剩余量、制导模式都属于不透明的信息指望完全精确获知是不现实的。处理这个问题博弈论里有一套标准做法就是海萨尼转换把不完全信息转化为不完美信息。引入一个虚拟的“自然”行动它先以某个概率分布抽取每个参与人的类型每个参与人知道自己的类型但只对手的类型有一个先验信念。在这个框架下面每个参与人最大化的是自己的期望支付均衡概念也随之从纳什均衡变成贝叶斯纳什均衡。这里有一个很容易走偏的思维有人会问为什么不用鲁棒控制那种最坏情况设计来做如果是鲁棒控制你会假设逃逸者永远以最强状态出现这样得到的追踪策略确实“绝对安全”但代价是极其保守可能在大多数实际场景里性能都很差。不完全信息博弈的思路不同它按照各类型的概率分布来权衡既不是对最坏情况投降也不是对平均情况无防备而是通过期望支付来引导策略选择。航天器末端博弈里燃料、推力都是有限资源过度保守等于浪费所以这个框架在实际场景中比纯鲁棒更合理。1.3 ε-Nash均衡和精确均衡差在哪教科书里的纳什均衡定义是没有任何一个参与人可以通过单方面偏离来获得更高支付。但在数值实现里这个标准几乎是不可达的。原因很简单连续控制空间、非线性动态方程、加上不完全信息的类型集合要在一个迭代算法里验证“任何偏离都无法获益”这句话代价非常高。ε-Nash均衡本质上是把这个严格条件松绑了每个参与人单方面偏离能获得的最大收益不超过ε。也就是说如果追方保持当前策略逃方无论怎么换策略最多比自己现在多拿ε的支付那这个策略组合就已经够用了。这个ε不是数值误差是你明确给出的容差线。工程上这个概念的吸引力在于它给了数值算法一个明确的停机条件。我跑迭代的时候不需要证明“真的没有更好的策略了”只需要验证“偏离的收益上界小于ε”就可以停手。复现论文的时候重点也不是去索要一个严格的均衡点而是按照论文的设定在合理的迭代次数内拿到一个满足ε条件的策略组合并验证它对应的终端距离是符合预期的。下面用一个表来区分这两类均衡在实际代码中的差别对比项精确纳什均衡ε-Nash均衡偏离收益条件任何偏离收益严格为零偏离收益小于等于ε数值停机条件几乎无法严格验证可量化、可验证迭代收敛难度高容易震荡相对宽松实用性好工程使用价值理论完美但难落地近似够用适合仿真分析对复现项目来说ε取多大直接决定策略质量和计算时间。我后面在调参部分会展开讲这个取值经验这里先记住一个原则ε不能太小否则迭代很难收敛也不能太大否则均衡没有任何意义一般要根据量纲先看终端距离的一个可接受偏差是多少再反推ε的量级。2. 建模细节从论文公式到可计算模型这一节是复现里花费时间最多的地方。论文里的公式是从物理规律出发写的但代码里的模型必须明确状态变量是什么、单位是什么、约束怎么表达。哪怕一个符号理解偏差都会导致结果完全对不上。2.1 相对运动方程与无量纲化航天器末端追逃通常选取参考轨道上的目标点建立当地轨道坐标系也就是LVLH坐标系。x轴指向远离地心的径向y轴沿迹向z轴垂直于轨道平面。在这个坐标系下两航天器之间的相对运动可以用C-W方程描述。近圆轨道假设下相对运动方程是线性的径向ẍ 3n²x 2n·ẏ a_x迹向ÿ -2n·ẋ a_y法向z̈ -n²z a_z其中n表示参考轨道角速度a_x、a_y、a_z是控制加速度分量。这里要注意C-W方程只适用于近圆轨道、两星距离远小于轨道半径的末端场景。如果参考轨道偏心率明显或者相对距离到了几公里的量级线性化方程误差会迅速放大那就要换成Tschauner-Hempel方程甚至J2摄动下的模型了。这次复现的论文场景在末端阶段距离短、时间短C-W方程是合理选择。但这里有个数值问题。航天器相对距离可能是几百米相对速度可能是每秒几米轨道角速度n对于近地轨道大概是0.001量级各项数值之间差了好几个数量级。直接拿这个单位体系去算最优控制优化器很容易因为条件数太大而出现病态。常规做法是做无量纲化。取特征长度L为初始相对距离取特征时间1/n把所有变量除以对应量纲的特征量让状态量全部落在0.1到10的尺度范围。无量纲化之后C-W方程形式基本不变但是控制加速度的约束上限也会相应缩放。这一步看起来不起眼实际上决定了后面fmincon或者打靶算法能不能稳定收敛。我复现时一开始没注意这个问题直接用的是物理单位迭代经常出现步长过小、目标函数不平滑的现象无量纲化之后基本就好了。状态方程对应的Matlab函数很直接function sdot cw_dynamics(s, u_p, u_e, n) % 输入 % s : 状态向量 [x; y; z; vx; vy; vz] % u_p: 追踪者控制加速度 [3x1] % u_e: 逃逸者控制加速度 [3x1] % n : 参考轨道角速度 % 输出 % sdot: 状态导数 [6x1] x s(1); y s(2); z s(3); vx s(4); vy s(5); vz s(6); ax 3*n^2*x 2*n*vy u_p(1) u_e(1); ay -2*n*vx u_p(2) u_e(2); az -n^2*z u_p(3) u_e(3); sdot [vx; vy; vz; ax; ay; az]; end这个方程是后面所有求解函数的核心。注意控制输入即使在物理上属于不同航天器但在相对运动方程里是两个独立控制向量的叠加。追逃双方分别在自己的推力约束范围内调整u_p和u_e博弈就在这两组控制输入之间展开。2.2 支付函数和终端约束怎么定学术论文里的支付函数经常写得很抽象落地到代码首先得明确是极大化还是极小化。零和博弈下追踪者极小化终端距离逃逸者极大化同一个量一个简单又常见的支付函数是J α·||r(t_f)|| β·∫(||u_p||² ||u_e||²)dt第一项是终端相对距离第二项是控制能量惩罚。α和β是权重。追踪者希望J越小越好逃逸者希望J越大越好。实际写代码时逃逸者一侧的求解通常是把追踪者的目标函数取负号然后做极小化。这里有一个我自己开始没想明白的点支付函数里的燃料项不是必须的。如果博弈时间很短燃料消耗相对次要可以直接只用终端距离。但加了燃料惩罚之后控制输入曲线不会出现那种剧烈的边界跳变数值上更平滑迭代更容易收敛。复现论文的时候最好先把支付函数写成带权重参数的形式先把β设成很小或者零跑通逻辑再加进去微调。控制约束是另一件必须明确的事。追逃双方推力有限每个轴的加速度都有上下界。常见写法是u_min ≤ u_i ≤ u_max, i x, y, z也就是分轴约束。也可以写总推力约束||u|| ≤ u_max但那种约束在直接法里是非线性约束求解成本高。我在复现里用了分轴约束因为它和论文描述一致而且在fmincon里只是边界条件不需要额外写非线性约束函数。末端时刻t_f也是关键参数。末端追逃博弈一般设定固定末端时刻也就是卫星发动机工作的一个时间窗口。这个窗口不能太短否则还没看出博弈效果就结束了也不能太长否则多轮迭代算得很慢。论文里通常会给一个参考值我复现的时候会在这个基准附近扫几组值看均衡策略下的距离曲线是否平滑。2.3 类型空间离散化处理不完全信息意味着逃逸者的某一个参数是随机的。在这个项目里通常是推力上限u_max_e。连续分布直接求解很难因为类型空间无限每个类型对应一个不同的博弈没法直接算。标准做法是离散化。把推力上限的取值范围切成几个代表点比如低、中、高三种类型并赋予对应的先验概率类型1推力上限低概率0.25类型2推力上限中概率0.5类型3推力上限高概率0.25为什么用三个类型而不是两个四个我试下来两个类型太少策略区分度不够四个类型会让迭代里的最优响应计算次数增加每个类型都要跑一遍完整的最优控制问题耗时上升明显。三个类型配合不对称概率分布已经能体现出不完全信息的核心效果。追踪者不知道逃逸者到底是哪个类型所以它的最优响应要基于期望支付来计算。期望的形式就是把三个类型对应的终端距离按概率加权求和。逃逸者知道自己的真实类型所以对每个类型单独求自己的最优响应。这种“追方按期望、逃方按真实类型”的不对称关系正是贝叶斯博弈在代码层面最核心的体现。类型空间概率分布的取值也会影响结果。概率越集中在某个类型均衡策略越接近该类型已知的完全信息博弈。概率分布均匀时追踪者会更稳妥不会对任何一种类型做过度反应。这一块属于参数敏感性分析的内容后面会有专门的小节。3. Matlab代码实现的核心模块模型一旦清楚代码结构其实很常规。整个项目拆成几层主脚本负责参数设置和流程控制函数模块负责状态方程、支付计算和最优响应求解最后是绘图模块。没有复杂的工程架构但每个文件职责要分明。3.1 程序骨架与数据流我用的目录结构是这样的main.m总入口设置所有物理参数、博弈参数、类型概率调用均衡求解流程最后出图dynamics/cw_dynamics.mC-W相对运动方程cost/payoff_total.m计算完整支付函数内部调用状态积分solver/solve_best_response.m给定对手策略用fmincon求当前参与人的最优控制solver/equilibrium_iteration.m迭代求ε-Nash均衡的主循环plot/plot_result.m绘制相对轨迹、控制曲线、终端距离收敛曲线数据流是main初始化参数把所有参数放进一个struct里传给equilibrium_iteration。迭代函数内部初始化一组控制策略然后交替调用solve_best_response每次求解完更新策略检查收敛条件。收敛后把最终策略和仿真轨迹返回给mainmain再做一次完整仿真并绘图。这里有一个经验不要把参数散在各个函数里一定要打包成一个params结构体。因为到后面调参的频率很高如果每个函数都有自己的硬编码参数改起来会非常痛苦。我后来把无量纲化比例系数也放进了params切换物理单位时只需要改main里的三个参数。3.2 核心函数状态方程、支付计算、最优响应求解状态方程已经写过了重点是payoff_total怎么实现。它做的事情就是给定双方控制策略从初始状态出发积分C-W方程得到终端相对位置然后算终端距离和控制能量项。积分我用的是ode45但需要注意ode45在博弈迭代的DIQK内部会被调用很多次速度并不快。如果时间节点很多可以考虑改成定步长RK4或者把状态方程向量化能省下不少时间。求解最优响应的时候可以用直接法或者间接法。直接法把控制序列离散成变量交给fmincon去优化好处是思路直观对初值不太敏感坏处是维度高、计算慢。间接法通过伴随变量把最优控制问题转成边值问题维度低、精度高但需要推导哈密顿函数和伴随方程初值猜不准就发散。这次复现代码里我采用的是直接法来保证稳定性框架大致如下function u_opt solve_best_response(opponent_u, player_type, params, is_pursuer) % 把控制序列离散为N个时间点 N params.N; u0 zeros(3, N); % 初始猜测 % 约束条件 lb repmat(params.u_min, 1, N); ub repmat(params.u_max, 1, N); % 定义优化目标函数 cost_func (u_vec) payoff_player(u_vec, opponent_u, player_type, params, is_pursuer); % 求解 options optimoptions(fmincon, Algorithm, sqp, ... MaxIterations, 300, OptimalityTolerance, 1e-6); u_opt fmincon(cost_func, u0, [], [], [], [], lb, ub, [], options); end这个函数里面有几个细节值得说。一是u0初始猜测如果对手策略是追踪者最优那逃逸者求解时初始猜测可以直接用上一轮的对手策略加一个小扰动收敛速度会快很多。二是fmincon的算法选sqp在约束较简单时比内点法稳定而且不容易产生不可行解。三是每次调用fmincon之后最好把返回的fval保存下来作为后续收敛判断的一部分。还有一个常见坑控制序列是矩阵但fmincon要求决策变量是向量。所以函数内部要做一次reshape把3×N的控制矩阵展开成长度为3N的列向量。计算支付时再reshape回去。这个看似琐碎但很多初学者在这里会写得很混乱。3.3 迭代求ε-Nash均衡的实现思路迭代求均衡思路不复杂本质上是交替求最优响应。我用的迭代流程写成伪代码是这样的% 初始化双方策略 u_p init_control(zeros(3,N), params); u_e_type {u_e1, u_e2, u_e3}; % 每个逃逸者类型各一个策略 for iter 1:max_iter % 第一步给定逃逸者各类型策略更新追踪者策略 u_p_new solve_pursuer_best_response(u_e_type, params); % 第二步给定追踪者策略更新每个逃逸者类型的策略 for type 1:3 u_e_type{type} solve_evader_best_response(u_p_new, type, params); end % 计算当前策略组合下的支付 [J_p, J_e] evaluate_payoffs(u_p_new, u_e_type, params); % 计算偏离收益用当前策略的支付减去最优响应的支付 dev_p J_p_best - J_p_current; dev_e max(J_e_best_type - J_e_current_type); % 判断是否满足ε-Nash if abs(dev_p) eps_nash abs(dev_e) eps_nash break; end % 更新策略加入松弛因子防止震荡 u_p alpha * u_p_new (1 - alpha) * u_p; end这个流程里最关键的是最后一步松弛。如果不加松弛直接让u_p等于u_p_new经常会出现来回跳变追踪者刚适应逃逸者逃逸者又换了一套策略循环永远收敛不了。我实测下来松弛因子取0.5到0.7之间最稳也就是每一轮只往最优方向移动一半左右。虽然迭代次数会多几轮但振荡会明显减少。收敛判定里面有一个容易被忽略的点要在每次迭代里保存当前策略与上一轮策略的差值作为另一个收敛指标。ε-Nash条件衡量的是“偏离没有好处”而策略差值衡量的是“策略不再变化”。两者配合使用能避免一种尴尬的情况策略几乎没动、但偏离条件还没满足于是白白多跑几十轮。逃逸者类型是离散的三个追踪者的最优响应必须考虑期望支付。这个期望支付可以写成J_p_expected w1·J_p_type1 w2·J_p_type2 w3·J_p_type3其中w是类型概率J_p_typei是在逃逸者为第i类时的终端支付。这个期望支付函数写进payoff_player里追踪者调用的时候传入整个逃逸者策略集合逃逸者调用时只需要传入自己对应类型的策略。3.4 仿真与可视化验证拿到均衡策略组合之后不能直接收工还要做一次完整的开环仿真来验证结果。把均衡策略代入状态方程从初始状态积分到终端时刻记录每一时刻的相对位置和相对速度然后画三张图相对运动轨迹、控制加速度曲线、终端相对距离随迭代的变化。相对运动轨迹图最能说明博弈效果。追踪者的轨迹会呈现出一种逐步逼近的形态逃逸者的轨迹则是尽量远离追踪者的方向。不完全信息的影响在这个图里很直观追踪者的轨迹不是针对某一个特定逃逸者类型设计的而是针对三种类型的一个折中所以它看起来不会特别激进初始阶段甚至会显得有点“犹豫”。控制加速度曲线可以验证约束是否起作用。理想情况下控制输入在大部分时间处于饱和状态末端时刻附近才有明显的切换。如果控制曲线出现高频抖动那说明fmincon没有收敛或者时间离散步长太大需要回到参数设置去检查。终端相对距离的收敛曲线是判断迭代是否成功的直接依据。稳定的迭代会呈现一种先快速下降、然后逐渐平稳的趋势。如果曲线一直在波动那基本可以断定是松弛因子太大或ε标准过严。另外我会做一个额外的对比实验把逃逸者类型当作已知信息重新求解一个完全信息博弈得到一组“理论最优”轨迹然后把不完全信息博弈的均衡策略放到三种逃逸者类型下分别仿真比较终端距离。这个对比能定量说明不完全信息带来的性能损失也是论文里常见的一张结果图。复现的时候这张图特别重要画出来之后基本可以对着论文里的数值来校验。4. 复现中踩过的坑与排查技巧复现论文的过程远没有公式推导那么顺畅。我卡住的时间主要花在收敛问题、参数敏感性以及与论文结果对不上这三个方面。这里把典型的坑和排查方法整理出来希望能省掉一些弯路。4.1 收敛性问题的常见根源fmincon不收敛是最常见的问题但原因往往不在fmincon本身。我排查收敛问题有一个固定顺序先看初值再看比例尺再看约束。初值问题是我踩得最多的。fmincon是个局部优化算法初始猜测离最优解太远就会陷进糟糕的局部极值。解决办法是分阶段求解先把支付函数里的终端距离项权重调大、燃料项调为零用一组合适的均匀控制作为初值求出结果后再把燃料项加上用上一轮的解作为新初值继续优化。这种“逐步升温”的思路在非线性最优控制里很常用。比例尺问题前面提过就是没有做无量纲化。如果你发现目标函数值在迭代过程中差不多是10^6量级而步长却是10^-4量级那么大概率是单位体系出了问题。把所有状态和控制输入无量纲化之后优化器的判别条件才有意义。还有一个隐蔽问题是时间离散步长。直接法里控制量是按时间节点离散的如果节点太少控制表达不出来如果节点太多fmincon的维度上升计算慢且容易不收敛。我自己常用的做法是先取20个节点跑通再看结果够不够平滑不够就加到50个。论文结果里的控制曲线通常是分段光滑的用50个节点基本能表达清楚。4.2 参数敏感性分析与调参经验ε-Nash均衡的容差ε是最关键的参数。根据公式量纲终端距离是以米为单位的如果允许的偏差是几十厘米量级那么ε取1e-2就是在合理范围内的。我复现的时候先用一个较大的ε快速跑出整体策略形态确认形态合理后再逐步减小ε。如果一上来就用1e-6很可能跑几十轮都收敛不了白白浪费时间。类型概率分布也是敏感参数。我试过把概率改成均匀分布和极不对称分布追踪者的初始策略差别很大。极不对称分布相当于追踪者认定逃逸者就是那个高概率类型策略会和完全信息博弈下的最优策略很接近。这个现象可以用来做验证如果代码逻辑正确那么把某个类型概率设为0.99时追踪者策略应该和完全信息下针对该类型的最优策略高度一致。用这个方法来检查代码是很有价值的。控制约束上限同样值得做几组扫描。推力上限越高博弈的对抗强度越大终端相对距离会明显下降。如果扫描时发现终端距离对推力上限不敏感那很可能是约束没有真正激活检查一下求解出来的控制序列是不是大部分时间都远低于上限。控制序列如果一直低于上限说明要么博弈时间过长、要么支付函数权重有问题。还有一个容易被忽略的参数是初始相对状态。末端博弈对初始相对速度非常敏感初始速度稍微偏向远离方向追踪者就要多花好几个时间步来“调头”。复现论文的时候初始状态一定要按论文表格里的数值精确设置哪怕差了1%结果图都可能对不上。4.3 与论文结果对不上的排查清单复现过程中最熬人的阶段就是仿真结果和论文对不上。我整理了一份排查清单按顺序走一遍能定位大多数问题。现象可能原因排查与解决办法终端距离整体偏大支付函数符号反了检查追踪者是否在极小化逃逸者是否在极大化轨迹形态完全不同初始状态单位错误对照论文的初始状态单位确认不是公里和米的混用控制曲线持续高频振荡时间步长太大或初值不光滑缩短离散步长增加时间节点数迭代发散松弛因子过大或ε标准过严降低松弛因子到0.5临时放宽ε追踪者策略过于激进类型概率设置不当检查概率分布尝试向中等类型集中逃逸者各类型响应几乎一样类型离散间隔太小增大推力上限的档位间隔数值规模异常大未做无量纲化引入特征长度和特征时间重新缩放状态量迭代策略来回跳变缺少松弛项按0.6比例混合新旧策略后更新排查的时候一定要保持一个原则每次只改一个参数。我一开始没守住这个原则同时改了权重、离散节点数和松弛因子结果根本判断不了是哪个改动让结果发生变化的。后来改成每次只动一个变量再配上结果图的对比定位问题的速度一下子快了很多。还有一个心得是复现论文时不要指望复现出“一模一样”的曲线。期刊论文里的结果图往往是经过大量调参后的代表性结果数字复现到同一量级、趋势一致、关键现象能重现就已经是合格的复现了。完全一致的数值结果只有在所有参数、初值、求解器设置都完全相同的情况下才可能做到而论文常不会公开这些细节。我个人在实际操作中还发现一个有用的做法把每次调参后的结果自动存档包括参数表、三个关键图和一个txt格式的日志文件。日志里自动记录时间和参数哈希这样哪天改乱了还能快速回滚到某个历史版本。这个习惯在后期参数敏感性分析阶段帮了我大忙没有它我会陷入反复对比旧图的混乱里。这个项目做完之后的扩展空间其实也不小。后续可以考虑把类型空间从离散类型推广到连续分布用数值积分替代加权求和也可以把固定末端时刻改成追方满足某个脱靶量阈值后自动终止的自由末端问题更进一步的玩法是把C-W方程换成更精确的二体相对运动模型验证近似均衡在非线性模型下是否仍然稳健。每一步扩展都能让这个复现项目在原有论文基础上再往前走一步对于需要发论文或者做课题积累的人来说都是值得投入的方向。