1. 项目概述为什么LMI工具箱是解决复杂优化问题的利器在控制系统设计、鲁棒滤波、信号处理乃至金融工程等领域我们常常会遇到一类“棘手”的优化问题目标函数或约束条件中包含了矩阵不等式。这类问题传统的基于梯度的优化算法如fmincon往往力不从心要么难以求解要么得到的解可靠性存疑。这时线性矩阵不等式Linear Matrix Inequality, LMI理论及其在MATLAB中的实现——LMI工具箱就成为了工程师和研究员手中的一把瑞士军刀。简单来说LMI就是一个关于矩阵变量的线性不等式其形式通常为 F(x) F0 x1F1 ... xnFn 0其中“ 0”表示矩阵是负定的。许多复杂的系统分析与设计问题如判断系统稳定性、求解线性二次型最优控制LQR的Riccati方程、设计H∞控制器等都可以转化为一个或多个LMI的可行性问题或优化问题。MATLAB的LMI工具箱正是为了高效、可靠地求解这类问题而生的。然而工具箱的强大也伴随着一定的入门门槛。其核心函数lmivar和lmiterm的语法相对独特需要使用者以一种“声明式”的思维来构建问题这与我们熟悉的命令式编程有所不同。网上很多资料要么过于理论化要么示例零散让初学者望而却步。本文将从一个实践者的角度手把手带你拆解这两个核心函数并通过一个完整的控制器设计案例展示如何从问题描述到代码实现最终求解出可靠的解。无论你是正在做课题的研究生还是需要解决实际工程问题的工程师这篇超详细的教程都将为你提供一条清晰的路径。2. LMI问题构建的核心思维从数学描述到工具箱语言在动手写代码之前我们必须先建立正确的思维模型。LMI工具箱的求解流程类似于我们向一个“求解器”提交一份“问题说明书”。这份说明书需要明确两点变量和约束。lmivar函数负责定义变量lmiterm函数负责描述约束中的每一项。整个构建过程是“分而治之”的先定义所有矩阵变量再针对每一个LMI约束逐一累加其构成项。2.1 理解LMI的标准形式LMI工具箱处理的问题一般形式为 找到矩阵变量 X1, X2, ..., Xk使得一组线性矩阵不等式成立 N1^T * L(X1, ..., Xk) * N1 M1^T * R(X1, ..., Xk) * M1 N2^T * L(X1, ..., Xk) * N2 M2^T * R(X1, ..., Xk) * M2 ... 其中 L(.) 和 R(.) 是矩阵变量的仿射函数。在大多数基础应用中我们遇到的是更简单的形式A^T * X * A - X Q 0或A*X X*A^T B*B^T 0等。工具箱要求我们将这些不等式都以“小于零”的形式写在等号的同一侧。例如不等式P 0P正定等价于-P 0。不等式A^T*P P*A Q 0则已经是在“小于零”的一侧了。2.2 变量定义函数 lmivar 详解lmivar用于定义LMI问题中的未知矩阵变量。其基本调用语法为[var_id, n, s] lmivar(type, struct)var_id: 输出参数是该变量的句柄ID号在后续的lmiterm中通过此ID来引用该变量。n: 输出参数当type1时表示变量的维度当type2或3时表示该变量在最终决策变量向量中的块大小可暂时忽略。type: 输入参数定义变量类型是关键所在。type 1: 对称块对角矩阵变量。这是最常见的情况用于定义对称矩阵如Lyapunov函数中的P矩阵。struct参数是一个 k×2 的矩阵每一行描述一个对角块。例如struct [n, 0]: 表示一个 n×n 的满对称块。struct [n, 1]: 表示一个 n×n 的标量矩阵即单位矩阵乘以一个标量变量。struct [n, -1]: 表示一个 n×n 的全零块占位无变量。type 2: 长向量变量即标准的决策向量 x。struct [m, n]表示一个 m×n 的矩形矩阵变量它将被按列拉直成一个长向量。这种类型在将传统优化问题转化为LMI形式时有用。type 3: 其他结构。struct是一个与变量维度相同的矩阵其中每个元素只能是 0, x, -x。0表示该位置固定为0x表示该位置是一个新的标量决策变量-x表示该位置是-1乘以一个已由x定义的变量用于强制对称性或特殊结构。这种类型最灵活但也最复杂。实操心得一变量定义的策略对于初学者我的建议是优先使用type1来定义所有对称矩阵变量。除非你非常确定需要type2或3否则它们容易引入错误。在定义type1时想清楚你的矩阵是“完全自由的对称矩阵”用[n, 0]还是“单位矩阵的倍数”用[n, 1]。后者可以显著减少变量个数加快求解速度适用于某些特定结构例如你需要寻找一个公共的Lyapunov矩阵比例因子。2.3 项描述函数 lmiterm 详解lmiterm用于向指定的LMI约束中添加一项。这是构建LMI最核心、也最容易出错的一步。其语法为lmiterm(termID, A, B, flag)termID: 一个四元向量[p, q, x, y]它指明了当前项添加到哪个位置。p: LMI的编号。正数表示第p个“小于零”不等式 0负数表示第p个“大于零”不等式 0内部会转换为-LMI 0。我们通常全用正数手动处理正负号更清晰。q: 项在LMI中的块行和块列索引。通常设为0表示该项作用于整个矩阵1x1块。只有在处理分块矩阵LMI时才会用到非零值。x和y: 矩阵变量ID。x指定包含变量的矩阵y指定其转置如果涉及。具体规则见下。A,B: 系数矩阵。它们与变量或其转置相乘。如果项中不包含变量则A或B中有一个是标量常数通常是1。flag: 可选字符串默认为s对称项。当添加形如X*A或A*X的项时如果希望自动添加其转置项以保持整体对称性则使用s。对于常数项或明显不对称的项使用n非对称。x和y的取值规则重中之重如果项是常数矩阵如Q则x0, y0。此时A就是该常数矩阵或标量。如果项是变量矩阵如P则x为该变量的IDvar_idy0。此时A通常是左乘的系数矩阵若A1则表示变量本身。如果项是变量矩阵的转置如P*A中的P如果写在左边但实际是A*P的一部分则x0, y为该变量的ID。此时B是右乘的系数矩阵。如果项是两个变量相乘较少见则x和y分别为两个变量的ID。核心技巧拆解项与符号处理lmiterm一次只添加一个“项”。一个复杂的表达式需要拆成多个项。例如表达式A*P P*A包含两个项A*P和P*A。我们需要分别添加它们。所有项都必须放在不等式的一侧。例如对于 Lyapunov 不等式A*P P*A -Q我们需要将其改写为A*P P*A Q 0。那么在代码中我们就需要添加三项A*P,P*A,Q。实操心得二lmiterm的“分项累加”思维我习惯像搭积木一样构建LMI。首先在纸上或注释里把目标不等式整理成Sum(terms) 0的形式。然后对求和号里的每一项单独写一个lmiterm。每写一项都检查其termID是否指向正确的LMI编号以及x,y赋值是否符合上述规则。一个非常有效的调试方法是在构建完所有LMI后用lmiedit命令打开GUI查看器直观地检查你构建的矩阵是否正确。3. 完整案例连续系统状态反馈控制器设计与求解现在我们通过一个经典问题——连续系统状态反馈镇定——来串联整个流程。问题描述给定一个线性时不变系统dx/dt A*x B*u设计状态反馈控制律u K*x使得闭环系统dx/dt (AB*K)*x渐近稳定。这可以转化为一个LMI可行性问题。步骤1问题转化为LMI根据Lyapunov稳定性理论寻找一个对称正定矩阵P和一个矩阵K使得闭环系统满足(AB*K) * P P * (AB*K) 0这个不等式关于P和K不是线性的因为包含了P*B*K和K*B*P。为了将其线性化我们引入一个变量替换。令Y K * X其中X inv(P)。对上式左右同时乘以X合同变换并代入Y K*X可以得到一个关于X和Y的LMIA*X X*A B*Y Y*B 0同时P 0等价于X 0。 因此我们的LMI问题为 寻找对称矩阵X和矩阵Y使得A*X X*A B*Y Y*B 0X 0(即-X 0) 求解得到X和Y后控制器增益为K Y * inv(X)。步骤2MATLAB代码实现% 案例连续系统状态反馈控制器设计 clear; clc; % 1. 定义系统矩阵 (不稳定系统) A [1 2; -1 0]; B [0; 1]; n size(A, 1); % 状态维度 m size(B, 2); % 输入维度 % 2. 初始化LMI系统描述 setlmis([]); % 开始描述LMI系统 % 3. 定义矩阵变量 % X 是 n x n 对称正定矩阵 (type1, 满块) [X_id, nX, sX] lmivar(1, [n, 1]); % struct[n,1] 表示标量矩阵这里错了应该是[n,0] % 更正对于自由对称矩阵P应使用 [n, 0] [X_id, nX, sX] lmivar(1, [n, 0]); % Y 是 m x n 的矩形矩阵 (type2) [Y_id, nY, sY] lmivar(2, [m, n]); % 4. 定义第一个LMI: A*X X*A B*Y Y*B 0 lmi_id 1; % 第一个不等式编号 % 项1: A*X lmiterm([lmi_id, 1, 1, X_id], A, 1); % [1,1, X_id]: 添加到LMI1的(1,1)块项为 A * X * 1 % 注意这里 A 左乘 X所以 A 是系数放在第二个参数位置。1表示右乘单位阵。 % 项2: X*A lmiterm([lmi_id, 1, 1, X_id], 1, A, s); % 使用s标志自动添加对称部分 (X*A) (A*X) % 实际上因为项1已经加了A*X这里用s会自动补上其转置(A*X) X*A X*A (X对称)。 % 更稳妥的写法是分别添加两项但s在这里更简洁且不易出错。 % 项3: B*Y lmiterm([lmi_id, 1, 1, Y_id], B, 1); % B * Y * 1 % 项4: Y*B lmiterm([lmi_id, 1, 1, Y_id], 1, B, s); % 使用s自动补全 (B*Y) (Y*B) % 5. 定义第二个LMI: X 0 (等价于 -X 0) lmi_id 2; lmiterm([lmi_id, 1, 1, X_id], -1, 1); % 添加项-1 * X * 1 0 X 0 % 6. 完成LMI系统描述 lmisys getlmis; % 7. 求解LMI可行性问题 [tmin, xfeas] feasp(lmisys); % 8. 检查可行性并提取解 if tmin 0 % tmin是可行性测度小于0表示可行 disp(LMI可行); % 从解向量xfeas中提取矩阵变量 X_sol dec2mat(lmisys, xfeas, X_id); Y_sol dec2mat(lmisys, xfeas, Y_id); % 计算控制器增益 K Y * X^{-1} K Y_sol / X_sol; % 等价于 Y_sol * inv(X_sol) disp(求解得到的 X:); disp(X_sol); disp(求解得到的 Y:); disp(Y_sol); disp(状态反馈增益矩阵 K:); disp(K); % 验证闭环系统稳定性 A_cl A B*K; eig_cl eig(A_cl); disp(闭环系统特征值:); disp(eig_cl); if all(real(eig_cl) 0) disp(闭环系统稳定); else disp(警告闭环系统不稳定请检查求解结果。); end else disp(未找到可行解。tmin ); disp(tmin); end步骤3代码关键点解析setlmis([])与getlmis: 这是构建LMI问题的固定框架。setlmis([])初始化一个空的LMI系统之后所有的lmivar和lmiterm都在向这个系统添加内容。getlmis最终获取完整的LMI系统描述对象lmisys用于后续求解。feasp函数: 这是求解LMI可行性问题的核心求解器。它尝试找到一组决策变量使得所有LMI约束都被满足。返回值tmin是最小化约束违背的量tmin 0意味着找到了可行解所有LMI严格小于零。xfeas是决策变量的解向量。dec2mat函数: 求解器返回的解xfeas是一个向量。我们需要用dec2mat函数根据之前定义的变量结构lmisys将其还原成我们熟悉的矩阵形式。参数依次为LMI系统、解向量、变量ID。符号处理: 注意第二个LMIX 0我们是通过添加项-X到 0的不等式中实现的。这是LMI工具箱的标准做法将所有约束统一为“小于零”形式。4. 高级应用与性能优化技巧掌握了基础构建方法后我们可以处理更复杂的问题并优化求解过程。4.1 处理多个LMI约束与分块矩阵有时一个问题包含多个LMI或者一个LMI本身是分块矩阵。lmiterm的termID中第二、三个参数p, q就派上了用场。示例带有性能约束的H∞状态反馈问题设计u Kx使得闭环系统满足||Tzw||∞ γ。这可以转化为以下LMI组简化形式 寻找X 0,Y, 使得[ A*X X*A B*Y Y*B Bw (Cz*XDzu*Y) ] [ Bw -γ*I Dzw ] 0 [ Cz*XDzu*Y Dzw -γ*I ]以及X 0。 这是一个2x2的分块矩阵LMI实际上左上角块是另一个矩阵不等式。% ... 假设已定义 A, B, Bw, Cz, Dzu, Dzw, gamma ... setlmis([]); [X_id, ~, ~] lmivar(1, [n, 0]); [Y_id, ~, ~] lmivar(2, [m, n]); lmi_id 1; % H∞性能LMI % 块(1,1): A*X X*A B*Y Y*B lmiterm([lmi_id, 1, 1, X_id], A, 1); lmiterm([lmi_id, 1, 1, X_id], 1, A, s); lmiterm([lmi_id, 1, 1, Y_id], B, 1, s); % s 会自动添加 B*Y Y*B % 块(1,2): Bw lmiterm([lmi_id, 1, 2, 0], Bw); % 常数项x0,y0 % 块(1,3): (Cz*X Dzu*Y) % 先添加 Cz*X lmiterm([lmi_id, 3, 1, X_id], Cz, 1); % 注意这项在(3,1)位置是(1,3)的转置 % 再添加 Dzu*Y lmiterm([lmi_id, 3, 1, Y_id], Dzu, 1); % 块(2,1): Bw (是块(1,2)的转置由工具箱自动保证对称性通常只需定义上三角或下三角) % 因为我们用s或完整定义这里通常只需定义(1,2)对称部分会自动处理。但为清晰也可定义 % lmiterm([lmi_id, 2, 1, 0], Bw); % 块(2,2): -gamma*I lmiterm([lmi_id, 2, 2, 0], -gamma, eye(size(Bw,2))); % 块(2,3): Dzw lmiterm([lmi_id, 2, 3, 0], Dzw); % 块(3,1): (已定义) % 块(3,2): Dzw lmiterm([lmi_id, 3, 2, 0], Dzw); % 块(3,3): -gamma*I lmiterm([lmi_id, 3, 3, 0], -gamma, eye(size(Cz,1))); % 第二个LMI: X 0 lmi_id 2; lmiterm([lmi_id, 1, 1, X_id], -1, 1); % ... 后续求解步骤同上在定义分块矩阵时关键是理清每个块的位置(p, q)。工具箱不强制要求定义对称位置但为了清晰和避免遗漏建议至少定义完上三角部分并利用s标志处理对称项。4.2 优化问题求解mincx 与 gevp除了可行性问题feaspLMI工具箱还能解决两类优化问题线性目标最小化 (mincx)最小化c * x满足LMI约束。其中c是用户定义的向量x是决策变量向量。这常用于最小化线性组合的变量例如最小化矩阵的迹trace(X)。关键需要用defcx函数来定义目标函数中的系数向量c。defcx(lmisys, k, X_id, i, j)用于设置与变量X_id的(i,j)位置元素相关的系数。对于最小化迹通常是对角线元素系数设为1。% 在 getlmis 之后定义目标为最小化 trace(X) n decnbr(lmisys); % 获取决策变量数量 c zeros(n, 1); for j 1:n [Xj, ~, ~] defcx(lmisys, j, X_id, j, j); % 获取X_id第(j,j)个元素在决策向量中的索引 c(Xj) 1; % 设置该位置系数为1 end options [1e-2, 100, 1e5, 10, 0]; % 优化选项 [copt, xopt] mincx(lmisys, c, options);广义特征值问题 (gevp)最小化标量 λ使得满足 LMI1(x) 0 且 LMI2(x) λ * LMI3(x)。这常用于求解 H∞ 范数最小化 γ。调用格式[lopt, xopt] gevp(lmisys, nlfc, options)。nlfc是第一个LMI约束的个数即LMI1。实操心得三求解器选项与数值稳定性默认的求解器选项可能不适用于所有问题。feasp,mincx,gevp都可以通过options向量调整参数如初始步长、最大迭代次数、精度等。options [1e-2, 200, 1e8, 10, 0]; % 参数含义[精度, 最大迭代次数, 可行性半径, 迭代显示频率, 目标值] [tmin, xfeas] feasp(lmisys, options);如果求解失败tmin 0可以尝试放宽可行性半径增大options(3)。检查LMI问题的标度。如果矩阵A,B的元素数量级差异巨大如1e-9和1e3会导致数值问题。尽量对系统模型进行归一化处理。确保问题本身是可行的。有时需要调整性能指标 γ 或初始猜测。5. 调试技巧与常见问题排查即使思路正确构建LMI代码时也极易出错。以下是我在实践中总结的排查清单。5.1 问题排查速查表现象可能原因检查与解决方法feasp返回tmin 0(不可行)1. 问题本身无解。2. LMI构建错误符号、项遗漏。3. 数值问题矩阵病态。1. 检查问题转化过程是否正确理论是否保证有解如系统是否可控。2.使用lmiedit可视化检查。输入lmiedit(lmisys)可以图形化查看每个LMI的每个块核对每一项是否正确添加。3. 尝试更宽松的可行性半径options(3)。4. 对系统矩阵进行缩放。求解结果K使系统不稳定1. 提取变量ID与定义顺序不符。2.dec2mat参数错误。3. 控制器计算公式错误。1. 确保dec2mat中的变量ID与lmivar返回的X_id,Y_id一致。2. 验证闭环矩阵AB*K是否等于A B*(Y_sol/X_sol)。3. 检查Lyapunov不等式是否用求解得到的X_sol和K验证成立。错误Index exceeds matrix dimensions或lmivar/lmiterm参数错误1.lmivar的struct参数格式错误。2.lmiterm的termID中变量ID超出范围。3. 系数矩阵A,B维度不匹配。1. 仔细核对lmivar的type和struct文档。2. 确保lmiterm中引用的var_id是已定义的。3. 计算每一项的维度size(A,2)应与变量X的行数匹配等。求解速度非常慢1. 问题规模太大变量多。2. LMI条件数差。1. 尝试减少变量用type1的标量块[n,1]代替满块[n,0]如果结构允许。2. 使用solver选项尝试不同的内点法算法如options(5)。3. 考虑是否有更简洁的LMI表述。5.2 必备调试工具lmieditlmiedit(lmisys)是LMI工具箱中最强大的调试工具没有之一。它会打开一个GUI窗口以矩阵块的形式展示你定义的所有LMI。你可以清晰地看到每个LMI的维度。每个块是由哪些项组成的。每个项是常数、变量还是变量乘积。变量的位置和系数。当你的求解结果不符合预期时第一件事就是打开lmiedit逐项核对是否与你纸上推导的公式一致。这能解决90%的构建错误。5.3 数值验证用解回代验证LMI求解得到xfeas后不要仅仅相信tmin 0。应该将解出的矩阵变量代回原始LMI计算其最大特征值确保其为负。% 假设有两个LMI解已存储在 xfeas 中 X dec2mat(lmisys, xfeas, X_id); Y dec2mat(lmisys, xfeas, Y_id); % 验证第一个LMI: A*X X*A B*Y Y*B 0 LMI1_left A*X X*A B*Y Y*B; eig1 eig(LMI1_left); max_eig1 max(real(eig1)); fprintf(第一个LMI左端最大特征值: %.4e\\n, max_eig1); % 验证第二个LMI: X 0 (即最小特征值0) min_eig_X min(real(eig(X))); fprintf(矩阵X的最小特征值: %.4e\\n, min_eig_X); if max_eig1 0 min_eig_X 0 disp(所有LMI约束均严格满足。); else disp(警告解不严格满足LMI约束); end这种验证能给你十足的信心确认整个建模和求解流程是正确的。6. 从理论到实践扩展应用场景与思路LMI的应用远不止于控制器设计。一旦掌握了lmivar和lmiterm的思维你可以将其应用到众多领域。1. 观测器滤波器设计对于系统dx/dt A*x B*u, y C*x设计全维状态观测器d\hat{x}/dt A*\hat{x} B*u L*(y-C*\hat{x})。误差动态系统为de/dt (A-L*C)*e。通过类似的变量替换W P*L可以转化为关于P和W的LMIA*P P*A - W*C - C*W 0, P0进而解得L P^{-1}*W。2. 饱和控制系统分析与设计考虑执行器饱和sat(u)。利用扇形条件或多面体描述可以将饱和非线性转化为一组线性微分包含进而用LMI来估计吸引域或设计抗饱和补偿器。3. 时滞系统稳定性分析对于含有时滞的系统利用Lyapunov-Krasovskii泛函或Lyapunov-Razumikhin函数可以将稳定性条件转化为包含时滞项的LMI通常需要引入一些松弛矩阵变量来降低保守性。4. 组合优化中的半定规划SDP图论中的最大割问题、传感器网络定位等可以建模为半定规划这正是LMI工具箱 (mincx) 可以求解的问题类型。个人体会与最后建议LMI工具箱的学习曲线前期比较陡峭核心难点在于从连续的数学等式/不等式到离散的、声明式的代码指令的思维转换。我的经验是永远从最简单的、有解析解的例子开始。例如对于稳定的矩阵ALyapunov方程A*PP*A-Q的解P一定是正定的。你可以先用lyap函数求解再用LMI工具箱去求解A*PP*A 0, P0对比结果是否一致。这种“双轨验证”是建立代码信心的最佳方式。另外不要畏惧lmiedit。在初期每写几行lmiterm就运行一次getlmis和lmiedit查看中间结果远比一次性写几十行代码然后面对一堆错误要高效得多。当你熟练之后构建一个中等复杂度的LMI问题就像搭积木一样自然这种将复杂系统问题转化为可计算模型的能力会让你在研究和工程中受益匪浅。