简介：《LMI工具箱手册》是一本关于线性矩阵不等式（LMI）技术的中英文双语专业参考资料，系统介绍LMI的基本理论、建模方法及在控制理论、优化和信号处理等领域的应用。手册详细讲解LMI工具箱在MATLAB环境中的使用，涵盖问题设定、求解函数（如lmiroot）、结果可视化（如lmiplot）以及与Simulink的联合仿真。内容涉及H_∞控制、鲁棒控制、状态/输出反馈设计、滤波器设计等核心应用场景，适合研究人员与工程师深入学习与实践。通过本手册，用户可快速掌握LMI建模与求解流程，提升复杂系统分析与设计能力。
线性矩阵不等式（LMI）：从理论到MATLAB实战的深度解析
在现代控制系统的工程实践中，我们常常面临一个核心挑战：如何在不依赖显式求解微分方程的前提下，判断一个复杂动态系统是否稳定？又该如何设计控制器，使其不仅能在理想条件下工作，还能抵御参数摄动、外部干扰甚至模型不确定性的影响？
传统方法如根轨迹法、Nyquist判据虽然直观，但在处理高维、非线性或含不确定性的系统时显得力不从心。而 线性矩阵不等式 （Linear Matrix Inequality, LMI）的出现，彻底改变了这一局面。🚀 它将原本抽象的稳定性分析和性能优化问题，转化为一组可计算、可编程且具备良好数值特性的凸优化问题——这不仅是数学工具的进步，更是控制工程师手中的一把“万能钥匙”。
想象一下，你不再需要手动猜测Lyapunov函数的形式，也不必为双线性项（BMI）的非凸性头疼；取而代之的是，只需定义几个矩阵变量，写下几条结构化的约束条件，剩下的搜索与验证过程全部交给高效的内点法求解器自动完成。这种“声明式建模 + 自动化求解”的范式，正是LMI的魅力所在。
本文将带你深入LMI的世界，从最基础的数学结构出发，逐步揭开它在稳定性分析、H₂/H∞控制、鲁棒设计中的强大能力，并结合MATLAB LMI工具箱的实际编码技巧，构建一条从理论推导到工程落地的完整路径。准备好了吗？让我们开始这场穿越矩阵世界的旅程吧！✨
什么是LMI？它的数学本质是什么？
简单来说， 线性矩阵不等式 就是形如：
$$
F(x) = F_0 + \sum_{i=1}^n x_i F_i \succ 0
$$
的矩阵不等式，其中 $ x \in \mathbb{R}^n $ 是决策变量向量，$ F_i $ 是对称常数矩阵，符号 $ \succ 0 $ 表示矩阵正定。
乍一看，这似乎只是一个关于变量 $ x $ 的仿射组合，但它的背后隐藏着深刻的几何意义： 可行解集是一个凸集 。这意味着如果你能找到两个满足条件的解，那么它们之间的任意线性插值也一定满足条件。这个“凸性”是LMI之所以强大的根本原因——它保证了任何局部最优解都是全局最优解，避免了陷入复杂的非凸陷阱。
举个经典的例子：考虑连续时间系统的Lyapunov稳定性条件：
$$
A^T P + P A \prec 0, \quad P \succ 0
$$
这是一个关于未知矩阵 $ P $ 的线性矩阵不等式。尽管 $ P $ 出现在乘积项中，但由于它是唯一的变量，整个表达式对 $ P $ 是线性的（即仿射）。因此，我们可以把它直接送入LMI求解器，让算法自动寻找是否存在这样一个正定矩阵 $ P $ 来证明系统稳定。
下面这段MATLAB代码展示了如何使用 lmivar 和 lmiterm 来构建一个简单的LMI：要求矩阵 $ P > I $。
% 示例：定义一个简单的LMI —— P > I
setlmis([]);
P = lmivar(1, [2,1]); % 定义对称矩阵变量 P ∈ S²⁺
lmiterm([1 1 1 P], 1, 1); % LMI: P
lmiterm([-1 1 1 0], 1, 1); % > I
LMISYS = getlmis;
这里的关键在于理解 lmiterm 的索引规则：
- [1 1 1 P] 表示这是第1个LMI的(1,1)位置，涉及变量P；
- [-1 ...] 中的负号表示该LMI应小于零，所以 -I < 0 等价于 P > I 。
是不是有点像搭积木？每一块 lmiterm 都在拼接最终的矩阵不等式结构。🧱
LMI如何重塑控制系统的设计逻辑？
如果说经典控制理论是“基于试错的艺术”，那么基于LMI的方法更像是“基于推理的科学”。它允许我们将多种性能指标统一在一个框架下进行联合优化。下面我们通过几个典型场景来看看LMI是如何工作的。
🔄 连续系统的稳定性验证：告别手算Lyapunov函数
对于一个二阶系统：
$$
\dot{x}(t) = A x(t), \quad A = \begin{bmatrix} -1 & 2 \ 0 & -3 \end{bmatrix}
$$
我们知道其稳定的充要条件是存在正定矩阵 $ P $ 满足：
$$
A^T P + P A \prec 0
$$
传统做法是假设 $ P = \begin{bmatrix} p_{11} & p_{12} \ p_{12} & p_{22} \end{bmatrix} $，然后展开不等式并尝试找到合适的参数。但这种方式效率低、易出错，尤其当维度升高时几乎不可行。
而用LMI方法，我们可以完全自动化这一过程：
% MATLAB代码：连续系统Lyapunov LMI求解
setlmis([]);
P = lmivar(1, [2,1]); % 定义对称矩阵变量P ∈ S²⁺
% 添加LMI项：A'P + PA < 0
lmiterm([1 1 1 P], 1, A', 's'); % s表示对称项 A'*P + P*A
lmiterm([-1 1 1 0], -1e-6, eye(2)); % 添加 -εI 确保严格小于零
% 构建LMI系统并求解可行性
lmisys = getlmis();
[tmin, feas] = feasp(lmisys);
if tmin < 0
disp('系统稳定，找到可行P');
P_sol = dec2mat(lmisys, feas, P);
disp('Lyapunov矩阵P:');
disp(P_sol);
else
disp('未找到可行解，系统可能不稳定');
end
💡 小贴士：为什么要在右边加上 -1e-6 * eye(2) ？
因为数值求解器很难精确判断“≤ 0”这样的半严格不等式。加入一个小扰动 $-\varepsilon I$ 可以确保我们得到的是严格的负定解，从而提升数值稳定性。
整个流程可以用Mermaid清晰地展现出来：
graph TD
A[开始] --> B[定义系统矩阵A]
B --> C[初始化LMI系统 setlmis]
C --> D[声明矩阵变量P = lmivar]
D --> E[构建LMI: A^T P + P A < 0]
E --> F[调用 lmiterm 添加各项]
F --> G[生成LMI系统 getlmis]
G --> H[调用 feasp 求解可行性]
H --> I{是否存在可行解?}
I -- 是 --> J[输出P，判定稳定]
I -- 否 --> K[系统可能不稳定]
J --> L[结束]
K --> L
这种模块化、可复用的建模方式，使得即使面对几十维的系统，我们也只需关注变量定义和约束构造，无需担心底层实现细节。
🔁 离散系统的稳定性：Schur补的妙用
离散时间系统的稳定性条件略有不同：
$$
x(k+1) = A_d x(k)
$$
渐近稳定的充要条件是存在 $ P \succ 0 $ 使得：
$$
A_d^T P A_d - P \prec 0
$$
虽然形式上是非线性的（出现了 $ A_d^T P A_d $），但我们可以通过 Schur补引理 将其转换为标准LMI：
如果 $ Z > 0 $，则 $ X - Y^T Z^{-1} Y < 0 $ 当且仅当
$$
\begin{bmatrix}
-X & Y^T \
Y & -Z
\end{bmatrix} < 0
$$
应用到这里，令 $ X = P $, $ Y = P A_d $, $ Z = P $，于是有：
$$
\begin{bmatrix}
P & A_d^T P \
P A_d & P
\end{bmatrix} > 0
$$
这就变成了一个关于 $ P $ 的线性矩阵不等式！
来看具体实现：
% 离散系统LMI稳定性验证
setlmis([]);
P = lmivar(1, [2,1]);
% 第一个LMI: P > 0
lmiterm([-2 1 1 P], 1, 1); % P > 0
lmiterm([-2 1 1 0], -1e-6, eye(2));% 小扰动保证严格正定
% 第二个LMI: [P, A'P; PA, P] > 0
lmiterm([3 1 1 P], 1, 1); % 左上角 P
lmiterm([3 1 2 P], 1, A_d'); % 上右 A_d' * P
lmiterm([3 2 2 P], -1, 1); % 右下 -P，注意符号反转
lmiterm([3 2 1 0], A_d); % 下左 P*A_d，常数项
lmisys = getlmis();
[tmin, feas] = feasp(lmisys);
if tmin < 0
disp('离散系统稳定');
P_sol = dec2mat(lmisys, feas, P);
eigvals = eig(P_sol);
fprintf('P的特征值: %.4f, %.4f\n', eigvals(1), eigvals(2));
else
disp('系统不稳定或无可行解');
end
你会发现，在 lmiterm 中右下角用了 -1 作为系数，这是因为我们要构造的是 大于零 的矩阵，而LMI工具箱默认处理的是 小于零 的形式，所以需要取反。
| 参数 | 含义 | 类型 |
|---|---|---|
A | 系统状态矩阵 | 数值矩阵 (n×n) |
P | Lyapunov矩阵变量 | 对称正定矩阵变量 |
lmivar | 定义矩阵变量 | 函数句柄 |
lmiterm | 插入LMI项 | 函数句柄 |
feasp | 可行性求解器 | LMI工具箱函数 |
⚠️ 面对不确定性：鲁棒稳定性不再是黑箱
现实中，系统参数往往存在摄动。比如某个电机的转动惯量可能因温度变化产生±5%的偏差。如果我们只针对标称模型设计控制器，一旦实际参数偏离，系统就可能失稳。
这时就需要引入 鲁棒稳定性分析 。假设系统矩阵具有范数有界不确定性：
$$
A = A_0 + \Delta A, \quad |\Delta A| \leq \delta
$$
更进一步，若 $ \Delta A = D F E $，其中 $ F^T F \leq I $，则可以通过S-procedure得到保守但可解的LMI条件：
存在 $ P \succ 0 $, $ \lambda > 0 $，使得：
$$
\begin{bmatrix}
A_0^T P + P A_0 + \lambda^{-1} E^T E & P D \
D^T P & -\lambda I
\end{bmatrix} < 0
$$
这个条件确保无论 $ F $ 如何变化（只要满足 $ F^T F \leq I $），闭环系统都保持稳定。👏
来看MATLAB实现：
% 鲁棒稳定性LMI建模
setlmis([]);
P = lmivar(1, [2,1]);
lam = lmivar(1, [1,0]); % 标量变量 λ > 0
% 构造复合LMI矩阵
% 第一块行第一列: A0'*P + P*A0 + (1/lam)*E'*E
lmiterm([1 1 1 P], 1, A0', 's');
lmiterm([1 1 1 lam], -1, 1); % -λ^{-1} via Schur补技巧
lmiterm([1 1 1 0], E'*E);
% 第一块行第二列: P*D
lmiterm([1 1 2 P], 1, D);
lmiterm([1 2 1 P], D', 1);
% 第二块行第二列: -λ*I
lmiterm([1 2 2 0], -eye(size(D,2)));
lmisys = getlmis();
[tmin, feas] = feasp(lmisys);
if tmin < 0
disp('系统对所有允许摄动均鲁棒稳定');
else
disp('无法证明鲁棒稳定性');
end
注意到我们用了一个小技巧： lmiterm([1 1 1 lam], -1, 1) 实际上是在利用Schur补处理 $ \lambda^{-1} $ 项。这是LMI建模中常见的“变通手法”——通过引入辅助变量绕过非线性运算。
性能指标也能放进LMI？H₂、H∞、极点配置全搞定！
稳定性只是第一步，真正的控制系统还需要满足各种性能需求。幸运的是，这些看似复杂的性能边界也可以被优雅地嵌入LMI框架。
📊 H₂性能：最小化平均能量增益
H₂范数衡量系统从白噪声输入到输出的平均能量响应。对于系统：
$$
\dot{x} = Ax + Bw, \quad z = Cx
$$
其H₂性能可通过以下LMI最小化问题求解：
$$
\min \gamma \quad \text{s.t.} \quad
\begin{bmatrix}
-P & P \
P & -Q
\end{bmatrix} < 0,\quad
\begin{bmatrix}
-Q & C^T \
C & -\gamma I
\end{bmatrix} < 0,\quad A^T P + P A + B B^T < 0
$$
简化后常用形式为：
存在 $ P>0, Q>0 $，使得：
$$
\begin{bmatrix}
A^T P + P A + C^T C & P B \
B^T P & -I
\end{bmatrix} < 0, \quad \text{tr}(Q) < \gamma
$$
MATLAB中可通过 mincx 实现目标函数最小化：
% H2性能优化示例
A = [-1 1; 0 -2]; B = [1; 1]; C = [1 0];
setlmis([]);
P = lmivar(1, [2,1]);
[gam2, X] = defcx(lmisys, P, 'trace'); % 定义trace(P)为目标
% LMI 1: A'P + PA + C'C < 0
lmiterm([1 1 1 P], 1, A', 's');
lmiterm([1 1 1 0], C'*C);
lmiterm([1 1 2 0], P*B);
lmiterm([1 2 1 0], B'*P);
lmiterm([1 2 2 -1], 1, 1); % -I
lmisys = getlmis();
% 最小化 trace(P)
c = defcx(lmisys, P, 'trace');
[copt, xopt] = mincx(lmisys, c);
fprintf('最优H2性能界: %.4f\n', sqrt(copt));
| 指标 | 物理意义 | LMI表达方式 |
|---|---|---|
| H₂ | 平均能量增益 | 迹最小化 + Lyapunov约束 |
| H∞ | 峰值增益 | Riccati型LMI或ARE松弛 |
| 极点区域 | 动态响应速度 | 区域映射LMI（D-稳定） |
🛡️ H∞性能：抑制最坏情况干扰
H∞控制的目标是使外部扰动对输出的影响尽可能小，即：
$$
|T_{zw}(s)|_\infty < \gamma
$$
根据 Bounded Real Lemma ，其充分条件为存在 $ P \succ 0 $ 使得：
$$
\begin{bmatrix}
A^T P + P A + C_z^T C_z & P B_w \
B_w^T P & -\gamma I
\end{bmatrix} < 0
$$
代码实现如下：
% H-infinity性能验证
gam = 1.0;
setlmis([]);
P = lmivar(1, [2,1]);
lmiterm([1 1 1 P], 1, A', 's');
lmiterm([1 1 1 0], Cz'*Cz);
lmiterm([1 1 2 0], P*Bw);
lmiterm([1 2 1 0], Bw'*P);
lmiterm([1 2 2 0], -gam*eye(size(Bw,2)));
lmisys = getlmis();
[tmin, feas] = feasp(lmisys);
如果求解成功，说明系统具有小于 $ \gamma $ 的H∞增益，具备良好的抗干扰能力。
🎯 D-稳定性：你想让极点待在哪？
有时候我们不仅希望系统稳定，还希望闭环极点落在特定区域内，比如左半平面某个圆盘内，以保证足够的阻尼比和响应速度。
例如，希望所有极点位于圆心 $ (-α, 0) $、半径 $ r $ 的圆内，则满足：
$$
\alpha I + r A^T P + r P A < 0
$$
这类条件可以与其他性能指标联合使用，形成多目标综合设计。🎉
如何搭建你的LMI开发环境？MATLAB配置全指南
再强大的理论也需要可靠的工具支持。目前， MATLAB的LMI Control Toolbox 是最成熟、最广泛使用的LMI求解平台。下面我们一步步教你如何配置一个高效的工作环境。
✅ 检查版本兼容性
首先确认你的MATLAB版本是否支持LMI工具箱。主流支持范围是 R2010a 至 R2023b 。自 R2020b 起，部分GUI功能（如 lmiedit ）已被弃用，推荐使用脚本化编程。
| MATLAB版本 | LMI工具箱状态 | 推荐使用场景 |
|---|---|---|
| R2010a - R2018b | 完全支持，功能完整 | 教学演示、传统项目维护 |
| R2019a - R2021b | 支持但提示弃用警告 | 过渡期项目，建议重构代码 |
| R2022a 及以上 | 新架构支持，仅限程序化建模 | 新项目开发、高性能计算 |
| R2024+ (预发布) | 不再支持GUI编辑器 | 需采用API方式建模 |
✅ 强烈建议使用 R2021b 或更高版本 ，以获得最新的求解器性能优化（如对Mosek、SDPT3的支持）。
🔧 安装与路径设置
% 检查是否已安装Control System Toolbox和LMI模块
if ~license('test', 'control_toolbox')
error('Control System Toolbox未安装，请先通过Add-ons安装');
end
if ~exist('lmivar', 'file')
warning('LMI Control Toolbox未检测到，请前往Home > Add-Ons > Get Add-Ons搜索"LMI Control Toolbox"进行安装');
else
disp('LMI工具箱已成功加载');
end
安装完成后，运行以下命令验证路径：
which lmivar
% 正常输出应类似：
% C:\Program Files\MATLAB\R2023b\toolbox\robust\lmi\lmivar.m
若路径错误，手动添加：
addpath(genpath(fullfile(matlabroot, 'toolbox', 'robust', 'lmi')));
savepath; % 永久保存路径更改
🧪 快速测试：验证工具箱是否正常工作
% 测试案例：验证连续系统A = [-1, 2; 0, -3] 的渐近稳定性
A = [-1, 2; 0, -3];
setlmis([]);
P = lmivar(1, [2 1]);
lmiterm([1 1 1 P], 1, A', 's');
lmiterm([1 1 1 0], -1e-6);
LMIs = getlmis();
[tmin, feas] = feasp(LMIs);
if tmin < 0 && feas ~= []
disp('✅ LMI工具箱工作正常：系统稳定，P存在正定解');
else
disp('❌ LMI工具箱存在问题：无法找到可行解');
end
这段代码不仅能检验安装完整性，也展示了基本建模流程，是部署新环境后的必备验证环节。
LMI建模的艺术：三大转化技巧揭秘
许多原始控制问题本质上是非凸的，尤其是当控制器增益直接作为变量时，会导致双线性矩阵不等式（BMI）。解决之道在于巧妙的变量重构。以下是三种核心技巧：
🔍 技巧一：Schur补引理的应用
Schur补是LMI建模中最常用的“魔法咒语”。它可以将形如 $ X - Y^T Z^{-1} Y < 0 $ 的非线性不等式转化为：
$$
\begin{bmatrix}
-X & Y^T \
Y & -Z
\end{bmatrix} < 0
$$
这在H₂性能、滤波器设计中极为常见。
🔄 技巧二：变量替换消除双线性项
以状态反馈为例，闭环稳定性条件：
$$
(A + BK)^T P + P(A + BK) < 0
$$
关于 $ K $ 和 $ P $ 是双线性的。引入新变量 $ Y = KP $，则变为：
$$
A^T P + P A + Y^T B^T + B Y < 0
$$
这是一个关于 $ P $ 和 $ Y $ 的线性不等式！恢复控制律时只需计算 $ K = Y P^{-1} $。
% 状态反馈 LMI 廢模（变量替换）
n = size(A,1); m = size(B,2);
setlmis([]);
P = lmivar(1, [n,1]); % P: n×n 对称
Y = lmivar(2, [m,n]); % Y = K*P, m×n 矩阵变量
% 构造 LMI: A'P + PA + Y'B' + BY < 0
lmiterm([1 1 1 P], 1, A, 's'); % A'P + PA
lmiterm([1 1 1 Y], 1, B', 's'); % Y'B' + B Y
% P > 0
lmiterm([-1 1 1 P], 1, 1); % -P < 0 ⇒ P > 0
% 求解可行性
lmisys = getlmis();
[tmin, feas] = feasp(lmisys);
if tmin < 0
P_val = dec2mat(lmisys, feas, P);
Y_val = dec2mat(lmisys, feas, Y);
K = Y_val / P_val; % K = Y * inv(P)
end
🧱 技巧三：分块矩阵提升结构性
在MIMO或多子系统耦合问题中，使用分块矩阵表示可以清晰展现系统结构，便于设计分散控制器或分析交互影响。
graph TD
A[原始控制系统] --> B[建立状态空间模型]
B --> C{是否存在非线性项?}
C -->|是| D[应用Schur补或变量替换]
C -->|否| E[直接构造LMI]
D --> F[引入松弛变量或新变量Y]
F --> G[重构为仿射矩阵函数]
G --> H[使用lmivar/lmiterm编码]
H --> I[LMI求解器判断可行性]
I --> J{可行?}
J -->|是| K[提取控制器增益K=Y*inv(P)]
J -->|否| L[调整权重或放宽约束]
实战演练：H∞控制器设计与Simulink联合仿真
最后，让我们走完从建模到验证的完整流程。
🛰️ 场景设定：飞行器纵向动力学模型
假设某飞行器名义系统为：
$$
A = \begin{bmatrix}-0.8&1.0\-1.5&-2.0\end{bmatrix},\
B_u = \begin{bmatrix}0.5\1.0\end{bmatrix},\
B_w = \begin{bmatrix}0.1\0.2\end{bmatrix},\
C_z = \begin{bmatrix}1&0\end{bmatrix}
$$
存在不确定性 $ \Delta A = H\Delta(t)E $，目标是最小化干扰增益 $ |T_{zw}|_\infty < \gamma $。
💻 MATLAB实现
setlmis([]);
n = 2; m = 1;
Q = lmivar(1, [n, 1]); % Q = P^{-1}
Y = lmivar(2, [m, n]); % Y = K*Q
gam = lmivar(1, [1, 1]); % γ²
% Term: AQ + QA' + Bu*Y + Y'*Bu'
lmiterm([1 1 1 Q], A, 1, 's');
lmiterm([1 1 1 Y], B_u, 1);
lmiterm([1 1 1 Y], 1, B_u', 2);
% Cross terms with Bw and weighting
lmiterm([1 1 2 Q], 1, B_w);
lmiterm([1 2 2 0], -eye(1));
% Output term
lmiterm([1 1 3 Q], C_z, 1);
lmiterm([1 3 3 0], -eye(1));
% Minimize gamma
c = mat2dec(lmisys, eye(n), zeros(m,n), 1);
options = [1e-5, 0, 0, 0, 0];
[copt, xopt] = mincx(lmisys, c, options);
Qsol = dec2mat(lmisys, xopt, Q);
Ysol = dec2mat(lmisys, xopt, Y);
K = Ysol / Qsol;
gamma_min = sqrt(dec2mat(lmisys, xopt, gam));
🧪 Simulink验证
将控制器增益导入Simulink：
save('controller_gain.mat', 'K');
在Simulink中使用“Gain”模块实现 $ u = -Kx $，配置仿真环境如下：
| 参数 | 数值/类型 |
|---|---|
| 求解器 | ode45 |
| 仿真时间 | 10 s |
| 扰动信号类型 | 白噪声 + 阶跃 |
| 输出观测变量 | $ x_1, x_2, z $ |
| 采样周期 | 0.01 s |
| 控制器更新方式 | 连续时间状态反馈 |
运行后绘制李雅普诺夫函数曲线验证单调递减：
V = sum((x_data * P) .* x_data, 2);
plot(t_data, V); grid on;
ylabel('V(x)'); xlabel('Time (s)');
title('Lyapunov Function Evolution under Fault Condition');
结语：LMI不只是工具，更是一种思维方式
当你掌握了LMI，你就不再只是一个“调参工程师”，而是成为了一名能够系统化思考、形式化表达、自动化验证的 控制架构师 。🧠
它教会我们：
🔹 复杂问题可以分解为结构性的约束；
🔹 非凸难题可以通过变量重构迎刃而解；
🔹 理论与实践之间，有一座名为“数值计算”的桥梁。
未来，随着机器学习与控制的深度融合，LMI的思想将继续延伸至数据驱动控制、安全验证、神经网络稳定性等领域。掌握它，意味着你站在了智能系统设计的前沿阵地。🚀
所以，别再犹豫了——打开MATLAB，写下你的第一个 setlmis([]) ，开启属于你的LMI之旅吧！💫
简介：《LMI工具箱手册》是一本关于线性矩阵不等式（LMI）技术的中英文双语专业参考资料，系统介绍LMI的基本理论、建模方法及在控制理论、优化和信号处理等领域的应用。手册详细讲解LMI工具箱在MATLAB环境中的使用，涵盖问题设定、求解函数（如lmiroot）、结果可视化（如lmiplot）以及与Simulink的联合仿真。内容涉及H_∞控制、鲁棒控制、状态/输出反馈设计、滤波器设计等核心应用场景，适合研究人员与工程师深入学习与实践。通过本手册，用户可快速掌握LMI建模与求解流程，提升复杂系统分析与设计能力。
扫一扫
3642
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
