1. 三自由度跑得通不代表六自由度能用我见过太多次同样的场景有人把纵向三自由度模型调得漂漂亮亮高度阶跃响应干净利落然后信心满满地把同一套控制参数搬到六自由度模型上仿真跑了 0.8 秒就发散了。第一反应都是控制器不行于是换滑模、换自抗扰、加观测器折腾两周问题依旧。最后打开日志一看发散的原因是滚转通道里一个耦合项把侧滑角顶到了 40 度以上而控制器的舵面早就饱和了。六自由度高超声速飞行器的建模与控制器设计核心难点从来不是写出十二个微分方程这件事本身而在于高超声速这个工况把所有在常规飞行器上可以忽略的耦合项、非线性项、参数时变项同时放大了一到两个数量级。你写的是同一套刚体动力学但 6 马赫和 0.6 马赫下模型表现出来的脾气完全是两回事。这篇文章我想按一个真实的工程流程来组织先把六自由度相对三自由度多出来的东西讲清楚再一步步把模型搭到能跑、能配平、能线性化然后才是控制器结构的选择和落地。整篇内容适合两类人看一类是做控制方向、准备拿高超声速飞行器当验证平台的另一类是做过数学建模、想从能列方程进阶到能跑出一条可信的仿真曲线的。中间的代码和检查清单可以直接抄去用我踩过的坑会明确标出来。2. 多出来的那三个自由度到底在系统里干了什么2.1 从质点弹道到刚体姿态状态量多了什么三自由度模型通常指的是质点弹道模型状态量是位置三个分量加速度三个分量总共六个状态实际上描述的是一个质点在空间中的运动。气动力只算升力和阻力方向由速度矢量决定飞行器的姿态信息完全被揉进了一个迎角假设里——你默认迎角是配平值或者干脆当成一个可以直接指令的输入。六自由度模型不一样它把绕质心的转动单独拿出来建模。典型的状态量组织是十二个速度相关三个速度大小 $V$、航迹倾角 $\gamma$、航迹偏角 $\chi$姿态相关三个迎角 $\alpha$、侧滑角 $\beta$、倾侧角 $\mu$或者直接用欧拉角 $\phi,\theta,\psi$角速率三个滚转率 $p$、俯仰率 $q$、偏航率 $r$位置相关三个地面系下的 $x_E, y_E, z_E$多出来的六个状态本质上是回答了两个三自由度回答不了的问题飞行器现在朝哪儿它转得多快这两个问题一旦被打开惯量张量、气动力矩、舵面力矩、陀螺耦合$p\times r$、$q\times r$ 那些交叉项就全进来了。我第一次做这个转换时犯的错很典型把三自由度里的迎角指令直接接到六自由度里当成状态量用结果控制器实际控制的是 $q$而 $\alpha$ 是通过 $\dot\alpha q - \dot\gamma$ 间接被影响的。这两者之间差了一个积分环节和一个重力投影项相位差直接从正变成负 90 度稳定裕度全没了。2.2 高超声速条件下被放大的耦合项有几个耦合项在低速时几乎没人管到高超声速会变成主要矛盾。第一是惯性耦合。飞行器通常关于纵向对称面对称所以 $I_{xy}I_{yz}0$但 $I_{xz}$ 一般不为零——尤其是细长体构型主轴和体轴之间有夹角。转动方程里 $I_{xz}$ 会让滚转和偏航互相串门你打副翼想滚转结果带出了偏航力矩。低速时这个耦合量级小姿态回路带宽低能压住高超声速时动压大、气动舵效强、姿态回路带宽被推高这个耦合就压不住了。第二是气动-推进耦合。吸气式构型的前体本身就是压缩面迎角一变进气道入口的流场就变推力跟着变推力方向又跟机体轴对齐推力的俯仰分量随着迎角变化。这是个闭环的正反馈路径建模时必须把推力写成 $T f(Ma, \alpha, \varphi)$ 的形式$\varphi$ 是燃料当量比而不是一个常数或者只跟油门开度有关。第三是运动学奇异。在 $\gamma$ 接近 ±90° 或者 $\beta$ 接近 ±90° 时风轴系下的运动学方程会出现 $\tan\beta$ 和 $1/\cos\beta$ 的奇点。三自由度模型里因为你根本不跟踪 $\beta$这个奇点是隐藏的六自由度一上来就暴露了。注意如果你的仿真在某个瞬间突然报出 $10^{15}$ 量级的状态量先别怀疑控制器去查一下 $\cos\beta$ 是不是过零了。2.3 什么阶段可以退回三自由度也不是所有工作都非得六自由度。我自己的判断标准是这样的工作阶段建议模型理由总体方案论证、弹道优化三自由度只关心能量和射程姿态是内环的事纵向控制器初版设计三自由度 等效舵效迭代快调参直观通道耦合分析、协调控制六自由度不引入耦合就看不到问题大包线鲁棒性验证六自由度 参数拉偏交叉耦合是鲁棒性的主要杀伤项弹性/热效应影响评估六自由度 弹性模态需要完整的力矩分解一句话概括验证控制律的时候必须上六自由度设计控制律初值的时候可以用三自由度省时间。但这个省时间的前提是你清楚自己在做简化而不是把三自由度的结论当成终局。3. 把方程搭起来坐标系、状态量和力的分解3.1 先钉死符号约定再动笔写方程这是我最想强调的一条经验六自由度建模之所以让人头大八成不是因为数学难而是因为符号约定没统一。同一套风轴系方程不同教材在三个地方会打架倾侧角 $\mu$ 以左倾还是右倾为正航迹偏角 $\chi$ 的零度方向是北还是东侧滑角 $\beta$ 正号对应速度矢量偏向机头左侧还是右侧。这三处只要有一处写反仿真里的表现都是配平能收敛但飞不稳排查起来极其痛苦。我的做法是在代码最前面写一段约定注释并且用退化检验锁死它。具体来说把六自由度方程令 $\beta0$、$\mu0$、$\chi$ 常数它必须严格退化成纵向方程组$$ \dot V \frac{T\cos\alpha - D}{m} - g\sin\gamma $$$$ V\dot\gamma \frac{T\sin\alpha L}{m} - g\cos\gamma $$$$ \dot\alpha q - \dot\gamma $$如果退化后对不上说明符号或者某一项漏了先别往下走。这一步花十分钟能省掉后面两天的调试。3.2 十二个状态量的动力学方程怎么分组我习惯把方程分成三组来写代码里也分成三个子函数这样调试时可以单独冻结某一组。第一组质心平动的力方程风轴系$$ m\dot V T\cos\alpha\cos\beta - D - mg\sin\gamma $$$$ mV\dot\gamma T(\sin\alpha\cos\mu \cos\alpha\sin\beta\sin\mu) L\cos\mu - Y\sin\mu - mg\cos\gamma $$$$ mV\cos\gamma,\dot\chi T(\sin\alpha\sin\mu - \cos\alpha\sin\beta\cos\mu) L\sin\mu Y\cos\mu $$第二组绕质心的力矩方程体轴系$$ \dot p \frac{I_{zz}\bar L I_{xz}\bar N - I_{xz}(I_{xx}-I_{yy}I_{zz})pq (I_{xz}^2 I_{yy}I_{zz} - I_{zz}^2)qr}{I_{xx}I_{zz} - I_{xz}^2} $$$$ \dot q \frac{\bar M - I_{xz}(p^2-r^2) - (I_{xx}-I_{zz})pr}{I_{yy}} $$$$ \dot r \frac{I_{xx}\bar N I_{xz}\bar L (I_{xx}^2 - I_{xx}I_{yy} I_{xz}^2)pq - I_{xz}(I_{xx}-I_{yy}I_{zz})qr}{I_{xx}I_{zz}-I_{xz}^2} $$第三组姿态运动学风轴系$$ \dot\alpha q - (p\cos\alpha r\sin\alpha)\tan\beta \frac{mg\cos\gamma\cos\mu\cos\alpha - mg\sin\gamma\sin\alpha - L}{mV\cos\beta} $$$$ \dot\beta p\sin\alpha - r\cos\alpha \frac{Y T\cos\alpha\sin\beta}{mV} \frac{g}{V}(\cos\alpha\sin\gamma - \sin\alpha\cos\gamma\cos\mu) $$$$ \dot\mu \frac{p\cos\alpha r\sin\alpha}{\cos\beta} \frac{g}{V}\left(\cos\gamma\cos\mu\tan\beta - \sin\gamma\right)\cdot(\cdots) $$第三组的最后一项各家写法差异最大我这里不展开写全因为不同文献对 $\mu$ 的定义确实不一致。判断自己写对没有最有效的办法是做一次零舵面自由响应检验给定一个初始侧滑角和初始滚转角速率飞行器应该在无控条件下做一次单调收敛的荷兰滚振荡而不是发散或者不振荡。提示如果你不打算做复杂的机动飞行只在配平点附近做小扰动任务可以先用体轴系欧拉角形式$\phi,\theta,\psi$ $p,q,r$避免 $\tan\beta$ 和 $1/\cos\beta$ 的奇异性。代价是姿态角在俯仰接近 ±90° 时也会奇异但对高超声速巡航这种小姿态变化的任务完全够用。3.3 气动、大气、推进三个必须外挂的模型动力学方程只是骨架真正决定仿真像不像的是外挂的三个模型。气动系数一般以表格形式给出维度是马赫数、迎角、侧滑角、舵面偏角。表格的插值方式很重要下面第 6 节会专门讲。系数形式通常是$$ C_L C_{L0}(Ma) C_{L\alpha}(Ma)\alpha C_{L\delta_e}(Ma)\delta_e $$我一般会保留到二次项甚至三次项因为高超声速大迎角下 $C_L$ 对 $\alpha$ 明显是非线性的线性外推会把配平点算偏。大气模型在 20~40 km 高度段用 1976 标准大气的分段线性温度模型就够了密度用理想气体状态方程推function [rho, a] atmos(h) % h 为几何高度(m)适用 0~40 km 段简化实现 R 287.05; g0 9.80665; p0 101325; h min(max(h, 0), 40000); if h 11000 T 288.15 - 0.0065*h; p p0 * (T/288.15)^(g0/(R*0.0065)); elseif h 20000 T 216.65; T11 288.15 - 0.0065*11000; p11 p0 * (T11/288.15)^(g0/(R*0.0065)); p p11 * exp(-g0*(h-11000)/(R*T)); else T 216.65 0.001*(h-20000); T20 216.65; p20 5474.9; p p20 * (T/T20)^(-g0/(R*0.001)); end rho p / (R*T); a sqrt(1.4 * R * T); end推进模型是吸气式构型里最难搞的部分。我一般用一个简化参数化模型function T engine(Ma, alpha, phi, cfg) % phi: 燃料当量比 0~1 if Ma cfg.Ma_ign || Ma cfg.Ma_max T 0; return; end T cfg.T_ref * (Ma/cfg.Ma_ref)^0.8 ... * (1 - cfg.k_alpha*(alpha - cfg.alpha_ref)^2) ... * phi; end这个式子没有任何物理严谨性但它抓住了三件事推力随马赫数增长、迎角偏离设计点后进气道性能下降、推力与当量比近似线性。做控制律验证足够了。真要用到定量结论得换成基于热力循环的计算模型。3.4 弹性模态加不加加几个细长体加上气动加热导致的结构刚度下降一阶弯曲频率可能落到姿态回路带宽的 2~3 倍以内这时候不加弹性模态仿真会给出过于乐观的稳定裕度。我的经验是如果姿态回路带宽超过一阶弯曲频率的 1/3就必须加。加的时候取前两阶弯曲模态一阶垂直、一阶扭转就够模态阻尼取 0.02~0.05频率做 ±20% 拉偏。加到方程里的方式是把模态广义坐标当成额外状态气动力和力矩分别对模态坐标求偏导形成耦合项。代价是状态量从 12 个涨到 16~18 个而且引入了两个高频极点。这时候定步长积分器要小心步长得压到 1 ms 以下才能把 10 Hz 以上的模态基频积准。4. 让模型跑起来配平、线性化、自检4.1 状态向量组织和求解器选择我习惯把十二个状态排成[V; gamma; chi; alpha; beta; mu; p; q; r; xE; yE; zE]高度用h -zE反推。这样排列有个好处前九个状态是飞机本体的运动后三个只跟导航有关做姿态回路分析时可以直接把后三个截掉。求解器选择上我基本不用ode45。原因是气动系数用的是表格插值插值节点的导数不连续ode45的变步长策略会反复在节点附近缩小步长仿真时间被拖得极长而且结果还带着插值噪声。改用定步长四阶龙格库塔ode4步长 0.002 s稳定得多。如果加了弹性模态步长取 0.0005 s。function dx hgv6dof(~, x, cfg) V x(1); gam x(2); chi x(3); alp x(4); bet x(5); mu x(6); p x(7); q x(8); r x(9); h -x(12); % 姿态限幅防止在配平搜索过程中 cos(beta) 越界 bet max(min(bet, 1.4), -1.4); a atmos(h); Ma V / a; qbar 0.5 * atmos(h) * V^2; de cfg.u(1); da cfg.u(2); dr cfg.u(3); phi cfg.u(4); [CX, CY, CZ, Cl, Cm, Cn] aero_lookup(Ma, alp, bet, de, da, dr); T engine(Ma, alp, phi, cfg); L -qbar * cfg.S * CZ; D -qbar * cfg.S * CX; Y qbar * cfg.S * CY; Mx qbar * cfg.S * cfg.b * Cl; My qbar * cfg.S * cfg.c * Cm; Mz qbar * cfg.S * cfg.b * Cn; % ... 力方程、力矩方程、运动学方程按 3.2 节写 ... end有个细节值得单独说限幅一定要写成max(min(...))而不是if分支。因为配平搜索和数值线性化会对状态量做有限差分扰动if分支在边界上会造成函数值跳变雅可比矩阵直接算废。4.2 配平是所有控制器设计的前置动作没有配平点后面所有事都做不了线性化没法做增益调度没有网格点控制器连初值都给不出来。配平的本质是解一个非线性方程组找到一组状态和控制输入让所有状态导数同时为零。工程上我会把问题降维——固定 $V$ 和 $h$ 作为包线网格点把 $\gamma$、$\chi$、$\beta$、$p$、$q$、$r$ 全部固定为零只让 $\alpha$、$\delta_e$、$\delta_a$、$\delta_r$、$\varphi$ 作为未知量然后要求$$ \dot V 0,\quad \dot\gamma 0,\quad \dot\alpha 0,\quad \dot q 0 $$四个方程、五个未知量欠定。解决办法是再加一个约束比如固定倾侧角为零或者固定当量比把未知量降到四个。function F trim_residual(u, V0, h0, cfg) cfg.u [u(1); u(2); u(3); u(4)]; % de, da, dr, phi x0 zeros(12,1); x0(1) V0; x0(12) -h0; x0(4) u(5); % alpha dx hgv6dof(0, x0, cfg); F [dx(1); dx(2); dx(4); dx(8)]; end u0 [deg2rad(2); 0; 0; 0.55; deg2rad(3)]; opts optimoptions(fsolve,TolFun,1e-10,TolX,1e-10, ... StepTolerance,1e-10); u_trim fsolve((u) trim_residual(u, V0, h0, cfg), u0, opts);这里有个坑我必须提fsolve默认的有限差分步长对舵面这个量级来说太大了。舵面配平值通常在 1~5 度也就是 0.02~0.09 rad如果差分步长取默认的 $\sqrt{\epsilon}\approx 1.5\times10^{-8}$那没问题但如果你的残差函数里存在插值微小扰动可能落在同一个插值格子里梯度算出来是零fsolve直接原地踏步。解决办法是把舵面的初始猜测放在插值网格的格子内部别贴着节点。4.3 小扰动线性化与状态空间提取配平点拿到之后线性化有两条路解析求导和数值差分。我强烈推荐数值中心差分因为你的气动模型是查表加插值解析求导根本不现实。function [A, B] num_jac(f, x0, u0, nx, nu) h_x 1e-6; h_u 1e-6; f0 f(x0, u0); A zeros(nx, nx); B zeros(nx, nu); for i 1:nx dx zeros(nx,1); dx(i) h_x; A(:,i) (f(x0dx, u0) - f(x0-dx, u0)) / (2*h_x); end for j 1:nu du zeros(nu,1); du(j) h_u; B(:,j) (f(x0, u0du) - f(x0, u0-du)) / (2*h_u); end end步长的选择有讲究。取太小会被浮点误差淹没取太大又会被插值的局部非线性污染。我的经验值是把步长设成该状态量典型变化范围的 $10^{-6}$ 倍迎角典型变化 0.1 rad步长就取 $10^{-7}$速度典型变化 1000 m/s步长取 $10^{-3}$。统一用一个绝对步长在十二个量级差异巨大的状态上是不行的。线性化之后一定要做特征值审查。高超声速飞行器在巡航点的纵向通常是一对长周期复极点 一对短周期复极点横航向是滚转收敛 荷兰滚 螺旋。如果你算出来的特征值里出现了明显的正实部先别急着说这个点本来就不稳定——大部分时候是配平没收敛干净或者线性化步长取错了。4.4 我自己在用的模型自检清单模型写完别急着接控制器。按下面这套顺序过一遍每一条都能挡住一类典型错误检查项操作方法通过判据退化一致性令 $\beta\mu0$对比三自由度结果误差 1%力的平衡配平后手算 $L$、$D$、$T$ 的合分量与 $mg$ 平衡到 3 位有效数字力矩平衡配平后检查 $\bar L, \bar M, \bar N$与惯性项之和 1e-6 rad/s²静稳定性固定舵面给 ±1° 迎角扰动俯仰力矩回中荷兰滚衰减给 3° 初始侧滑无控响应振荡收敛周期合理能量守恒关闭推力和阻力纯抛体总机械能恒定到 1e-8惯性矩阵正定特征值检查最小特征值 0能量守恒这一条特别好用。把气动力和推力全部置零飞行器就是一个自由落体$V^2/2 gh$ 应该是常数。如果你的方程里某个投影项写错了这个量就不会守恒误差会随时间线性漂移。很多隐蔽的坐标变换错误都是这么被查出来的。5. 控制器结构怎么选从增益调度到非线性5.1 单点 LQR 为什么撑不起全包线有人第一次做这个题会在 8 马赫、30 km 的配平点设计一个 LQR然后拿它去跑 6~10 马赫、20~40 km 的全包线扫掠。结果必然是在包线边缘发散。原因很直白动压在这段包线里变化接近一个数量级舵效跟着动压走建模增益也就跟着变。同一个 $\delta_e$ 在 20 km 产生的俯仰力矩可能只有 40 km 的三倍以上。固定的状态反馈增益不可能同时匹配两端。更重要的是构型本身在全包线内就不是同一架飞机。迎角配平值随马赫数和高度变化导致前体压缩流场变化气动焦点位置移动俯仰静稳定裕度可能从小正变成了小负。这种情况下单点 LQR 的相位裕度估计完全失真。5.2 内外环分离 增益调度工程上最稳的路我推荐的结构永远是时标分离把状态按响应速度分成三层每层用不同带宽的回路控制。快回路内环$p, q, r$ 角速率。带宽最高用动态逆或 LQR。因为这层动力学近似线性增益调度网格可以粗一点。中回路$\alpha, \beta, \mu$ 姿态角。带宽取内环的 1/5~1/3输出是角速率指令。慢回路外环$V, h$或者 $V, \gamma$。带宽取中环的 1/5输出是迎角/倾侧角指令。这样分层的最大好处是每一层的增益调度表都是低维的。比如内环只需要按马赫数和动压做二维调度其实按动压一维就行中环按马赫数和高度做二维外环按高度做一维。总共三张表比一个十二维增益调度好维护一万倍。% 内环增益调度按动压一维分段 qbar_grid [5e3, 1e4, 2e4, 4e4, 8e4]; Kq_tab zeros(1, numel(qbar_grid)); for k 1:numel(qbar_grid) [A, B] linearize_at(qbar_grid(k), cfg); Kq_tab(k) lqr(A, B, Qq, Rq); end Kq_now interp1(qbar_grid, Kq_tab, qbar_now, linear, extrap);这里有个实际会遇到的问题LQR 算出来的增益矩阵是多维的interp1只能插一维标量。解决办法是逐元素插值把 $K$ 矩阵拆成 $n\times m$ 个标量各自插值再拼回去。我知道这听起来有点笨但工程上就是这么做稳定可靠。提示调度量一定选动压而不是高度。动压综合了高度和速度的影响用它做调度量增益表的非线性程度低得多网格可以更稀。5.3 动态逆和反步法拿模型精度换性能如果模型精度足够动态逆NDI能在整个包线内给出几乎一致的闭环响应。它的思路很直接把非线性项全部搬到控制律里抵消掉剩下一个线性系统用线性方法设计。以俯仰通道为例角速率方程是$$ \dot q \frac{\bar M}{I_{yy}} f_{nl}(p, r) $$控制律取成$$ \delta_e \frac{I_{yy}}{\bar M_{\delta_e}}\left(\nu - \frac{\bar M_{other}}{I_{yy}} - f_{nl}\right) $$$\nu$ 是新定义的线性控制量可以随便用一个 PI 或者 LQR 来生成。这样闭环就变成 $\dot q \nu$与工作点无关。反步法Backstepping思路类似但是递归地设计虚拟控制不需要精确求逆对模型误差的容忍度略高。代价是每一步都要对虚拟控制求导如果你的模型是查表形式这个导数要解析算出来非常麻烦一般用数值微分代替会引入噪声。我踩过的坑是动态逆对建模误差极其敏感尤其是气动系数和惯量的误差。我给一个 10% 的 $C_{m\delta_e}$ 误差闭环带宽就掉了 8%相位裕度从 50 度掉到 30 度。所以如果用动态逆务必在最终验证时做参数拉偏。5.4 滑模、自抗扰、自适应什么时候值得上这几种方法我都试过说点实话。滑模SMC鲁棒性确实好但抖振是个绕不开的问题。用饱和函数替代符号函数能压住抖振代价是牺牲了一部分鲁棒性用边界层闭环精度跟边界层厚度正相关带宽又跟噪声耦合。我的经验是合理取值下SMC 相比增益调度 LQR 的优势主要体现在大扰动场景正常巡航段优势不明显。自抗扰ADRC核心是扩张状态观测器ESO把所有未建模动态和扰动打包估计出来再补偿对模型精度要求最低。问题在 ESO 的带宽观测器带宽通常要取闭环带宽的 3~5 倍才有估计精度而高超声速的姿态回路带宽已经不低再加 5 倍观测器的噪声放大就很可观了。如果你的角速率信号是带 1e-3 rad/s 量级噪声的ESO 出来的补偿量会很跳。自适应L1 等理论保证漂亮但工程落地时参数多。自适应增益设小了收敛慢设大了在饱和边界容易振荡。我一般只在有明确参数不确定性比如气动系数随烧蚀变化的场景用它。下面这张表是我这几年攒下来的选型参考方法模型精度要求计算量全包线适应性调试耗时我实际会用的场景增益调度 LQR中极低好低首选基线方案动态逆高低极好中模型可信、包线跨度大反步法高中好中高有明确级联结构时滑模低中好中存在大扰动、不确定项ADRC极低中高好中模型粗糙、信号干净自适应低高好高参数慢时变明显我的实践结论是先用增益调度 LQR 搭一个能跑通全包线的基线把整个仿真链路和验证流程都跑顺然后再考虑上非线性方法做对比。跳过基线直接上高级算法你连是算法不行还是模型错了都分不清。6. 仿真现场最常见的几类翻车6.1 积分发散八成不是控制器的锅这是最高频的问题。我的排查顺序是固定的看第一个发散的状态量是哪个。如果是角速率查力矩方程如果是迎角或侧滑角查运动学方程如果是速度查力方程。看发散的时间尺度。1 ms 内爆掉通常是代数环或者除零0.1 s 内爆掉通常是符号错误或者增益符号反了几秒后慢慢漂通常是配平没收敛或者有积分漂移。把控制器换成开环配平值再跑一次。如果开环也发散那肯定是模型问题跟控制器无关。这一步能省掉大量无谓的调参。把所有气动系数乘以 0.9 再跑。如果原本发散的情况变好了说明是某个通道增益过高、相位裕度不足如果更差了说明问题在别处。我印象最深的一次仿真在 3.2 s 突然爆掉查了半天控制器最后发现是大气模型里高度插值在 32000 m 处有个台阶——因为我在那个高度上换了一段温度梯度公式两段之间温度跳了 0.6 K密度跳了一点动压跟着跳舵效突变环路瞬间失稳。大气模型的分段一定要保证连续最好用pchip而不是interp1(...,linear)。6.2 气动插值不连续造成的假抖振气动系数表格通常是按马赫数、迎角、舵偏角三维给出的。如果用interp1的线性插值函数本身连续但一阶导数不连续用interp3的默认线性插值也一样。这在定步长积分里表现得特别明显舵面偏转一点点气动力矩的斜率突变角加速度出现尖峰控制器以为是真实扰动反打一下结果形成一个自激的高频振荡。你会以为是控制器的抖振实际上问题在插值上。解决办法有两个一是把interp1/interp3换成保形三次插值MATLAB 里用pchip或者splineinterp3可以指定spline或cubic二是干脆把表格在离线阶段拟合成多项式或者解析函数运行时不算插值。我一般选第二个因为拟合之后模型可微做线性化和动态逆都方便代价是极端工况下拟合精度会掉一点。6.3 步长、速率限幅和执行机构建模很多人做控制器设计时默认舵面可以瞬间从 0 打到 20 度。真实舵机做不到一般速率限幅在 30~60 度/秒带宽 10~20 Hz。你设计的控制律如果要求的舵面速率超过了这个值仿真里会出现极限环而你完全不知道为什么。把执行机构的一阶惯性 速率限幅加进仿真是让控制器从仿真好看变成工程可用的关键一步。我的做法很简单function [de_out, de_rate] actuator(de_cmd, de_prev, dt, cfg) % 一阶惯性 速率限幅 tau 1/(2*pi*cfg.act_bw); de_tmp de_prev (de_cmd - de_prev) * dt / (tau dt); rate (de_tmp - de_prev) / dt; rate max(min(rate, cfg.rate_max), -cfg.rate_max); de_out de_prev rate * dt; de_out max(min(de_out, cfg.de_max), -cfg.de_max); end加上这个之后你会发现原本响应很快很爽的控制器变得温吞了需要重新设计带宽。这就是真实的样子。6.4 配平点漂移与控制权限饱和增益调度的前提是每个包线点都能配平到一组合理的状态和舵面值。但实际扫包线时会发现某些点上舵面配平值已经打到了机械限幅的 80% 以上这时候再给它任何指令都没有剩余控制权限了。我的处理原则是在包线网格扫描阶段就把这些点标出来控制权限不足的区域直接裁掉不要在不可行的工作点上设计控制律。具体判据是任一舵面配平值超过机械限幅的 60%或者迎角配平值超过配平允许上限的 70%这个点就不进入增益调度表。另外倾侧角 $\mu$ 的控制权限通常是最紧张的。高超声速构型一般靠副翼和方向舵差动产生滚转力矩$\mu$ 的可用范围可能只有 ±30 度。如果你的任务要求大倾侧机动得提前在包线评估里算清楚。6.5 角度制、弧度制和单位制这条听起来最基础但栽在上面的人最多。我自己的规则是内部计算全部用弧度只有在日志输出和界面显示的时候才转成度。一旦在代码里混着用sin(3.14)和sin(3.14°)差三个数量级仿真直接爆掉而且爆的方式和控制器发散长得一模一样。还有一个隐蔽的惯量矩阵的单位。有些气动数据手册里 $I_{yy}$ 给的是 kg·m²有些给的是 slug·ft²中间差 1.36 倍。拿错一次俯仰响应频率差 16%你看曲线是看不出来的。注意在hgv6dof函数的入口和出口各加一个 assert检查所有角度量是否在 [-π, π] 之内。花两行代码能挡掉很多问题。7. 把模型和控制器串起来之后我会怎么验收模型能跑、控制器能稳住只是及格线。真正决定这套东西能不能拿出去用的是下面这套验收流程。我一般分四轮每轮加一个维度的不确定性。第一轮是包线网格扫掠。马赫数从 5 到 10 每隔 1 取一个点高度从 20 km 到 40 km 每隔 5 km 取一个点一共 30 个工作点。每个点做一次配平、线性化、闭环特征值检查。判据是所有极点实部小于 -0.5阻尼比大于 0.3舵面配平值不超过限幅的 60%。这一轮下来通常能砍掉 20%~30% 的包线点剩下的才是真正可用的工作区域。第二轮是参数拉偏。气动系数整体 ±15%单个舵效系数 ±20%惯量 ±10%大气密度 ±10%。我用拉丁超立方采样抽 200 组每组跑一次全包线的阶跃响应统计超调量和调节时间的分布。如果 95% 分位数下的超调超过 20%说明控制律的鲁棒性不够得回头调增益或者换方法。第三轮是时域大扰动。给 5 度的迎角阶跃、3 度的侧滑脉冲、200 m/s 的速度扰动看恢复过程。这一轮主要暴露饱和问题——很多在频域里看着很稳的设计在大扰动下会因为舵面饱和产生退绕响应变得很慢甚至出现二次超调。处理办法是加抗积分饱和或者用指令限幅。第四轮是带执行机构和不连续性的联合仿真。把定步长缩到 0.0005 s加上舵机模型、气动插值、传感器噪声角速率加 1e-3 rad/s 的高斯白噪声跑 200 秒的长时仿真看有没有累积漂移或者极限环。这一轮最容易发现的问题是窄带滤波器引入的相位滞后和执行机构带宽不足之间的耦合振荡。四轮下来如果都过了这套模型和控制器才算是能交付的状态。我一般会把每一轮的配置、随机种子、统计结果都存成.mat因为后面写报告和复发问题时一定会用到重跑一遍的成本太高了。最后再分享一个我觉得挺有用的小技巧把配平、线性化、控制器综合这三步全部脚本化并让脚本在包线网格上自动跑。人工点一个点做一次的时代已经过去了30 个工作点手动做一遍至少要一天脚本化之后十几分钟就跑完而且中间任何一次参数修改都能一键复现。这套脚本我写了大概 400 行比控制器本身长得多但它省下来的时间早就超过写它的成本了。
