简介:围绕上海交通大学仪器科学与技术研究生课程“线性系统分析与设计”的大作业,提供了一份完整的MATLAB实现与报告参考。该作业面向需要夯实线性系统理论并掌握控制仿真工具的研究生和工程师,覆盖系统建模、传递函数推导、稳定性判据应用、控制器校正设计以及Simulink动态仿真等关键环节,既可用于课程项目参考,也能作为控制类课题的入门练习。压缩包包含3个文件,以2个m脚本和1个docx说明文档组成,整体大小仅193KB。两个m脚本分别对应不同分析模块,实现系统构建与仿真计算;docx文档则整理推导步骤、结果图表和设计结论,便于对照代码理解完整流程。目前已有680人学习该资源,适合正在准备线性系统课程大作业或希望快速上手MATLAB控制工具箱的读者。通过研究代码结构、仿真逻辑和报告写法,可以节省大量调试时间,同时获得从理论建模到工程实现的全链路思路。
1. 线性系统大作业到底在卷什么:从仿真发散到一次性过审的完整路径
凌晨两点,仿真曲线第三次冲向无穷大,你以为是自己运气差,其实只是状态反馈里某个符号写错了。线性系统大作业的核心套路其实非常固定:给一个物理对象写出状态空间模型,判断能控性与能观性,设计状态反馈与状态观测器,最后用仿真验证闭环性能并写成报告。表面上是“设计一个控制器”,实际考查的是你对状态空间、极点配置、观测器对偶原理的完整把握。下面就把这条路完整走一遍,把参数怎么定、仿真怎么搭、报告里哪些地方最容易被扣分,一个一个说清楚。适合正在写这门课大作业、准备复试项目或想补控制基础的从业者。
2. 从物理对象到状态空间:建模与能控能观性检查
2.1 先写出最小实现:机理建模的两种路线与四矩阵
常见的大作业对象无非那么几类:倒立摆、双水箱、直流电机、质量弹簧阻尼系统。我一般会优先选质量弹簧阻尼,因为它的物理直觉清楚、矩阵写出来不绕,出问题也好查。你看大部分评分点其实不在对象多复杂,而在状态空间写没写对、后面的设计步骤是不是完整闭环。
建模有两种路线。机理建模是直接对对象受力分析,写出微分方程再转状态空间;实验辨识则是给输入输出数据,用系统辨识工具箱拟合一阶或二阶模型。做课程大作业,机理建模是主流,因为老师要看到你从物理概念到数学描述的推导过程。辨识那条路留给实际工程里没法列方程的对象。
以质量弹簧阻尼系统为例,运动方程是:
m·ẍ + c·ẋ + k·x = F定义状态变量为位移和速度,即 x₁ = x、x₂ = ẋ,那么一阶微分方程组就是:
ẋ₁ = x₂ ẋ₂ = -(k/m)·x₁ - (c/m)·x₂ + (1/m)·F这一步是整个大作业的地基。很多人在后续控制器设计时发现算出来的 K 阵符号不对,回头一查,往往就是这里某个状态变量的正方向定义和题目不一致。我习惯把所有符号约定写在脚本注释里,避免做着做着把自己绕进去。
对应到 MATLAB 代码,最小实现是这样写的:
% 对象参数:质量 m(kg)、阻尼 c(N·s/m)、刚度 k(N/m) m = 2.0; c = 0.8; k = 5.0; % 状态空间四矩阵:x = [位移; 速度] A = [0 1; -k/m -c/m]; B = [0; 1/m]; C = [1 0]; % 输出位移 D = 0; % 建立连续时间系统模型 sys = ss(A, B, C, D);这段代码的逻辑很直白:A 矩阵的第一行对应“位移的导数等于速度”,第二行对应“速度的导数等于加速度,即合力除以质量”。B 矩阵只有第二行非零,因为外力 F 直接改变的是加速度,不直接改变位移。C 矩阵取第一列,表示我们传感器只能测到位移。D 为零,因为没有输入到输出的直通路径。
参数层面的讲究在量的匹配。m 变大,系统自然频率变低、响应变慢;c 变大,衰减变快;k 变大,振荡频率变高。如果题目给的参数单位是厘米而不是米,你又不统一,后面算出的极点、增益全部是错的,而且这种错非常隐蔽,仿真曲线看起来“有振荡、有收敛”,但数值对不上评分标准。
2.2 能控性与能观性:用秩判断,但别只信秩
能控性说的是系统里所有模态能不能被输入信号“够到”,能观性说的是所有模态能不能从输出信号“看穿”。这两个判断直接决定你后面能不能用状态反馈、能不能设计观测器。判断工具就是能控性矩阵 Qc 和能观性矩阵 Qo 的秩。
% 能控性矩阵与秩:rank = n 表示完全能控 Qc = ctrb(A, B); n = size(A, 1); fprintf('能控性矩阵秩 = %d (系统阶次 n = %d)\n', rank(Qc), n); % 能观性矩阵与秩:rank = n 表示完全能观 Qo = obsv(A, C); fprintf('能观性矩阵秩 = %d\n', rank(Qo));这段代码的要点是,rank 的结果必须等于系统阶次 n,才允许你做完整的极点配置和观测器设计。对二阶系统,Qc 是 2×2 方阵,很多人会顺手写 det(Qc),用行列式是否为 0 判断,但这里有个明显的坑:行列式对数值尺度极其敏感,当参数量级差异大时,明明能控的系统也可能算出接近 0 的行列式,把你带进错误结论。rank 背后基于奇异值分解,比 det 稳健得多。
如果发现秩不够,先别急着改代码。检查两件事:第一,输入 B 矩阵是否真的作用到了所有模态上,比如外力只加在位移方程里而速度模态完全不受影响,那就不可控;第二,模型是否是非最小实现,即存在可以被零极点对消掉的模态,这种情况要先做 minreal(sys) 再重新取四矩阵。能控能观性分析在报告里必须出现,因为它直接决定了你后面设计方法的合法性,少了这一段,评审老师一眼就能看出来你对状态空间的理解不到位。
3. 状态反馈与极点配置:把 K 阵算出来并当场验证
3.1 极点配置的两种落地写法:acker 与 place 怎么选
状态反馈的闭环形式是 u = -Kx,代入原系统后闭环矩阵变成 A - BK,闭环极点就是它的特征值。所谓极点配置,就是找一组 K,让这些特征值落在你希望的位置上。这里有两个现成函数:acker 和 place。单输入单输出系统,两个都能用;多输入系统,place 更稳定。
% 性能指标:自然频率 wn、阻尼比 zeta wn = 2.0; zeta = 0.5; % 主导极点:一对共轭复数极点 p_dom = roots([1 2*zeta*wn wn^2]); % 远极点:放在实轴更远处,保证不影响主导动态 p_far = -10 * wn; % 期望极点向量:注意系统是 3 阶,需要 3 个极点 p = [real(p_dom(1)) + 1i*imag(p_dom(1)); ... real(p_dom(2)) - 1i*imag(p_dom(2)); ... p_far]; % 单输入用 acker 直观,多输入或数值敏感的用 place K = place(A, B, p); % 验证:闭环极点与期望极点的最大偏差 cl_poles = eig(A - B*K); fprintf('最大极点偏差 = %e\n', max(abs(cl_poles - p)));这段代码的逻辑是:把性能指标转成极点位置,再让 place 反解出 K。对数域,主导极点决定超调量和调节时间,远极点负责处理未建模动态和执行器带宽。验证那一步非常关键,eig(A - B*K) 必须非常接近你设定的 p,但凡偏差在 1e-6 以上,就要怀疑是 A、B 维数不匹配,或者期望极点数不等于系统阶数。
参数上,zeta 取 0.4~0.7 对应超调量约 10%~25%,工程上常用 0.5 左右。wn 的选取要参考对象本身的响应速度,拿质量弹簧阻尼系统来说,wn 远大于 sqrt(k/m) 时,执行器需要输出的力会超出物理限制。这一点在后面仿真里会直接表现为控制量曲线封顶削波。还有个容易被忽视的细节:place 对期望极点的顺序不敏感,但如果有两个完全相同的极点,部分版本会给出 NaN。遇到这种情况,可以把重复极点拆成两个靠得很近但不相等的实极点。
3.2 极点位置不是越大越好:控制量与执行器带宽的边界
很多新手以为极点越往左移,系统收敛越快,就一定能拿高分。这是大作业里最容易翻车的心态。极点左移的本质是提高闭环带宽,带宽一高,K 阵数值变大,控制量 u = -Kx 的峰值会跟着涨。执行器一旦饱和,系统就从线性变成非线性,前面算的极点全部失效。
性能指标与极点参数的对应关系,可以直接套这两个公式:
| 指标 | 公式 | 说明 |
|---|---|---|
| 超调量 MP | exp(-ζπ / √(1-ζ²)) | ζ 越小,超调越大 |
| 调节时间 ts | 4 / (ζ·ωn) | ζ·ωn 越大,收敛越快 |
设计时先定 ζ,再定 ωn,最后才生成极点。如果只盯着收敛速度把 ωn 抬高,仿真画出来的控制量幅度会让你吓一跳。我一般会在脚本里顺手算一下控制量的峰值:
% 闭环系统与零输入响应 sys_cl = ss(A - B*K, B, C, D); t = 0:0.01:10; x0 = [1; 0; 0]; % 初值:位移 1,速度 0,附加状态 0 [y, t, x] = lsim(sys_cl, zeros(size(t)), t, x0); % 控制量 u = -Kx,取峰值 u_series = -K * x'; u_peak = max(abs(u_series)); fprintf('控制量峰值 = %f\n', u_peak);这里的逻辑是先用闭环系统做零输入响应,把状态轨迹 x 完整算出来,再通过 K 反推出控制量。初值选 [1; 0; 0] 表示给系统一个单位位移扰动,这是检验控制器最直接的方式。如果你发现 u_peak 已经超过执行器限幅,比如电机最大力矩是 10 而算出来是 30,那就说明极点拉太狠了,必须把 ωn 降下来,或者接受一个稍慢的响应。这段验证不是可选项,报告里放上控制量曲线,评审会认为你考虑过工程可实现性,这在评分里是明显的加分项。
4. 状态观测器与闭环仿真:从龙贝格增益到分离定理落地
4.1 龙贝格观测器:用对偶原理把 place 再抄一遍
实际工程里状态往往不能全部直接测量,比如质量弹簧阻尼系统里位移好测,速度要额外加传感器,成本立刻上去。观测器的意义就是用输出 y 和输入 u 把全部状态估计出来。龙贝格观测器的方程是:
x̂̇ = A·x̂ + B·u + L·(y - C·x̂)误差 e = x - x̂ 的动力学是 ė = (A - L·C)·e。你看这个形式和状态反馈的闭环矩阵 A - B·K 几乎一样,只是 K 换成了 L,B 换成了 C。所以设计方法直接复用对偶原理:把 A 转置,把 C 转置成“输入”,再用 place 求 L,最后再转置回来。
% 观测器极点:比控制器主导极点快 3~5 倍 p_obs = [3*real(p_dom(1)); 3*real(p_dom(1)); 3*real(p_dom(2))]; % 对偶设计:place(A', C', p_obs) 后转置 L = place(A', C', p_obs)'; % 验证观测器收敛动力学 eig_obs = eig(A - L*C); fprintf('观测器极点 = %s\n', mat2str(eig_obs)); % 增广闭环:真实状态在上,估计状态在下 A_aug = [A, -B*K; L*C, A - L*C - B*K]; x0_aug = [1; 0; 0; 0; 0; 0]; % 前三个是真实状态初值,后三个是估计初值 t = 0:0.01:10; % 零输入响应,增广系统 sys_aug = ss(A_aug, zeros(6, 1), zeros(1, 6), 0); y_aug = lsim(sys_aug, zeros(size(t)), t, x0_aug);这里最核心的逻辑是分离定理:闭环系统的极点等于状态反馈极点与观测器极点的并集。也就是说,你可以放心地先设计 K,再独立设计 L,两者不会互相干扰。但有个度要把握:观测器极点如果比控制器快太多,比如 10 倍以上,L 阵数值会变得很大,把测量噪声成倍放大到状态估计里,控制量曲线会毛刺频出。
p_obs 这个参数的设计习惯是取控制器主导极点实部的 3~5 倍。注意不要直接写 3p_dom,因为 p_dom 是复数,你乘出来还是复数,观测器用复极点本身没问题,但工程上为了简单,常用实极点。我这里取了两个相同的实部数值和一个稍远的实极点,确保观测器收敛速度高于闭环响应,同时 L 不会太大。增广矩阵 A_aug 的拼法也值得多看两眼:左上角是真实对象,右下角是观测器自身动力学,右上角的 -BK 表示控制器用的是估计状态,左下角的 L*C 是观测器从输出误差得到的修正。整个矩阵拼对了,仿真才能同时看到真实状态和估计状态的收敛过程。
4.2 拼闭环仿真:把真实状态与估计状态分开看
仿真时最容易犯的错误是让控制器直接用了真实状态 x,那观测器就成了摆设。正确做法是:被控对象模块里跑真实状态,控制器模块里跑估计状态 x̂,两者通过 u 连接,形成闭环。Simulink 里搭的话,被控对象是一个状态空间模块,观测器是一个自己搭的积分器加增益回路,控制器是一个增益矩阵 K。但脚本方式更快,也更适合放进报告。
% 从增广响应中拆出真实状态与估计状态 x_real = y_aug(:, 1:3); x_est = y_aug(:, 4:6); % 状态估计误差 err = x_real - x_est; % 画真实状态与估计状态:位移通道 plot(t, x_real(:, 1), 'b-', t, x_est(:, 1), 'r--'); legend('真实位移', '估计位移'); xlabel('时间 (s)'); ylabel('位移 (m)');这组曲线是报告里的核心证据。观测器初值设为零,真实状态初值有偏差,误差曲线从一开始的明显差值指数衰减到接近零,这个过程能直观说明观测器设计有效。参数上,初值偏差越大,误差消散过程越明显,但也别大到离谱,否则观测器的瞬态会被评分老师质疑。
如果你用的是 lsim,有个细节:增广系统是零输入系统,因为状态反馈 u = -K·x̂ 已经体现在 A_aug 的耦合项里了。所以第二个参数传 zeros(size(t)),意思是外部输入为零,整个动态全靠初值激发。这一步很多人想不通,以为还要把 u 单独加进去,结果矩阵维数对不上,报错后乱改一通。记住增广矩阵已经把 u 消掉了,你看到的响应完全来自初值 x0_aug。
5. 大作业避坑指南:评分点、数值病态与常见翻车现场
5.1 仿真一开始就发散成 NaN:先查极点,再查步长
现象:曲线前几秒直接飞上天,数值变成 NaN 或者 Inf。
原因:大多数时候是期望极点放太左,比如 -100,MATLAB 的默认变步长求解器在初始时刻无法分辨这么快的动态;另一种可能是 A-BK 算出来的特征值根本不是你要的极点,说明 K 设计有误,系统实际上处于正反馈状态。
解决:先把 eig(A - B*K) 打印出来和期望极点对比,排除设计错误。然后调小仿真步长或改用刚性求解器 ode15s。极点的实部绝对值决定了系统最快的动态时间常数,步长至少要小于这个时间常数的十分之一。如果这两个都查完仍然发散,再看看是不是观测器增益 L 和控制器 K 的符号接反了。血泪经验是,90% 的发散源于符号错误,而不是数值问题。
5.2 place 函数报错或返回全 NaN:先核对维数,再核对可控性
现象:place(A, B, p) 直接报错,或者返回的 K 全是 NaN。
原因:A、B 的维数不匹配;期望极点个数不等于系统阶数;系统本身有不可控模态;极点里有重根而矩阵又接近奇异。
解决:打印 size(A)、size(B)、length(p),确认三者关系是 size(A,1) = size(A,2) = length(p)。然后用 ctrb 和 rank 核对能控性,若秩不足,说明模型存在不可控模态,先做 minreal 或检查输入作用位置。如果只是重根问题,把两个相同的实极点改成相差 0.1 的两个实极点即可。注意 acker 在接近不可控的数值病态区域会给出极其夸张的 K,这时候宁可报错也比算出一个“看似能收敛”的假 K 好。
5.3 观测器把噪声放大到控制量里:极点数倍不是越大越好
现象:在输出端加一个小的测量噪声,控制量曲线立刻变成锯齿状,状态估计也有明显毛刺。
原因:观测器极点取太快,L 阵数值过大,测量噪声通过 L 直接注入状态估计回路。
解决:把观测器极点倍数从 5 降到 3,重新算 L。如果噪声仍然明显,在 Simulink 中对输出 y 加一个低通滤波器,但注意滤波器会引入相位滞后,可能让观测器收敛变慢,需要重新评估。另一种思路是改用卡尔曼滤波,它本身就是带噪声统计信息的最优观测器,只是大作业里通常作为加分项而非必须项。
5.4 能控性矩阵秩看着不够:数值病态与非最小实现
现象:手推明明能控的模型,rank(ctrb(A,B)) 就是比 n 小。
原因:参数量级差异过大,比如 m 是 1000、c 是 0.001,能控矩阵里元素跨越好几个数量级,奇异值分解判定时把极小值当成零;另一个常见原因是模型里有可约的零极点对,系统实际阶次比表象低。
解决:用 minreal 先消去可约模态,再取四矩阵;或者用 svd(Qc) 查看最小奇异值,若它只是比最大奇异值小 1e-12 这种级别,数值上可以认定能控,报告里附上奇异值说明即可。单位统一也是这步的隐形要求,角度用弧度别用度,长度用米别用厘米,否则病态问题会加重。
5.5 观测器方程符号写反:报告里最容易被抓的错误
现象:报告里写观测器方程时,把符号写成 x̂̇ = A·x̂ + B·u - L·(y - C·x̂),少了一个关键负号。
原因:直接从记忆里抄观测器形式,没有从误差动力学推导。既然误差动力学 ė = (A - L·C)·e,那么观测器必须带正号修正项,即 +L·(y - C·x̂),符号一错,整个稳定性分析全部站不住。
解决:写报告时先推导误差动力学,再反推观测器方程,顺序固定下来就不会错。这个细节在很多大作业验收现场被问到过,答不上来比写错更尴尬。状态反馈 u 的符号同理,u = -K·x̂,负号丢了就是正反馈,仿真直接跑飞。
6. 让大作业从“能跑”到“能答辩”:三个马上可用的收尾技巧
第一个技巧是把所有参数收敛到一个脚本头部,物理参数、期望极点、观测器倍数、仿真时长全部集中定义。这样做的好处是,验收时老师问“你把阻尼系数改成 0.1 会怎样”,你改一个数重跑一次就能当场回答。我见过太多人把参数散落在命令行历史里,验收时自己都说不清上次跑的是哪组参数,那种场面比答不出原理还难看。养成这个习惯后,报告里的每张图都能追溯到一个确定的脚本状态,这个可信度在评分里非常值钱。
第二个技巧是固定画三张图:状态响应、控制量曲线、真实与估计状态的误差曲线。这三张图正好对应设计目标、工程约束、观测器性能三个评分维度。画图时用统一的色号和线型,坐标轴标注清晰,把 K 阵和 L 阵的数值放在图注里。评分老师扫一眼就能看出你完整走通了“设计—验证—分析”链路,这比在报告里多写两页理论推导更有效。要注意控制量曲线必须画出峰值并注明是否超过执行器限幅,哪怕没有真执行器,这个讨论也能体现你的工程意识。
第三个技巧是把极点配置扩展成 LQR 对比。大作业做到极点配置这步,已经是完整闭环,但如果你想在答辩时多一个“为什么不用 LQR”的对话空间,可以顺手加一组对比:用 Q 对角阵和 R 调出另一组闭环极点,和手动极点配置结果放在同一张表里对比超调量与调节时间。LQR 给的是某种指标下的最优解,但极点位置不直观,这个对比本身就是很好的讨论素材。前提是你把 Q、R 取值也写进脚本头部,别让对比变成玄学。
我自己那门课拿到的教训是:能观性检查当时因为单位混用被判了秩不足,查了大半天才发现是长度单位不一致。从那以后,所有脚本第一行永远是参数声明和单位注释。希望这个习惯和上面这套流程能帮到你,让你的线性系统大作业一次跑通、顺利过审。
本文还有配套的精品资源,点击获取