简介这是基于扩展卡尔曼滤波(EKF)的航天器姿态轨道耦合控制系统Matlab实现面向航天器动力学与控制方向的高年级本科生、研究生及科研工作者适用于课程设计、专题研究与毕业设计采用模块化参数化编程覆盖姿态确定、轨道动力学、EKF状态估计、PD/LQR控制等核心环节并适配轨道维持、交会对接等任务场景支持在2014a、2019b、2024b等版本直接运行。资源包共89个文件包含30个m源码、19个mat数据、15个bmp仿真曲线、2个slx/Simulink模型及14个zbak备份、README说明等整体16.53MB已有43人学习下载。通过随附案例数据集与仿真图像读者可深入理解EKF在非线性耦合系统中的状态估计实现、控制律调试与结果分析涵盖姿态误差、轨道偏差、控制力矩等关键指标完整源码和注释为后续修改扩展提供基础可直接支撑课程项目或学位论文的仿真平台是航天器控制教学与科研实践的高性价比工具。1. 基于EKF的航天器姿态轨道耦合控制这个问题为什么值得自己做一遍做航天器GNC的人大概率遇到过这个场景轨道机动时推力器点火星体被偏置力矩顶歪姿态机动时推力方向跟着偏轨道又被扰动。分开设计的姿态环和轨道环在仿真里各自收敛一联调就互相甩锅。把姿态和轨道放进同一个状态向量做EKF估计再用耦合补偿的控制律收回来是解决这个问题的正路。这篇笔记从耦合建模、EKF实现、LQR控制律到Matlab代码一步步拆适合正在做预研或毕业设计的工程师看完可以直接把代码跑起来。2. 姿态轨道耦合动力学建模状态方程、坐标约定与Matlab里的离散化2.1 三个必须进模型的耦合项航天器在轨的工作点通常是“姿态对地定向轨道保持”这种情况下姿态和轨道之间的耦合相对弱把两个回路分开设计还能勉强工作。一旦进入轨道机动、交会接近或者大角度姿态机动场景下面三个耦合项就不能忽略了。第一推力矢量与质心不重合。每个推力器都有安装位置和方向推力并不严格通过质心会产生固定力臂的干扰力矩 M_dist d × F。这个力矩在轨道机动时持续存在姿态控制系统必须实时补偿否则星体姿态被推跑。第二姿态偏差会改变推力方向。推力指令在体坐标系下表达如果姿态有误差 δθ惯性系里实际推力方向就和期望方向差了 δθ。对轨道控制来说这会直接造成速度增量方向偏差是典型的“姿态误差转化为轨道误差”。第三重力梯度力矩与轨道位置相关。非球形引力场里轨道高度和经纬度变化都会改变重力梯度力矩而姿态控制又把这种力矩当干扰。耦合趋势在低轨道、大惯量差的卫星上更明显。把这三个耦合项写进模型EKF和控制律才能看到“一个卫星”而不是“两个独立回路”。2.2 状态向量为什么是13维笛卡尔坐标与四元数的选择常用状态向量取 x [r; v; q; ω]其中 r、v 是惯性系下的位置和速度各3维q 是星体姿态四元数4维ω 是体坐标系角速度3维合计13维。轨道用笛卡尔坐标而不是轨道根数是因为EKF预测方程在笛卡尔坐标下写起来直观、雅可比好算轨道根数在接近圆轨道时存在奇异没必要给自己找麻烦。四元数有单位模约束不能像普通向量那样做加权平均和协方差运算。这也是后续用误差状态EKFMEKF而不是直接状态EKF的前提。用欧拉角做状态可以绕开约束但大角度机动时会遇到万向锁姿态估计直接翻车工程上基本没人拿欧拉角跑EKF。四元数这里统一采用标量在前约定 [q0; q1; q2; q3]与Matlab自带工具包的顺序不一致跑之前用单位旋转验证一下。2.3 耦合动力学Matlab函数从惯性系到体坐标系的完整代码function xdot coupled_dynamics(x, u, p) % 航天器姿态轨道耦合动力学 % x [r(3); v(3); q(4); w(3)] % u [F(3); M(3)] 惯性系推力(N)、惯性系控制力矩(N*m) % p 参数结构体: mu, J, d_cm(推力作用点相对质心) r x(1:3); v x(4:6); q x(7:10); w x(11:13); q q / norm(q); % 四元数保持单位模 R quat2dcm(q); % 体坐标-惯性坐标按工具包约定调整 F_body R * u(1:3); % 推力转体坐标系 M_body R * u(4:6) cross(p.d_cm, F_body); % 轨道: 二体引力(可扩展J2项) g -p.mu * r / norm(r)^3; % 姿态运动学: qdot 1/2 * Omega(w) * q Omega [0, -w(3), w(2); w(3), 0, -w(1); -w(2), w(1), 0]; qdot 0.5 * [-q(2:4) * w; q(1)*eye(3)*w Omega * q(2:4)]; % 姿态动力学: J wdot w x (J w) M_body wdot p.J \ (M_body - cross(w, p.J * w)); xdot [v; g; qdot; wdot]; end逻辑说明代码分四段。轨道加速度目前只用了二体引力需要J2摄动时在 g 那一行追加即可不影响EKF接口。姿态运动学用的是四元数微分方程其中 Omega 矩阵实现叉乘运算。姿态动力学里 M_body 除了控制力矩还加了 cross(p.d_cm, F_body)这就是推力偏置力矩是整个“耦合”在模型里的具体体现删掉这一项后面所有控制补偿都无从谈起。参数说明p.mu 是地球引力常数常用值 3.986004418e14 m^3/s^2p.J 是 3×3 转动惯量矩阵p.d_cm 是推力作用点相对质心的位置向量取值建议在几厘米到几十厘米量级——太小会让偏置力矩淹没在数值噪声里太大则质心标定可能有问题。实际任务里这个值来自地面质量特性测试误差通常在2mm以内仿真里给固定值即可。2.4 用RK4把连续动力学离散成仿真步EKF预测步需要求解状态转移仿真主循环也需要对动力学积分。固定步长四阶Runge-Kutta是航天仿真最常用的做法在秒级步长时精度足够、计算量可控。function x_next rk4_step(dynfun, x, u, p, dt) % 固定步长RK4传播一个仿真步 k1 dynfun(x, u, p); k2 dynfun(x dt/2*k1, u, p); k3 dynfun(x dt/2*k2, u, p); k4 dynfun(x dt*k3, u, p); x_next x dt/6 * (k1 2*k2 2*k3 k4); endRK4截断误差是 O(dt^5)对LEO轨道秒级步长完全够用。注意 u 在整个步长内保持不变这是零阶保持假设如果控制律更新频率低于动力学积分频率就需要在主循环里管理好几个时钟后面第4章的主循环会演示怎么处理。大偏心率轨道或者近地点剧烈机动时固定步长会吃亏解决办法是改 ode45 跑真值、EKF预测步仍然固定步长这一点在第5章避坑详细展开。3. EKF状态估计设计误差状态、预测步与观测更新的Matlab实现3.1 误差状态EKF回避四元数约束的工程做法直接对13维状态跑EKF会遇到数学问题四元数有单位模约束EKF的预测和更新都假设状态服从高斯分布加入四元数后分布不再高斯协方差更新会把单位模约束破坏掉。预测完归一化一次、更新完再归一化一次短期跑起来还行一旦观测噪声变大或机动增强P矩阵就会失去正定性滤波性能断崖式下降。工程标准做法是误差状态EKF主状态仍用13维但P矩阵只作用在12维误差状态 δx [δr; δv; δθ; δω] 上。姿态误差用三维角度向量 δθ 而不是四元数增量协方差更新始终在无约束向量空间里进行四元数只在主状态里传播和注入约束与更新分开既满足高斯假设又保住四元数几何结构。误差状态EKF又叫MEKFMultiplicative EKF因为姿态误差用四元数乘法定义真实姿态 q δq ⊗ q_predδq 是小角度四元数与 δθ 的关系是 δq ≈ [1; δθ/2]。航天工程里基本都用这套。3.2 预测步实现数值雅可比是偷懒但稳妥的选择误差状态EKF预测分两步先让主状态沿非线性动力学走一步再让误差状态协方差沿线性化矩阵传播。线性化矩阵 F ∂f_err/∂δx 理论上可以推解析雅可比但姿态轨道耦合下推导容易出错、模型一改就要重推我更常用数值雅可比。差分变量是12维误差状态四元数方向做扰动时要先转成小角度四元数再乘到主状态上这是和普通数值差分最大的区别。function [x_pred, P_pred] mekf_predict(x, P, u, p, dt, Qd) % MEKF预测步: 标称状态传播 误差状态协方差传播 % Qd是离散化后的过程噪声协方差(12x12) f (xx) rk4_step(coupled_dynamics, xx, u, p, dt); % 1. 标称状态传播 x_pred f(x); x_pred(7:10) x_pred(7:10) / norm(x_pred(7:10)); % 2. 对误差状态做中心差分求雅可比 F zeros(12, 12); eps_j 1e-6; for i 1:12 dxe zeros(12, 1); dxe(i) eps_j; dxm zeros(12, 1); dxm(i) -eps_j; dx_after_e propagate_error(dxe, x, x_pred, f); dx_after_m propagate_error(dxm, x, x_pred, f); F(:, i) (dx_after_e - dx_after_m) / (2*eps_j); end P_pred F * P * F Qd; end function dx_after propagate_error(dx, x, x_pred, f) % 在标称状态x处加误差扰动, 传播后再取回误差状态 xp x; xp(1:6) xp(1:6) dx(1:6); % 位置/速度直接加扰动 dq angle_to_quat(dx(7:9)); % 姿态误差角 - 小角度四元数 xp(7:10) quat_mult(dq, xp(7:10)); xp(11:13) xp(11:13) dx(10:12); xp_pred f(xp); xp_pred(7:10) xp_pred(7:10) / norm(xp_pred(7:10)); dx_after zeros(12, 1); dx_after(1:6) xp_pred(1:6) - x_pred(1:6); dx_after(7:9) small_angle(xp_pred(7:10), x_pred(7:10)); dx_after(10:12) xp_pred(11:13) - x_pred(11:13); end逻辑说明主状态传播只调用一次 rk4_step中心差分里每次循环调用 propagate_error 两次12个维度一共24次动力学积分步长内计算量可接受。函数 angle_to_quat 把三维误差角转成小角度四元数quat_mult 做四元数乘法small_angle 取两个四元数之间的误差角向量这三个自定义函数各自十几行按前面约定的四元数顺序实现即可。参数说明差分步长 eps_j 选1e-6对SI单位制合适如果状态里混入了角度量级差异很大的变量按各分量单独定步长更稳。Qd 是12×12过程噪声协方差初值可按下表量级给然后根据NEES检验微调。误差状态分量Qd对角线初值参考物理含义δr位置(1e-6)^2 * dt未建模轨道摄动δv速度(1e-4)^2 * dt推力误差、摄动加速度δθ姿态误差角(1e-5)^2 * dt重力梯度力矩等未建模δω角速度(1e-6)^2 * dt干扰力矩造成的角加速度3.3 观测更新GPS、星敏感器与陀螺的H矩阵和残差观测更新的核心是把残差投影到误差状态空间。位置观测残差直接是位置差角速度观测残差直接是角速度差姿态观测残差要先算测量四元数与预测四元数的误差四元数再取小角度向量不能直接做四元数减法。function [x_upd, P] mekf_update(x_pred, P, z, obs_type, R) % MEKF更新步: 支持 r(位置) / q(四元数姿态) / w(角速度) 三类观测 % z: 观测向量, R: 对应观测噪声协方差 switch obs_type case r H [eye(3), zeros(3,9)]; y z - x_pred(1:3); case q dq quat_mult(z, quat_conj(x_pred(7:10))); if dq(1) 0, dq -dq; end % 消去符号歧义 y 2 * dq(2:4); % 小角度姿态误差 H [zeros(3,6), eye(3), zeros(3,3)]; case w H [zeros(3,9), eye(3)]; y z - x_pred(11:13); end S H * P * H R; K P * H / S; dx K * y; % 误差状态注入主状态 x_upd x_pred; x_upd(1:6) x_upd(1:6) dx(1:6); x_upd(11:13) x_upd(11:13) dx(10:12); dth dx(7:9); if norm(dth) 1e-12 dq_upd [cos(norm(dth)/2); sin(norm(dth)/2) * dth/norm(dth)]; x_upd(7:10) quat_mult(dq_upd, x_pred(7:10)); x_upd(7:10) x_upd(7:10) / norm(x_upd(7:10)); end P (eye(12) - K * H) * P; P (P P) / 2; % 强制对称 end逻辑说明P 是12×12H 按观测类型取对应行K 维度随之变化。姿态观测残差 y 2*dq(2:4) 用的是小角度近似残差超过5度时建议换成 atan2 形式的精确公式。误差状态注入时姿态误差角要转回四元数再乘到预测四元数上这个顺序不能颠倒。参数说明R 矩阵按传感器手册标定星敏感器一般给 (1e-4)^2 rad^2 量级GPS位置给 1^2 m^2陀螺角速度给 (1e-4)^2 rad^2/s^2。注意星敏给出的姿态通常已经包含安装矩阵校准仿真里做测量生成时也要先转成同一坐标系。3.4 Q和R整定这不是调平滑度的旋钮Q和R整定是EKF使用里最像玄学的部分。原则是R描述传感器噪声的物理精度应该通过标定或手册确定不要乱动Q描述模型误差的大小模型漏了什么就把对应分量Q调大。仿真里常见错误是把Q当“滤波强度”旋钮Q调得越大曲线越平滑但真实误差反而变大。我一般这样起步先用开环真值仿真控制量给定、状态真值已知单独验证EKFQ初值按第3.2节表格给跑一次蒙特卡洛看平均NEES。NEES偏低说明滤波器过自信把Q调大NEES偏高说明过程噪声给大了把Q调小。调参调到NEES落进置信区间后再闭环接控制律此时如果曲线毛刺大优先检查R是不是给小了而不是继续动Q。4. 姿态轨道耦合控制律设计LQR反馈、前馈补偿与Matlab主循环4.1 控制目标与误差状态定义控制目标是让卫星跟踪一条期望轨道和期望姿态。期望轨道通常是机动目标轨迹期望姿态要么对地定向要么保持推力方向在惯性系下固定。误差状态取 δx [δr; δv; δθ; δω]与EKF的误差状态一致控制律直接吃估计器输出的误差。注意位置误差 δr 是在轨道控制坐标系常见取 LVLH 当地轨道坐标系下定义而姿态误差 δθ 是在体坐标系下定义两者混在一个向量里没问题但反馈增益K的坐标一定要对应好。仿真里如果真值和估计值都在惯性系传播控制律内部要先做一次坐标旋转把惯性系误差转到控制坐标系再乘K。4.2 LQR线性化模型C-W方程加刚体姿态方程LQR需要一个线性时不变模型。轨道部分在圆轨道工作点用Clohessy-Wilshire方程姿态部分在零角速度工作点取欧拉方程小扰动。两部分做分块对角合成耦合项当作扰动由反馈裕度吸收这是工程简化做法效果靠闭环仿真验证。function A build_linear_model(p) % 12阶线性模型: 轨道C-W 姿态小扰动 n p.n; % 轨道角速度 rad/s A_cw [0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; 3*n^2 0 0 0 2*n 0; 0 0 0 -2*n 0 0; 0 0 -n^2 0 0 0]; % 工作点取零角速度, 重力梯度力矩当外部扰动 A_att [zeros(3), eye(3); zeros(3), zeros(3)]; A blkdiag(A_cw, A_att); end function K design_lqr(p) % LQR权重: Bryson规则按物理容差给初值 A build_linear_model(p); B [zeros(3,3), zeros(3,3); eye(3)/p.m, zeros(3,3); zeros(3,3), zeros(3,3); zeros(3,3), inv(p.J)]; r_tol 10; % 位置误差容差 10 m v_tol 0.5; % 速度误差容差 0.5 m/s th_tol 0.01; % 姿态误差容差 0.01 rad w_tol 0.01; % 角速度误差容差 0.01 rad/s F_max 20; % 最大推力 20 N M_max 0.2; % 最大控制力矩 0.2 N*m Q_ctrl diag([1/r_tol^2*ones(3,1); 1/v_tol^2*ones(3,1); ... 1/th_tol^2*ones(3,1); 1/w_tol^2*ones(3,1)]); R_ctrl diag([1/F_max^2*ones(3,1); 1/M_max^2*ones(3,1)]); [K, ~, ~] lqr(A, B, Q_ctrl, R_ctrl); end逻辑说明C-W方程前三行是运动学后三行是动力学包含 3n^2 径向项和 2n 科里奥利耦合项这是圆轨道相对运动的固有特性不能删。B矩阵的控制输入是控制坐标系下的推力与力矩与第2章动力学函数的u接口一致。参数说明Bryson规则用物理容差倒数的平方作为权重初值目的是让不同量纲的状态贡献归一化。r_tol、v_tol、th_tol、w_tol 表示你允许的误差幅度F_max、M_max 表示执行机构极限。这一组设置完K矩阵的量级基本合理后续按仿真表现微调。4.3 前馈补偿与执行机构分配耦合控制的关键一笔LQR反馈算出的控制量需要叠加前馈才能处理耦合。推力要跟踪期望加速度同时推力偏置会给姿态一个干扰力矩姿态控制必须预先抵消它。这就是耦合控制里最关键的“前馈补偿”一笔。function [F_cmd, M_cmd] coupled_lqr_correction(err_x, ref, p, K) % err_x: [δr; δv; δθ; δω] 控制坐标系下 % ref : 期望加速度 a_des(3) 等参考量 % 反作用反馈 u_corr -K * err_x; F_corr u_corr(1:3); M_corr u_corr(4:6); % 轨道机动前馈: 期望加速度直接转推力 F_ff p.m * ref.a_des; F_cmd F_corr F_ff; % 推力偏置力矩前馈补偿 M_ff cross(p.d_cm, F_cmd); M_cmd M_corr - M_ff; % 执行机构饱和 F_cmd sat_vec(F_cmd, p.F_max); M_cmd sat_vec(M_cmd, p.M_max); end逻辑说明M_ff d × F 是推力不通过质心产生的干扰力矩姿态控制要抵消它所以从反馈力矩里减去。符号约定要跟第2章动力学里 M_body 一致否则补偿变成叠加姿态发散得更快。饱和函数 sat_vec 按向量模长限制超出上限时等比例收缩而不是逐分量截断避免改变推力方向。参数说明期望加速度 ref.a_des 来自轨道制导层常见做法是跟踪一个机动规划器输出的轨迹。d_cm 在这里同时出现在动力学模型和控制器里如果实际值和模型值不一致控制效果会变差这正好引出第5章的避坑案例。4.4 主仿真循环把估计、控制、积分三个频率串起来dt_sim 0.01; % 动力学积分 100 Hz dt_est 0.1; % EKF估计 10 Hz dt_ctrl 0.1; % 控制 10 Hz T_end 600; t 0; x_true x0; x_est x0; P P0; u_cmd zeros(6,1); for k 1 : T_end / dt_sim % 1. 估计周期 if mod(t, dt_est) 0 z gen_measurement(x_true, R); % 产生观测并加噪声 [x_est, P] mekf_predict(x_est, P, u_cmd, p, dt_est, Qd); [x_est, P] mekf_update(x_est, P, z, r, R_r); [x_est, P] mekf_update(x_est, P, z_q, q, R_q); end % 2. 控制周期 if mod(t, dt_ctrl) 0 err_x get_control_error(x_est, x_ref(t), p); [F_cmd, M_cmd] coupled_lqr_correction(err_x, ref(t), p, K); end % 3. 仿真积分 x_true rk4_step(coupled_dynamics, x_true, [F_cmd; M_cmd], p, dt_sim); t t dt_sim; log_data(t, x_true, x_est, F_cmd, M_cmd); end逻辑说明估计频率和控制频率可以不同这里都取10Hz。EKF预测步的 dt 是0.1s但动力学积分步长是0.01s两者严格保持一致不要用积分步长去做预测。观测生成要在估计时刻进行控制量用的是上一个控制周期的值也就是零阶保持。5. 避坑EKF耦合系统里的四个翻车现场与排查方法5.1 P矩阵失去正定性现象、原因与解法现象仿真跑着跑着P矩阵对角线出现负值卡尔曼增益变成奇怪的数值估计值突然跳变。原因直接状态EKF把四元数当作普通向量更新协方差更新破坏了单位模约束即使做了归一化P矩阵在强观测下也会迅速失去正定性。另一个常见诱因是数值雅可比的差分步长不适合量级差异大的状态。解决切换到误差状态EKF让P只作用于12维误差状态。同时每个观测更新完做一次 P (PP)/2 强制对称必要时加微小对角项保证正定。如果不想换框架至少要在预测和更新后都归一化四元数但这不是长久之计。5.2 姿态残差为0但滤波仍然漂移隐藏的四元数符号歧义现象q 观测残差一直是零向量EKF的姿态估计误差却持续增长位置和角速度都很正常。原因四元数 q 和 -q 表示同一个姿态。如果测量生成时四元数符号翻转姿态残差计算里 dq 接近 180° 旋转小角度近似失效或者说残差被错误地算成接近零。Matlab自带工具包某些函数对四元数符号敏感很多人在这里浪费过一整天。解决在姿态残差计算里判断 dq(1)小于零就取负强制把符号唯一的表示挑出来。代码已在第3.3节给出。另一个辅助手段是观测更新前先打印 dq 的四元数模长和第一个分量人工确认残差量级在预期范围内。5.3 推力偏置力矩正反馈补偿符号搞反的经典翻车现象闭环仿真里轨道控制收敛姿态误差越来越大检查控制器逻辑没发现明显bug。原因前馈补偿 M_ff 的符号和动力学模型里 M_body 的符号约定不一致。动力学里是 M_body M_cmd d×F控制器里如果写成 M_cmd M_corr M_ff干扰力矩变成两倍而不是抵消姿态必然发散。解决把动力学和控制器里的力矩符号写在同一个注释块里仿真前做一个纯开环检查常值推力下控制器只给 M_corr 0M_cmd -d×F看姿态是否保持稳定。这一步验证通过再闭环符号问题基本绝迹。5.4 大偏心率轨道固定步长积分发散预测步的按需做现象圆形轨道仿真一切正常换成椭圆轨道后EKF速度估计发散真值轨迹也出现可疑振荡。原因RK4固定步长在近地点速度大、状态变化快时截断误差急剧增大EKF预测步跟着把错误当成过程噪声吸收。另一个隐藏问题是数值雅可比用固定 eps_j1e-6对近地点位置量级变化很大的状态不够自适应。解决真值生成用 ode45 变步长积分EKF预测步继续用固定步长但缩短到近地点步长比如0.01s或者对差分步长按分量自适应取 eps_j 1e-6 * max(1, abs(x(i)))。不要用同一套固定步长硬跑椭圆轨道这是血泪经验。6. 用蒙特卡洛与NEES检验滤波器的可信度验证技巧与进阶做法EKF调完参数后判断滤波器是否可信的标准不是姿态曲线漂不漂亮而是NEESNormalized Estimation Error Squared。NEES的定义是 x_err^T P^{-1} x_err / nx其中 x_err 是估计误差P 是滤波器输出的协方差nx 是误差状态维数12。对每一次蒙特卡洛仿真逐时刻计算NEES再对多次仿真的同一时刻做平均得到平均NEES。function nees_mean calc_nees(x_err_all, P_all) % x_err_all: nx-by-N-by-M 误差状态序列 % P_all: nx-by-nx-by-N 协方差序列(每次蒙特卡洛相同) % 平均NEES逐时刻计算 [nx, N, M] size(x_err_all); nees_t zeros(N, M); for m 1:M for k 1:N xk x_err_all(:,k,m); nees_t(k,m) xk / P_all(:,:,k) * xk / nx; end end nees_mean mean(nees_t, 2); end判定标准12维误差状态95%置信区间对应的平均NEES约在0.72到1.32之间。平均NEES远大于1.3说明滤波器过于自信P给得太小估计误差比滤波器声称的大得多平均NEES远小于0.7说明保守P偏大滤波器没有充分利用观测信息。这里不能只盯曲线否则滤波器的过自信会在真实任务里变成灾难。验证流程我一般分成三步。第一步开环验证控制量固定或按脚本给定状态真值来自仿真EKF只吃观测跑20次蒙特卡洛看平均NEES同时观察位置和姿态误差的均值是否在零附近。第二步注入参数失配把模型里的 d_cm 偏移2mm、J 矩阵加5%误差再看NEES是否仍在区间内。失配后NEES轻微上升是正常的如果急剧超限说明Q给得太小模型误差没被编码进去。第三步闭环保验证接上第4章LQR跑完整轨道机动场景检查轨道重访误差和姿态残差是否收敛在容差内。一个提升效率的技巧把 Qd 和 R 的整定写成脚本函数蒙特卡洛结果自动返回平均NEES的均值与区间外比例用表格对比不同参数组。这样调参从“看曲线猜”变成“量化决策”。我做这类系统时的习惯是先把开环NEES调到带内再闭环保验证不然两个环互相怪罪翻车概率高得多。希望帮到你。本文还有配套的精品资源点击获取