飞翼无人机鲁棒控制设计：H∞、滑模与自适应方法及Matlab实现

简介面向航空航天工程与无人机控制研究者的Matlab飞翼无人机鲁棒控制实现资源以H∞控制、滑模控制、自适应控制为主线针对飞翼布局下气动不确定性与模型摄动所导致的飞行控制难题给出从系统建模、不确定性分析到控制器设计、仿真验证的完整链路。压缩包共13个文件包括10个m脚本源码、2张开闭环响应对比图像和1份图文并茂的学术文档包体仅174KB便于快速下载与本地复现。已有154人学习下载适合具备自动控制理论基础并希望深入研究无人机鲁棒控制的科研人员与工程师在校学生也可作为案例教学材料使用。内容充分体现‘最坏情况设计’思想其中频域整形与控制器综合脚本可直接复用相关模块有助于开展非线性仿真与参数摄动试验读者可据此构建自己的验证流程也可将脚本迁移至工程项目初步设计参考。1. 为什么飞翼无人机的鲁棒控制不是调参能解决的飞翼无人机取消了水平尾翼和垂直尾翼整架飞机就是一个升力面。气动中心与重心的相对位置几乎决定了纵向稳定性而重心会随燃油消耗、载荷抛放明显移动。巡航点调试好的PID控制器在重心后移3%MAC后相位裕度可能从45°掉到12°再叠加阵风和舵机延迟迎角就开始振荡。这不是调参能解决的问题而是控制器本身没有把不确定性纳入设计。鲁棒控制把模型误差当作对手在设计阶段就要求控制器在给定的参数摄动、未建模动态内都保持稳定。最简单的说法是为最坏情况做设计。这套Matlab工程从飞翼构型的线性化模型出发先用makeweight构造加权函数再用hinfstruct做H∞综合同时给出滑模和自适应控制器的实现最后用Monte Carlo仿真量化不确定性影响。对做飞控算法预研的工程师、上过自动控制原理课又想在无人机上落地的研究生都是一个可复现的参考框架。2. 飞翼无人机建模与不确定性分析从状态方程到加权函数2.1 飞翼构型的纵向线性化状态空间模型飞翼的气动布局决定了它横航向和纵向耦合强但控制设计第一步仍然是从纵向通道开始。在巡航点做小扰动线性化状态向量一般取x [ΔV, Δα, Δq, Δθ]^T输入是升降舵偏角Δδe和推力ΔδT。状态空间形式是Δx_dot A Δx B Δu Δy C Δx D ΔuA矩阵里最关键的是Cm_α俯仰力矩对迎角的导数。常规飞机Cm_α为负值是稳定的飞翼为了巡航效率常把静稳定度设计得很低甚至放宽静稳定性导致Cm_α接近零或者为正这就是开环不稳定的根源。fland.m在资源里承担模型定义的工作。常见写法是传入飞行状态和机体参数返回A、B矩阵和配平点function [A, B, C, D, trim] fland(V0, alpha0, S, m, Iy, cbar, rho) % 气动导数这里用一组飞翼常用量级作示意 CL_alpha 4.9; Cm_alpha 0.03; % 注意飞翼Cm_alpha可以大于0 Cm_q -5.2; Cm_delta -0.9; % 舵面效率 CD_alpha 0.18; CL_delta 0.32; qbar 0.5 * rho * V0^2; % 力与力矩对状态和输入的偏导数 Xu -0.03; Xalpha 0.2; Xq 0; Xdelta 0; Zu -0.15; Zalpha -qbar*S*CL_alpha/(m*V0); Zq 0; Zdelta -qbar*S*CL_delta/(m*V0); Mu 0; Malpha qbar*S*cbar*Cm_alpha/Iy; Mq qbar*S*cbar^2*Cm_q/(2*Iy*V0); Mdelta qbar*S*cbar*Cm_delta/Iy; A [Xu, Xalpha/57.3, Xq, -9.8*cosd(alpha0); Zu/V0, Zalpha/V0, Zq/V01, -9.8*sind(alpha0)/(V0); Mu, Malpha/57.3, Mq, 0; 0, 0, 1, 0]; B [Xdelta, 0.1; Zdelta/V0, 0; Mdelta, 0; 0, 0]; C eye(4); D zeros(4,2); trim struct(V0, V0, alpha0, alpha0, elev, 0.02); end这个函数里A阵的Malpha和Mq直接决定短周期模态。飞翼的Cm_alpha是正时开环极点会出现在右半平面直接仿真发散。C取eye(4)意味着四个状态都可测如果实际只有部分状态比如没有空速需要把C矩阵改成可测量组合后面H∞设计也要相应调整。这里x的单位用角度和弧度混着注意在注释里标明。提示fland.m返回的矩阵单位不一致会导致极点数量级错误。一般做法是角度全部转成弧度再进入A阵最后输出到曲线时再转回角度。2.2 不确定性来源与乘性不确定性集合飞翼无人机的不确定性有几个来源气动导数的风洞误差、重心位置随飞行变化、舵机动态随温度变化、机翼弹性模态。鲁棒控制不把这些问题逐个做“保守的修正”而是用一个摄动块Δ把它们全部包进去。最常用的是乘性不确定性描述G_p(s) G_0(s) [1 W_t(s) Δ(s)], ||Δ||_∞ ≤ 1G_0是标称模型W_t是加权函数它的幅值表示不同频率下模型误差的相对大小。低频段重心变化导致的增益偏差高频段弹性模态和舵机延迟导致的相位偏差都可以体现在W_t的幅值曲线上。这样设计出来的控制器不需要知道Δ的具体值只要它小于等于1就能保证稳定。不确定性来源主要影响通道权函数特征气动导数摄动短周期极点偏移W_t在低频段约0.1~0.3重心前/后移俯仰力矩系数低频幅值随移动量增大舵机延迟与速率饱和输入相位滞后高频幅值超过1机翼弹性模态输入输出耦合在弹性频率处出现峰值实际操作时先用风洞数据或试飞数据把不同状态下的频率响应画出来再看标称模型相对真实模型的最大误差包络用一阶或二阶传递函数去拟合这个包络。makeweight就是干这件事的函数。2.3 用makeweight构造性能与不确定性加权函数增广对象里需要两类权函数性能加权W_p和不确定性加权W_t。W_p的物理意义是要求系统对低频扰动阵风、重心偏移有足够大的开环增益使闭环灵敏度函数S1/(1GK)小W_t则要求在高频段限制补灵敏度函数T防止未建模动态被激活。二者是矛盾的H∞设计就是折中。makeweight.m提供了一种快速构造权函数的方式。三个参数分别是低频增益、转折频率、高频增益% 性能权低频要求1000倍扰动抑制5 rad/s处交越高频放松到0.1 Wp makeweight(1000, 5, 0.1); % 不确定性权低频误差10%50 rad/s之后误差超过200% Wt makeweight(0.1, 50, 2); figure; bode(Wp, b, Wt, r, {0.01, 1000}); grid on; legend(Wp (性能), Wt (不确定性));makeweight生成的传递函数在低频段增益趋向第一个参数在转折频率附近穿过0dB如果参数设置合适高频段趋向第三个参数。注意第一个参数是直流增益不是0dBWp的1000意味着低频开环增益要达到60dB才能满足性能要求。如果计算出来的控制器奇异值曲线压不住Wp的倒数就要降低Wp低频增益或放宽转折频率。3. H∞、滑模、自适应三种鲁棒控制器在Matlab中的实现3.1 基于hinfstruct的固定结构H∞综合经典H∞用hinfsyn会得到全阶控制器阶数与增广对象相同对机载计算机来说往往太高。hinfstruct的优势是可以指定控制器结构比如PID、二阶状态空间然后做结构化H∞综合得到的控制器可以直接落代码。资源里的学术论文.docx把增广对象的推导写得很清楚Matlab侧调用方式如下% 构造增广对象输入为 G输出为 [Wp*e; Wt*u; y] 的形式 P augw(G, Wp, Wt, []); % 可调控制器2输入2输出4个状态 C0 ltiblock.ss(C, 2, 2, 4); % 固定结构H∞综合 [C, CL, info] hinfstruct(P, C0);需要注意几个点。第一P的输入是G的控制输入和外部扰动输出是加权误差、加权控制量和反馈测量量augw会自动按照这个规则拼接具体增益和输入输出顺序可以通过P.InputName/P.OutputName验证。第二hinfstruct是非凸优化依赖多起点搜索运行一次可能落在局部最优工程上要多次改变C0的初始值保留info.gamma最小的一组。第三返回的C是一个可调控制块对象可以用getValue取出当前控制器也可以直接和G做feedback得到闭环。如果最后需要连续状态空间矩阵用ss(C)转换。要验证闭环稳定性和鲁棒裕度把C的频响与Wt、Wp的倒数画在一起看bode(ss(C), {0.01, 1000}); grid on;hinfstruct属于Robust Control Toolbox如果用的是Matlab在线版或较老版本需要确认工具箱已经授权否则会报undefined function。这与安装不相关是工具箱许可问题。3.2 滑模控制器设计从切换面到趋近律滑模控制对匹配不确定性不敏感意思是当不确定性满足B范-matched条件时只要切换增益足够大系统状态就会被吸引到切换面s0上。飞翼内环姿态控制习惯用滑模做角度和角速率控制。对于二阶系统x[α, q]取滑模面s c1(α - α_ref) (q - q_ref)控制律由等效控制加切换控制组成c1 5; % 滑模面权重 x_ref [alpha_ref; q_ref]; s c1*(x(1)-x_ref(1)) (x(2)-x_ref(2)); % 等效控制令 ds/dt 0 解出的连续控制 u_eq -( (c1*B(1) B(2)) )^(-1) * (c1*A(1,:) A(2,:)) * (x - x_ref); % 切换控制epsilon抵抗不确定性k加速趋近 epsilon 0.2; k 3; u_sw -((c1*B(1) B(2)))^(-1) * (epsilon * sign(s) k * s); u u_eq u_sw;这里的sign会带来抖振飞翼舵机受不了高频切换。常见处理是用饱和函数sat(s/φ)代替sign边界层厚度φ取0.05左右。epsilon如果给得小于实际不确定性幅值滑模面会存在稳态误差这不是参数整定问题而是可达性条件没满足。要检查s*s_dot是否小于0在仿真里把这个量打出来看。滑模控制器对建模依赖比H∞弱但没有直接给出稳定裕度量化所以工程上常用它做内环、H∞做外环或者用滑模控制某一通道。3.3 自适应控制与参数更新律自适应控制针对的是慢变参数比如气动导数随马赫数的变化、重心随燃油消耗的移动。最简单的是模型参考自适应控制MRAC。让实际系统跟踪一个参考模型xdot_ref Aref x_ref Bref r误差e x - x_ref参数θ按梯度法更新。% MRAC参数更新Gamma为学习率P为Lyapunov方程的解 A_ref [-2, 1; -10, -3]; B_ref [0; 1]; P lyap(A_ref, eye(2)); Gamma 0.5; Phi [x(1); x(2); r]; % 回归向量 theta_dot -Gamma * Phi * (e * P * B_ref); theta theta - theta_dot * dt; % 实际离散积分 u theta * Phi;注意这里的Phi和theta维度需要与控制通道一致。MRAC的缺点是瞬态响应可能有大超调且对高带宽扰动敏感。飞翼这个场景里我更倾向把自适应用来在线补偿气动导数误差而不是取代H∞主回路。资源中的PID.m和Duffing.m分别对应基准控制器和被控对象的非线性振动测试可以在同一个框架下对比。3.4 三种方案的选型对比控制器设计工具不确定性处理阶数工程难点H∞hinfstruct / hinfsyn显式建模加权函数中到高权函数选择滑模手写s函数Simulink匹配不确定性有界低抖振与可达性自适应手写MRAC慢变参数在线估计低瞬态超调和激励条件选型上没有绝对最优。H∞适合给出可承诺的稳定裕度滑模适合快速内环自适应适合在线配平修正。实际项目通常混用而不是只用一种。4. 从main.m到Monte Carlo开环闭环对比与不确定性统计4.1 仿真主流程与noise.m噪声注入main.m是整个工程的主入口工程中一般把模型初始化、控制器生成、仿真与绘图分成四个部分。一个比较标准的流程是先用fland得到标称模型再用makeweight构造权函数然后调用hinfstruct得到控制器之后在simulink或循环脚本里做时域仿真。噪声部分单独放在noise.m里。Duffing.m通常会往被控对象里注入一个非线性扰动项用来观察控制器在极限环或结构非线性下的鲁棒性。% main.m 入口脚本节选 rng(2024); dt 0.01; T 10; t (0:dt:T); % 1. 标称模型 G0 ss(A, B, C, D); % 2. 不确定性摄动对气动导数做±20%的随机变化 [A_pert, B_pert] fland_perturbed(0.2); G ss(A_pert, B_pert, C, D); % 3. 外部扰动noise.m 生成阵风序列 gust noise(t, gust, 3); % 第三个参数表示阵风强度 dist_u zeros(length(t), 2); dist_u(:, 2) gust; % 把阵风作用在推力通道noise.m的实现可以是高斯白噪声通过一阶成形滤波器也可以直接生成Dryden风模型。常见做法是bandlimited white noise加低通滤波这样更接近大气紊流功率谱。如果发现仿真曲线高频抖得厉害先看noise.m生成的噪声是否已经经过合适的滤波不要在控制器里强行加低通掩盖数值问题。4.2 开环不可飞、闭环可增稳的响应对比飞翼开环在扰动下的响应是发散的因为Cm_alpha为正或静稳定度过低。用lsim分别对开环和闭环打同样的阶跃输入能够直观看到控制器是否把极点拉回左半平面。代码% 开环响应 y_open lsim(G, dist_u, t); % 闭环响应控制器需事先用 ss(C) 转成状态空间 sys_cl feedback(G * ss(C), eye(2), 1); y_cl lsim(sys_cl, dist_u, t); % 画迎角通道 figure; subplot(211); plot(t, y_open(:,2)); title(Open-loop alpha); subplot(212); plot(t, y_cl(:,2)); title(Closed-loop alpha);feedback里的1是正反馈标志。飞控系统的舵面极性如果定义成“负反馈与正反馈”这里最容易出错开环发散而闭环也会因为符号不对直接正反馈发炸。判断方法很简单把舵面偏角输出打出来看它是不是在抵抗扰动。飞翼舵效本来就低正反馈几秒钟就会让姿态发散。4.3 Monte Carlo参数摄动下的统计验证单次仿真的时域曲线只能说明一个摄动样本下控制器有效不能说明鲁棒性。Monte Carlo方法是在给定的摄动范围内随机采样N组A、B矩阵每组都跑一次闭环仿真统计关键指标。这样能把“最坏情况”从经验判断变成数据。N 300; OS zeros(N,1); Ts zeros(N,1); MaxDev zeros(N,1); for i 1:N sigma 0.15 * randn; % 标准差15% [Ai, Bi] fland_perturbed(sigma); Gi ss(Ai, Bi, C, D); yi lsim(feedback(Gi * ss(C), eye(2), 1), dist_u, t); OS(i) max(yi(:,2)); % 迎角超调峰值 idx find(abs(yi(:,2) - ref) 0.1*abs(ref), 1, last); if isempty(idx), Ts(i) T; else Ts(i) t(idx); end MaxDev(i) max(abs(yi(:,2) - ref)); end fprintf(超调峰值 mean%.3f std%.3f\n, mean(OS), std(OS));指标含义设计建议迎角超调峰值最大迎角与参考值之差通常要求小于1.5°调节时间进入5%误差带的时间飞翼不宜超过3s最大偏差对参考曲线的最大偏离反映阵风抑制能力闭环特征值实部稳定性余量实部应小于-0.5如果Monte Carlo中有5%的样本超调翻倍不要直接增大控制增益那样只会让高频段变差。正确做法是回去调整Wt或者Wp的交越频率重新做H∞综合。这也体现了鲁棒控制“设计在频域、验证在时域”的闭环流程。5. 工程落地技巧用bode图和控制器降阶守住裕度5.1 用bode图和灵敏度函数验证鲁棒裕度仿真通过不代表能装机。飞翼无人机对相位裕度敏感工程上通常要求开环增益裕度大于6dB、相位裕度大于45°。用bode命令把控制器和对象串联后的开环传递函数画出来直接读裕度L G0 * ss(C); % 开环注意G0是标称模型 bode(L); grid on; [Gm, Pm] margin(L); fprintf(增益裕度 %.2f dB相位裕度 %.2f°\n, 20*log10(Gm), Pm);如果是多输入多输出系统margin对每个通道分别给出结果还要看灵敏度函数S的最大奇异值是否小于Wp倒数的包络。使用sigma(S)可以检查多变量系统的鲁棒性能单一Bode图在多通道耦合时不够充分。5.2 对高阶H∞控制器降阶并检查逼近误差hinfstruct如果不限定状态数得到的控制器阶数一般比标称模型高。机载飞控计算机的主频有限而且高阶控制器在高频段容易放大测量噪声。可以用balred做降阶然后对比降阶前后的闭环频响。C2 balred(ss(C), 2); % 降到2阶 bode(ss(C), C2); grid on; % 看中低频段是否贴合 S1 feedback(1, G0 * ss(C)); S2 feedback(1, G0 * C2); sigma(S1, S2); grid on; % 降阶后的灵敏度不能明显抬升降阶不是越低越好。工程中先保留4阶再逐步减观察灵敏度函数的峰值是否超过2dB。飞翼的操纵面少控制器内部还可能要加陷波器所以实际部署前要把降阶后的控制器连同预滤波器一起做硬件在环仿真。5.3 连续控制器离散化与部署验证Matlab里用c2d对控制器离散化采样时间要低于飞行控制频率的5~10倍。飞翼通常用200Hz或400Hz控制频率Tustin变换比ZOH更适合保留相位特性Ts 1/200; Cd c2d(C2, Ts, tustin); bode(C2, Cd); grid on;离散化之后要重新检查闭环稳定裕度因为Tustin变换在高频段会略微改变相位。如果高频段裕度下降可以考虑把控制器和对象放在同一个离散框架里做H∞综合而不是连续综合后离散化。最后把Cd转成差分方程放进飞控代码时注意输入输出定标和保护限幅飞翼舵面效率低限幅太小会降低扰动抑制能力限幅太大又可能激励弹性模态。本文还有配套的精品资源点击获取