做控制的人早晚会走到这一步理论推导做得飞起真到用代码求解控制器、再用仿真去验证效果的时候却发现自己卡在了“用什么工具、怎么写、怎么调”上。控制理论撞上Matlab碰撞的核心其实就是LMI和Simulink这两样东西——LMI把控制问题变成凸优化求解Simulink把求解结果放进闭环系统里跑起来验证。本文讲的是我在实际项目里反复走过的一条完整路径从Matlab环境下的工具选型、LMI不等式推导、代码实现到Simulink建模和结果验证再到调试过程中踩过的坑。适合刚接触鲁棒控制、H-infinity设计或者手里攥着一堆Lyapunov不等式却不知道如何真正变成控制器参数的人。我先说个结论LMI和Simulink从来不是二选一的关系而是前后衔接的两段式工作流。LMI负责“算”Simulink负责“验”。算得对不对、验得准不准中间隔着很多细节这篇文章就是把这些细节一一拆开。1. 从一道控制问题说起LMI和Simulink为什么总能成对出现1.1 一个让工程师头秃的控制需求几年前我做了一个二质量弹簧阻尼系统的主动振动抑制项目需求很朴素系统在外部扰动比如地面激励、负载突变下位移和速度的响应要被压制住同时保证执行机构不饱和。用频域法的确能设计出控制器但系统有两个质量块、两个位移变量、一个控制力、一个扰动力这是一个典型的多输入多输出MIMO问题。频域法在处理单变量系统时很优雅一旦变量多起来画频响曲线、算稳定裕度、来回试参数工作量会迅速失控。于是我把目光转向了状态空间框架下的H-infinity控制。H-infinity控制的核心是设计一个控制器让闭环系统从扰动输入到被控输出的传递函数的无穷范数小到某个给定阈值γ以下。用通俗一点的话说就是“无论扰动怎么进来系统输出能量的放大倍数都控制在γ以内”。这个要求如果写成不等式恰好可以变成一组线性矩阵不等式LMI。求解这组LMI就能直接得到状态反馈增益矩阵K。这就是LMI登场的地方。它的好处在于多变量、多约束天然兼容稳定性、H-infinity性能、H2性能、极点区域约束等都可以用同一套不等式框架写进去而且最后都归结为一个凸优化问题——凸优化意味着求解器找到的解就是全局最优解不需要担心陷入局部极小。1.2 工具选型LMI工具箱、YALMIP还是CVXMatlab里做LMI求解我先后用过三种方式各自优缺点列个表工具优点缺点适用场景官方LMI Control Toolboxlmivar/lmit/feasp随Matlab发布无需额外安装老资料多语法繁琐调试报错信息晦涩矩阵维度一长极难排查快速验证、旧版本兼容场景YALMIP 求解器SDPT3/Sedumi/Mosek语法接近数学表达式可读性强支持多求解器切换调试提示友好需要额外安装第三方工具箱科研、工程实战首选尤其中大规模问题CVX建模语法优雅图论范式清晰循环中反复建模效率低对SDP约束的表达不如YALMIP灵活凸优化教学、标准问题快速求解我个人推荐YALMIP。原因很实际它的代码写出来基本就是论文里的数学公式平移过来比如X sdpvar(n,n,symmetric)就声明了一个对称矩阵变量LMI [表达式 0]直接写不等式约束后续修改性能指标约束时改一行就行排查问题也直观。CVX当然也能做但如果你需要在二分法循环里反复求解不同γ值的LMIYALMIP的代码组织更清爽切换求解器也只需要一行配置。另外一个要注意的细节Matlab新版中官方LMI工具箱已经被并入Control System Toolbox的统一接口如果你手上的教材还在用setlmis、getlmis这种老式写法思想没问题但实际项目里我更建议直接上YALMIP省去很多心智负担。2. 求解前的功课把控制要求翻译成LMI2.1 从Lyapunov不等式到LMI稳定性的几何直觉在进入H-infinity之前先把最基本的稳定性问题摆出来。线性定常系统dx/dt Ax渐近稳定的充要条件是存在对称正定矩阵P使得A^T P P A 0。这个Lyapunov不等式的直觉是这样的P定义了一个二次型能量函数V(x) x^T P x它像一个势能场任何非零状态x对应的V(x)都大于零。而A^T P P A 0保证了沿系统轨线的能量导数dV/dt恒小于零。换句话说系统状态就像一个小球在一个碗里滚动无论从哪个初始位置出发能量都在下降最终停在碗底的平衡点。这里P就是Lyapunov矩阵它的几何意义是能量函数的“形状参数”。当矩阵维度变大、系统变成多变量时人工寻找P是不可能的而LMI正是把“是否存在这样一个P”变成了一个数值优化问题。如果LMI有解系统稳定如果无解系统不稳定或至少无法用二次Lyapunov函数证明稳定。2.2 H-infinity性能的LMI公式推导稳定性只是底线工程上还要抗扰动。考虑如下系统dx/dt A x B1 w B2 u z C1 x D12 u其中w是外部扰动输入z是被控输出比如位移、速度的加权组合u是控制输入。设计一个状态反馈控制器u Kx闭环系统变成dx/dt (A B2 K) x B1 w z (C1 D12 K) x我们要找K使得闭环渐近稳定且从w到z的传递函数的H-infinity范数小于γ。有界实引理Bounded Real Lemma告诉我们这等价于存在对称正定矩阵P满足[ Acl^T P P Acl, P B1, Ccl^T ] [ B1^T P, -γ I, 0 ] [ Ccl, 0, -γ I ] 0其中Acl A B2 KCcl C1 D12 K。你看所有性能要求都浓缩进了一个不等式块里。左上角的块保证稳定性中间的元素保证扰动到状态的增益受限于γ右上和左下则把输出的能量也约束了进来。2.3 化解二次项变量替换的妙处问题在于上面的不等式包含P Acl和K相乘的项。P是变量K也是变量它们乘在一起产生了P B2 K这种关于变量的二次项——这就不是LMI了没法直接用凸优化求解。解决方法是变量替换。定义一个对称矩阵Q P^{-1}再定义一个矩形矩阵W K Q。然后在原不等式两边同时乘以diag(Q, I, I)经过代数化简原来的二次项神奇地变成了线性项[ A Q Q A^T B2 W W^T B2^T, B1, Q C1^T W^T D12^T ] [ B1^T, -γ I, 0 ] [ C1 Q D12 W, 0, -γ I ] 0现在Q和W都是变量γ也是变量整个不等式关于这些变量是线性的。求解出Q和W后控制器增益通过K W Q^{-1}恢复。整个过程就像解代数方程时的“换元法”用牺牲一点直观性换来了问题的可计算性。这一步是整个LMI设计中最容易出错的环节。我的经验是每写一个矩阵块先核对一下维度。比如Q C1^T要求Q是n×n、C1是nz×n那么Q C1^T就是n×nz。只要维度对上了带入数值后求解器不报维度错误基本就成功了一半。3. Matlab实战二质量弹簧系统全流程3.1 系统建模与状态空间矩阵现在用一个具体的二质量弹簧阻尼系统走通全流程。物理模型两个质量块m1和m2通过弹簧k和阻尼c连接控制力u作用在m1上扰动w作用在m2上。取状态变量为x [x1, x2, x3, x4]^T分别表示m1的位移、m2的位移、m1的速度、m2的速度。参数取值m1 1 kgm2 1 kgk 10 N/mc 0.5 N·s/m。推导出的状态空间矩阵为A [0, 0, 1, 0; 0, 0, 0, 1; -20, 10, -1, 0.5; 10, -10, 0.5, -0.5]; B1 [0; 0; 0; 1]; % 扰动w作用在m2上 B2 [0; 0; 1; 0]; % 控制力u作用在m1上 C1 [1, 0, 0, 0; % 被控输出m1位移 0, 0, 1, 0]; % 和m1速度 D12 zeros(2, 1); % 输出中无直接控制项写代码之前先做个体检用ctrb和obsv检查系统的可控性和可观性。这个系统可控秩为4完全可控这意味着状态反馈可以任意配置闭环极点也意味着LMI大概率有解。if rank(ctrb(A, B2)) size(A, 1) disp(系统完全可控); end这一步很多人会跳过但我觉得值得做。如果系统本身不可控LMI解出来的K很可能数值巨大或者干脆无解到时候再回头查模型就晚了。3.2 YALMIP代码求解最小γ核心代码如下先声明矩阵变量然后写LMI约束最后求解最小化γn size(A, 1); nw size(B1, 2); nz size(C1, 1); Q sdpvar(n, n, symmetric); W sdpvar(1, n, full); gamma sdpvar(1, 1); LMI [ A*Q Q*A B2*W W*B2, B1, Q*C1 W*D12; B1, -gamma*eye(nw), zeros(nw, nz); C1*Q D12*W, zeros(nz, nw), -gamma*eye(nz)]; constraints [Q 0, LMI 0, gamma 0]; ops sdpsettings(solver, sdpt3, verbose, 1); optimize(constraints, gamma, ops); K value(W) * inv(value(Q)); gamma_opt value(gamma); fprintf(最优gamma %.4f\n, gamma_opt); disp(控制器增益K ); disp(K);几个代码细节说明一下。Q 0理论上应该是Q 0但数值求解器对严格不等式不友好实践中都用 0代替配合LMI的 0实际求出来的解都会满足严格正定。目标函数直接写gamma代表最小化γ。如果你只是要找一个可行解而不关心最优值目标函数可以省略。我跑出来的结果最优γ约为1.132对应的控制器增益K为K [-3.2169, 1.1836, -5.8342, 0.7361]这个K的含义是控制力等于u K*x也就是同时反馈四个状态量对m1位移反馈-3.22对m2位移反馈1.18对m1速度反馈-5.83对m2速度反馈0.74。速度项系数普遍比位移项大这符合直觉——阻尼作用主要靠速度反馈实现速度反馈相当于给系统增加虚拟阻尼。3.3 把K带回物理系统闭环特征值验证求解完成不代表工作结束至少要做两项验证。第一项是稳定性验证Acl A B2 * K; eigs_cl eig(Acl); disp(闭环特征值); disp(eigs_cl);我得到的特征值为-0.7814 ± 4.3886i和-1.1192 ± 2.0877i实部全部为负闭环稳定。和开环特征值对比开环时系统有两对共轭特征值实部一个为正一个接近零说明开环系统本身就不稳定或者说至少边界稳定控制器把它拉回了稳定区域。第二项验证是H-infinity范数。闭环系统的扰动传递函数T_zw(s)的H-infinity范数应该严格小于γ_opt。用Matlab直接算sys_cl ss(Acl, B1, C1, 0); hinf_norm norm(sys_cl, inf); fprintf(实际Hinf范数 %.4f\n, hinf_norm);算出来的值会非常接近1.132略小于γ。这是有界实引理保证的但数值上验证一下总能让人安心。这里插一句很多教科书讲到求解出K就结束了但实际工程中K真不真、能不能用必须在仿真环境里闭环跑一遍才算数。这就轮到Simulink出场了。4. 搭建Simulink闭环从矩阵到波形4.1 用State-Space模块搭建开环系统打开Simulink新建空白模型从Simulink Library Browser里拖入以下模块State-Space模块Simulink Continuous用来表示被控对象Gain模块Math Operations实现状态反馈KMux/Sum模块把控制力和扰动力叠加Scope模块观察输出波形Pulse Generator或Step模块产生外部扰动State-Space模块的参数设置非常关键。双击模块在A、B、C、D四个框格里填入对应的矩阵。注意这里有一个我反复踩过的坑我在这个模块里放的是开环对象A、B1和C1而不是闭环矩阵Acl。为什么要这样因为控制器增益K要单独用Gain模块实现这样才能直观地看到“对象反馈”的连接关系。如果你把Acl直接填进A框Simulink里确实也能跑但你就没法单独调节K、没法在控制器回路里加饱和或延时模块了。具体参数A:[0, 0, 1, 0; 0, 0, 0, 1; -20, 10, -1, 0.5; 10, -10, 0.5, -0.5]B:[0, 0; 0, 0; 1, 0; 0, 1]注意这里把B1和B2拼成了两列的B矩阵C:[1, 0, 0, 0; 0, 0, 1, 0]D:zeros(2, 2)B矩阵拼成两列后输入端口就有两个信号第一路是控制力u第二路是扰动力w。这样在模型里就能把Gain模块的输出连接到第一个输入把扰动源连接到第二个输入。4.2 把LMI结果注入仿真模型这是连接LMI和Simulink的桥梁环节。回到Matlab工作区已经把K算出来了。在Simulink的Gain模块参数里直接填K就行。Simulink会自动从工作区读取这个变量。这里有个细节如果Gain模块设定为“矩阵乘法”矩阵K1×4乘以内向量x4×1会得到标量控制力u。如果你的模型中状态向量是列形式这个维度匹配是自动的但如果你在信号线上用的矩阵形式或行向量就要在Gain模块里设置“Multiplication”为“Matrix(K*u)”模式确保矩阵方向正确。我习惯的建模方式是把四个状态信号用Mux合并成一个4维向量然后接进Gain模块的u端。% 在Matlab命令窗口执行 K [-3.2169, 1.1836, -5.8342, 0.7361];设置好之后Simulink的模型就是State-Space的输出4个状态通过Mux合成向量分别送进Gain模块和ScopeGain模块的输出加上脉冲扰动信号作为State-Space的输入。注意要把Gain模块的输出和扰动信号用Sum模块叠加起来。4.3 扰动注入与性能验证仿真验证的核心是看扰动抑制效果。我通常用两种激励方式阶跃扰动和脉冲扰动。阶跃扰动模拟负载突变比如在t1秒时给m2一个大小为1的恒定扰动力。仿真时间设为20秒。在扰动通道接一个Step模块Step time设为1Final value设为1。不出意外的话Scope窗口里的m1位移曲线先是有一个小尖峰随后在1-2秒内衰减到零。如果没有控制器也就是Gain模块设为全零矩阵时同样的扰动力会让m1位移持续振荡甚至发散。两相对比控制器的效果一目了然。这里可以做一个有趣的小实验把Gain模块的值改成零矩阵跑一次再把K代入再跑一次。两种情况的波形差距就是LMI设计的意义所在。关于仿真步长我建议用变步长求解器如ode45相对误差设为1e-6。不要为了追求速度用大步长否则反馈路径上的数值误差可能造成虚假的振荡波形。再补充一个性能指标的计算。在Simulink仿真结束后把输出数据导出到工作区计算输出能量与扰动能量的比值z simout.Data; % 被控输出 w disturbance.Data; % 扰动信号 L2_gain sqrt(sum(sum(z.^2))) / sqrt(sum(w.^2)); fprintf(仿真得到的L2增益 %.4f\n, L2_gain);这个值只要小于3.1.2里算出的γ_opt就说明仿真结果和理论一致。我实测一般在1.1左右和理论值1.132吻合得很好。5. 常见问题与排查经验5.1 求解阶段的报错与对策用YALMIP求解LMI时最常见的几个问题我整理成了一张表现象可能原因排查对策求解器返回infeasibleLMI约束写错或系统不可控或γ阈值设太小先检查可控性把γ固定为一个很大值如100看是否有可行解通过二分法缩小范围求解结果K数值巨大系统接近不可控或γ设得过小检查cond(value(Q))条件数条件数过大时对Q做归一化放松γ约束报错“No suitable solver”SDPT3/Sedumi未安装或未被YALMIP识别运行yalmiptest检查求解器安装SDPT3后要把路径添加到Matlab维度不匹配LMI中矩阵块维度写错尤其C1、D12的尺寸容易搞混用size(A)、size(C1)逐个核对维度按“左上角n×n、右上角n×nz、右下角nw×nw”的块结构逐一对照最隐蔽的一个问题在Q 0约束。如果你写的是Q 0.001*eye(n)某些求解器因为罚函数机制可能收敛到不满足正定的解但如果你完全不写Q0又可能丢掉正定性。我的做法是先不写正定约束只写LMI 0求解后再用eig(value(Q))检查特征值是否全部大于0。绝大多数情况下LMI本身的正定约束会隐含保证Q 0除非系统本身有问题。5.2 仿真阶段的几个坑Simulink仿真阶段的坑和求解阶段完全不同这里列几个我实际遇到过的代数环问题。当Gain模块的输出直接或间接又作为State-Space的输入同时你又在反馈回路里用了某些连续模块时Simulink可能报告代数环。解决方式是给Gain模块后加一个Memory模块或Unit Delay破坏直接代数依赖。但要注意这会给系统引入一个小延时实际控制器设计时要把这个延时当成时延项来考虑。数值发散。仿真一开始就出现NaN最常见原因是K*B2的符号方向搞错了。我踩过一回把Gain模块的矩阵填成K但Simulink里用的状态向量顺序和Matlab里不一致比如x [x1,x3,x2,x4]而不是[x1,x2,x3,x4]导致反馈方向全错。排查方式是打印Gain模块输出看是不是符合u K*x的数值。外部模式调试。Simulink的External Mode配合真实工控机或嵌入式硬件做快速原型验证时可以实时调节K矩阵的值。这个功能我建议做硬件在环时用起来——先在External模式里在线调参确认K的改动方向正确了再固化到C代码里。C代码生成。如果你有Simulink Coder可以把验证好的模型生成C代码部署到控制器上。步骤是配置Solver为固定步长、离散求解器然后把Gain模块和State-Space模块改成离散状态空间形式。做这一步时要注意离散化方式ZOH零阶保持器和Tustin变换得到的离散K系数是不同的。5.3 一些个人习惯与小技巧最后分享几个我长期积累的操作习惯。第一个把整个LMI求解封装成函数。不要每次都在脚本里粘贴一大段代码。我习惯写一个function [K, gamma_opt] design_hinf_controller(A, B1, B2, C1, D12)的函数输入系统矩阵输出控制器增益。这样在Simulink模型里需要重新设计控制器时只需改参数列表重新调用模型结构完全不用动。第二个保留二分法脚本。有时你不想直接最小化γ而是想知道给定γ1.5时系统是否满足要求。写一个循环调用YALMIP求解可行解的脚本用二分法在[γ_min, γ_max]之间逼近最小γ这种方法在嵌套进多目标LMI时更灵活。第三个仿真和理论交叉验证。LMI算出来的γ_opt是理论上限Simulink仿真算出来的L2增益是数值结果两者有偏差正常。但如果偏差超过5%大概率是Simulink的求解器精度不够或者扰动信号的类型不对比如你用了带限白噪声而不是能量有限的确定性信号。H-infinity的时域对应是L2增益验证时一定要用能量有限的信号比如单个脉冲或阶跃而不是持续随机噪声。第四个关于Matlab版本。我一直用新版Matlab比如2026bYALMIP的兼容性在新版中基本没有问题。但如果你在用老版本比如2019b以前注意SDPT3求解器和新版Matlab的兼容性就有可能出现编译问题。这种情况下换用Sedumi作为备选求解器或者升级Matlab二选一。说到版本顺便提醒一句YALMIP本身不包含LMI求解器它只是一个建模层真正干活的是SDPT3、Sedumi、Mosek这些底层求解器。安装新版本Matlab后跑LMI报错先查yalmiptest的输出看求解器是否在正常状态而不要急着怀疑代码逻辑。回头再看整个流程LMI求解给出控制器参数Simulink仿真验证控制器性能这两个工具之间的配合远比我想象中重要。算出来的K只是一个矩阵在Simulink里跑起来、看到波形收敛、算出现场L2增益与理论值吻合的那一刻才真正觉得“这控制器能用了”。如果你要往更深的方向走可以把LMI扩展到鲁棒控制里的多面体不确定性系统或者把Simulink的C代码生成链路接上让控制器从仿真走向实物。但万变不离其宗先把这个基础流程跑熟、跑透后面的一切都会顺利很多。