简介围绕上海交通大学仪器科学与技术研究生课程“线性系统分析与设计”的大作业提供了一份完整的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 的峰值会跟着涨。执行器一旦饱和系统就从线性变成非线性前面算的极点全部失效。性能指标与极点参数的对应关系可以直接套这两个公式指标公式说明超调量 MPexp(-ζπ / √(1-ζ²))ζ 越小超调越大调节时间 ts4 / (ζ·ω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 换成了 LB 换成了 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。原因大多数时候是期望极点放太左比如 -100MATLAB 的默认变步长求解器在初始时刻无法分辨这么快的动态另一种可能是 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 取值也写进脚本头部别让对比变成玄学。我自己那门课拿到的教训是能观性检查当时因为单位混用被判了秩不足查了大半天才发现是长度单位不一致。从那以后所有脚本第一行永远是参数声明和单位注释。希望这个习惯和上面这套流程能帮到你让你的线性系统大作业一次跑通、顺利过审。本文还有配套的精品资源点击获取