简介上海交大仪器系研究生课程“线性系统分析与设计”的期末大作业以MATLAB为主要实现工具围绕连续/离散线性系统的建模、分析与校正展开。整个项目打包为rar格式共3个文件包括2个m脚本与1份docx报告压缩后仅193KB体量轻但内容完整便于快速检索。m代码侧重系统建模与仿真覆盖传递函数构建、Routh-Hurwitz与根轨迹稳定性分析、控制器参数整定以及阶跃/脉冲响应对比等典型环节docx报告则详细阐述了问题背景、设计流程和结果分析代码与文字相互印证便于复现和举一反三。已有680人学习下载对于正在修读线性系统理论、需要参考完整MATLAB课程设计或希望熟悉Control System Toolbox的读者是一份很实用的学习素材。1. 从公式到能跑的Matlab仿真这份线性系统大作业到底解决什么问题拿到某高校线性系统课程的大作业题目时我手里只有几个状态矩阵、一个“设计状态反馈使系统稳定”的需求和一张A4纸的评分标准。真正动手才发现公式推导是一回事把这些公式变成一段能在Matlab里跑通、能画图、能应付答辩演示的完整程序完全是另一回事。这份课程大作业资源就是这么一套东西把状态空间建模、能控能观性分析、状态反馈、观测器设计和离散仿真串成了一条完整链路适合正被作业卡住的学生和想快速上手状态空间方法的工程师。下文按复现顺序讲每段代码可独立执行。2. 建模与状态空间化写Matlab之前先想清楚这三件事状态反馈听起来高大上但第一步是把被控对象变成A、B、C、D四个矩阵。这一章先把对象模型选型讲清楚再落到Matlab代码最后说采样周期的坑——这三件事没理顺后面算出来的K再漂亮也是空中楼阁。2.1 为什么是带执行机构动态的三阶模型而不是教材双积分器很多课本例题喜欢用二阶双积分器因为数学上清爽能控能观矩阵看一眼就能求画个根轨迹就收工。但真实的大作业场景里这种模型太“娇贵”它只有两个积分环节状态反馈一加上去闭环特征方程简单到让人看不出设计空间。实际上课程大作业为了考察完整的现代控制流程通常会引入执行机构动态把对象变成三阶也就是双积分器串联一个一阶惯性环节。我用的对象模型是G(s) 1 / [s²(τs 1)]其中时间常数τ1s代表一个响应偏慢的执行机构。选这个结构出于三点考虑第一它包含积分环节能模拟位置伺服系统这类真实对象后边算状态反馈时能直接对应到“控制器输出到执行器再到机械位置”的物理链条第二三阶系统开环不稳定必须靠反馈拉稳能明显体现状态反馈的作用第三三个状态对应位置、速度、执行机构输出答辩时每个状态都能讲出物理意义不至于被追问时露怯。对应的状态空间表示如下状态变量定义为 x₁位置x₂速度x₃执行机构输出的加速度指令ẋ Ax Bu y Cx Du矩阵取值矩阵取值含义A[0 1 0; 0 0 1; 0 0 -1]位置→速度、速度→加速度、执行机构一阶滞后B[0; 0; 1]控制电压直接作用在执行机构输入端C[1 0 0]测量输出为位置D0无直接馈通这里x₃是执行机构输出不是加速度本身也能讲但把它当成“经过一阶滞后后的控制作用”最直观u经过一个时间常数为1s的惯性环节后变成x₃再积分两次得到位置。这种状态选取不是唯一答案但它是与物理过程对应最直接的一组后边设计观测器时也容易设置初值。2.2 把题干写进Matlab状态空间对象与交叉验证建模这一步我习惯直接写矩阵而不是用传递函数转换因为矩阵写出来能自己核对每一行物理意义。最小工作代码% 定义三阶对象模型 A [0 1 0; 0 0 1; 0 0 -1]; % 系统矩阵最后一行为执行机构动态 B [0; 0; 1]; % 输入矩阵控制量作用在第3个状态 C [1 0 0]; % 输出矩阵测量位置x1 D 0; % 直接馈通为零 % 生成连续时间状态空间对象 sys ss(A, B, C, D); % 检查维度一致性并显示对象信息 n size(A, 1); % 状态数应为3 m size(B, 2); % 输入数应为1 p size(C, 1); % 输出数应为1 disp([系统阶数 n , num2str(n)]);这里的size检查很多人会跳过但我建议保留后续代码里如果出现B矩阵维度写错比如忘了给第三行写1ss会因为维度不匹配直接报错到时候排错比这多花十分钟。sys对象建立后后续阶跃响应、bode图、c2d离散化都可以直接复用它不需要重复写矩阵。除了直接写矩阵我还会做一个交叉验证把同一对象写成传递函数再用tf2ss转换对比特征根是否一致。这一步不能证明状态选取唯一但能确认建模没有笔误。% 用传递函数交叉验证 G tf(1, [1 1 0 0]); % 分母: s^3 s^2对应 1/(s^2(s1)) sys_tf ss(G); % 转换为状态空间 disp(eig(sys.A)); % 应该为 0, 0, -1 disp(eig(sys_tf.A)); % 与上行一致即可说明一下tf的分子为1分母写成 [1 1 0 0]对应 s³s²展开就是 s²(s1)。eig比较的是状态矩阵A的特征值正常情况下都是0、0、-1。如果两次结果不一致说明矩阵符号或行序有误。这种“从传递函数换算一遍”的交叉验证是建模阶段最便宜的后悔药。2.3 连续域还是离散域采样周期的约定状态反馈算出来之后很多人习惯直接在连续域里跑Simulink然后把控制器原样搬进数字系统结果走了个大弯路。我的做法是一开始就定好两个仿真层级连续域用于验证控制算法的数学正确性离散域用于验证工程可实现性。两者差一个环节零阶保持器。离散化这一步的核心参数是采样周期Ts。对三阶对象我一般先按系统带宽的10到20倍来选。以本对象为例开环有一个零特征根闭环期望带宽大约在25rad/s那么Ts取0.01s到0.05s都是合理的。太大会引入明显相位滞后导致离散控制器与连续设计脱节太小则仿真步数多运行慢且没有实际意义。% 采样周期选择与离散化 Ts 0.01; % 采样周期10ms工程中常见 sys_d c2d(sys, Ts, zoh); % 零阶保持离散化 % 观察离散能控性矩阵是否保持满秩 Co_d ctrb(sys_d.A, sys_d.B); disp([离散系统能控性秩: , num2str(rank(Co_d))]);关于zoh和foh的选择控制系统仿真里默认用zoh因为它模拟的是真实数字控制器的信号保持行为。foh在某些信号重建场景更平滑但大作业没必要引入这种复杂度。离散化后必须重新检查能控性这一点与连续系统不是一回事具体原因在第3章末尾细说。3. 能控能观性检验rank命令算出1之前先理解矩阵在测什么状态反馈设计的前提是系统能控观测器设计的前提是系统能观。这一步教材上有现成判据但Matlab里rank函数的结果经常让人摸不着头脑明明数学上满秩显示出来却是2。这章先把矩阵判据讲透再给出带容差处理的实操代码。3.1 能控矩阵的列向量到底在测什么能控性矩阵定义是 Co [B AB A²B ...]对三阶系统就是 [B AB A²B]。每一列代表的含义B列是控制量直接作用的状态方向AB列是控制量经过一步系统动态后影响的状态方向A²B列是经过两步后影响的方向。如果这个矩阵秩等于状态数n说明控制量通过最多n步的动态传播能覆盖整个状态空间。反过来秩不足就意味着存在某个状态方向永远不在控制量的可达范围内。物理上常表现为某个状态变量与控制通道之间存在结构性的解耦是零极点对消造成的。比如对象传递函数分子上有一个与系统极点相同的零点那个模态就被输入“隔离”了。能观性矩阵 Ob [C; CA; CA²]测的是另一个方向从输出y的当前值和各阶导数值中能否反推出全部状态初始值。对三阶系统如果Ob满秩就能从位置测量值及其一阶、二阶导数信息中恢复速度与执行机构状态。这在后边设计观测器时是必要前提。3.2 Matlab实操判定代码与阈值容差直接写rank(ctrb(sys.A, sys.B))是最常见的做法但我不推荐只用这一句交差。原因是Matlab的rank默认带了一个与机器精度相关的容差默认约等于 max(size(A)) * eps * norm(A)。若系统矩阵数值刚好让某个奇异值落在容差附近rank会返回一个偏小的值造成误判。更稳妥的写法是同时计算秩和最小奇异值并做一个显式的阈值判断% 能控能观性完整检验 Co ctrb(A, B); % 能控性矩阵 3x3 Ob obsv(A, C); % 能观性矩阵 3x3 % 分别计算奇异值 s_Co svd(Co); s_Ob svd(Ob); % 相对阈值小于最大奇异值的1e-10则视为数值零 tol 1e-10; r_Co sum(s_Co max(s_Co) * tol); r_Ob sum(s_Ob max(s_Ob) * tol); disp([能控性秩: , num2str(r_Co), / , num2str(size(A,1))]); disp([能观性秩: , num2str(r_Ob), / , num2str(size(A,1))]);这里tol用的是相对阈值比rank默认容差更直观理解如果一个奇异值比最大的小10个数量级以上说明对应方向上的能控/能观性弱到工程上可以忽略。对本文三阶模型Co和Ob的理论秩均为3上面代码输出应为3/3和3/3。如果你计算出的最小奇异值在1e-3量级而最大奇异值是1e2比值1e-5系统虽然数学满秩但工程意义上已经接近不可控控制器增益会很大——这是后边第5章要展开的坑。3.3 Gram矩阵与PBH判据答辩备用验证当rank结果模棱两可时另一个更平滑的判断工具是能控性格拉姆矩阵。连续系统能控性等价于格拉姆矩阵正定它能给出“控制能量”的量化视角最小特征值小说明让某个状态方向收敛需要很大的控制能量。% 能控性格拉姆矩阵 Wc gram(sys, c); disp(eig(Wc)); % 三个特征值均应显著大于0 % PBH判据: 对每个特征值 lambda检查 [lambda*I-A B] 行满秩 lambda eig(A); for i 1:length(lambda) M [lambda(i)*eye(n) - A, B]; if rank(M) n disp([特征值 , num2str(lambda(i)), 不可控]); end endPBH判据的实用价值在于定位不可控模态哪条特征值对应的矩阵缺秩就是哪个模态失控。应付答辩时老师问“你怎么证明极点都能任意配置”用PBH判据作答比只贴rank结果更有说服力因为结论直接与特征值绑定。末尾提醒一个容易忽略的点系统在连续域能控离散化后不一定能控。当采样周期Ts恰好使某个连续模态变成离散域中的采样零点离散能控矩阵秩就会下降。这就是上一章结尾我坚持要重算ctrb(sys_d.A, sys_d.B)的原因。对本文Ts0.01s不存在这个问题但若你把Ts换成1s会发现原本可控的系统离散后接近失控。4. 状态反馈与观测器极点怎么选、增益怎么算、分离定理怎么用能控能观满秩后真正的大作业重头戏才开始极点配置。大多数人卡在不是不会调acker而是不知道怎么选极点以及为什么选了观测器极点后整个系统就稳定了。这章把选择依据、代码实现和分离定理串在一起讲。4.1 极点位置与系统响应的权衡状态反馈的本质是把闭环极点移动到设计者指定的位置。极点实部绝对值越大响应越快但控制量峰值也越大虚部不为零时产生振铃阻尼比越小超调越大。我整理三组参考极点极点组合阻尼比/形式预期效果控制量特点-2±2j, -50.707超调约5%响应中等K元素适中执行器压力小-3±1.5j, -60.894超调3%以内调节时间更短K元素较大控制量峰值翻倍-1±3j, -40.316明显振荡第一个超调可能上20%反馈增益偏小但系统动态差初学者最容易犯的错是把极点全部设在实轴上并且离原点很远比如[-10, -11, -12]。这样闭环响应确实“快”但K矩阵各分量常常达到几十上百控制量远超执行器限幅仿真图上响应快得离谱实际系统一台机都动不了。我一般先按阻尼比0.707、自然频率23rad/s来定主导极点第三个非主导极点取-5保证距离主导极点约2倍以上避免第三极点影响主导模式。4.2 用acker还是place计算状态反馈增益的姿势Matlab提供两个核心函数acker单输入系统和place通用多输入。对单输入大作业acker够用且代码短但它在极点非常接近时数值上容易崩。place在MIMO和重极点场景更稳但要求期望极点不能重复得太严苛。我的习惯是先acker出现数值告警立刻切place。% 期望极点主导极点阻尼0.707自然频率3rad/s p_dom -3*0.707 1i*3*sqrt(1-0.707^2); % -2.12 2.12j p3 -5; % 非主导极点 poles [p_dom; conj(p_dom); p3]; % 状态反馈增益计算 K acker(A, B, poles); % 单输入系统直接算行向量 % K place(A, B, poles); % 若出现重极点警告时换成这个 % 闭环稳定性验证 Acl A - B*K; disp(闭环特征值:); disp(eig(Acl)); % 检查K元素幅度是否合理 disp([反馈增益 max(abs(K)) , num2str(max(abs(K)))]);代码中conj(p_dom)是补齐共轭极点这点不能漏否则特征多项式系数是复数物理不可实现。K的计算结果是一个1×3的行向量对应三个状态各自的反馈通道。注意eig(Acl)必须与poles一致差一位说明符号写错A-BK和ABK只差一个负号输出完全不同。K元素幅度是常被忽略的量级检查。对本文模型合理期望是max(abs(K))不超过20。如果你算出来K里有三位数多半是极点放得过远或系统病态回头重新权衡响应速度。4.3 观测器设计与分离定理的工程含义状态反馈设计假定全部状态可测但大作业通常只会告诉你输出是位置速度和执行机构状态要靠状态观测器重建。观测器增益L的求法本质上与K对称只不过把极点配置用在误差动态上% 观测器极点比反馈极点快3倍左右 obs_poles 3 * real(poles); % 实部放大3倍 obs_poles obs_poles - 1i * imag(poles)*0.5; % 适当减小虚部防噪声放大 % 求观测器增益列向量 L place(A, C, obs_poles); % 转置是经典坑少一个转置直接算错 % 验证观测器误差动态 A_obs A - L*C; disp(观测器误差特征值:); disp(eig(A_obs));place用于观测器时第一个参数是A第二个是C算完结果再转置回来得到L。原因是误差动态方程的特征值与A-LC的转置A-CL相同写成转置形式才能复用place接口。obs_poles取反馈极点的3倍经验上能保证观察误差收敛速度明显快于控制器响应但又不会把测量噪声放大到不可接受。分离定理说组合系统的闭环特征值等于反馈极点与观测误差极点的并集设计上可分别独立进行。为了向自己证明这一点可以把两个子系统拼成四阶组合矩阵验证% 组合系统矩阵验证分离定理 A_full [A - B*K, B*K; zeros(n,n), A - L*C]; eig_full eig(A_full); % 分离后对比 eig_sep [eig(A - B*K); eig(A - L*C)]; disp(组合特征值误差:); disp(max(abs(sort(eig_full) - sort(eig_sep))));数值上两者必须逐位相同误差通常在1e-12左右。这个大括号里的验证代码答辩时直接展示比空口解释分离定理有说服力得多。4.4 备选方案时间紧直接上LQR如果不想花时间挑极点LQR是一条兜底路线。LQR本质上是通过最小化二次型代价函数来确定K不需要指定极点位置但需要设定Q和R矩阵。我建议Q取状态权重R取控制权重% LQR 设计备选 Q_lqr C*C 0.1*eye(n); % 输出权重轻微状态惩罚 R_lqr 1; % 控制能量权重 [K_lqr, ~, ~] lqr(A, B, Q_lqr, R_lqr); disp(LQR反馈增益:); disp(K_lqr);Q矩阵中C*C是输出项yy的扩展意思是重点惩罚输出偏差加0.1eye(n)是让三个状态均有一定惩罚防止某个状态被“无视”。R越小控制越激进R越大增益越小响应越慢。LQR的好处是不用担心极点选错缺点是答辩时难以讲清“为什么Q取0.1而不是0.2”所以做演示我还是推荐用ackerplace的主线LQR留作补充对比。5. 大作业最常翻车的五个坑仿真发散不等于你算错了把状态反馈和观测器放进仿真之前先把自己从“下一步做什么”切换到“哪里可能出错”。这五个坑是我反复在类似任务里遇到的每一条都真实发生过且每条都有明确的排查路线。5.1 现象rank(ctrb(A,B)) 显示2但理论上明明满秩查代码半天发现矩阵没写错能控性矩阵算出来也是3×3可是rank返回2。原因在于Matlab的rank默认容差是相对机器精度的当A矩阵较大或状态变量间数值相差悬殊时最小奇异值落在容差附近会被当成零。解决方法是显式指定容差或像第3章那样用svd看完整奇异值分布rank(ctrb(A,B), 1e-8) % 显式指定容差 % 或 s svd(ctrb(A,B)); % 看全部奇异值 disp(min(s) / max(s)); % 比值小于1e-12再考虑判定不可控如果比值在1e-4到1e-8之间虽然rank满秩也要小心对应的状态方向能控性很弱反馈增益会异常大。这个坑的教训是矩阵判据在纸上是个0/1问题在Matlab里是个数值精度问题不能只看rank返回值。5.2 现象acker 报错“极点太近”place 返回的 K 大得离谱直接套用期望极点集合[-5, -5.0001, -6]时acker会因极点接近无法求解place能用但返回的K分分钟上百。原因是极点接近会使得特征多项式在所求区域ill-conditioned闭环系统对参数扰动极度敏感。解决思路是把极点距离拉开要么均匀分布在一条圆弧上要么保留至少0.51的实部间隔。% 拉开间距的极点设置 r 3; % 半径 p r * exp(1i*(pi (0:2)*2*pi/3)); % 在半径3的圆上均布三个极点 K place(A, B, p); disp(max(abs(K)));均布极点的好处是任意两个极点之间角度差120°最小间距大数值求解稳定。算完务必看一眼max(abs(K))超过30就该减小半径重新来。别为了追求“理想超调”牺牲控制量。5.3 现象Simulink 一运行就刷 NaN频率高到屏幕闪仿真发散最常见原因不是公式错而是初始条件不一致观测器初值取0但真实系统初值在[1;0;0]误差动态虽稳定但初值误差大触发非线性环节数值溢出。另一个高频原因是A-B*K本身不稳定极点没配上或状态反馈加进了Simulink时把K矩阵的乘号顺序写反。排查路径按顺序做先在命令行跑step(ss(A-B*K, B, C, D))如果这里也发散是极点配置的问题先回头查K如果这里稳定而Simulink发散再去Simulink里核对State Space模块的A、B、C、D参数和初值。我的习惯是Simulink里的状态空间初值与真实系统设成一致观测器初值单独设到另一个Constant模块避免手滑。5.4 现象离散化后控制电压高频抖动连续域仿真却一切正常连续域仿真里K算好、响应平滑换上c2d离散对象后控制量u-K*x在每个采样点高频振荡。原因是采样周期Ts太大离散控制器的相位滞后让闭环在高频段失去稳定裕度另一个常见原因是观测器极点选得太快比如取反馈极点的10倍以上导致离散化后观测器放大高频噪声经K反馈后出现抖振。解决方向把Ts从0.05s降到0.01s或0.005s观测器极点与反馈极点的倍数控制在35倍以内并给测量输出加一个一阶低通滤波器。离散化后一定要检查离散闭环极点是否都在单位圆内并且离单位圆周有足够余量eig_d eig(sys_d.A - sys_d.B*K); disp(max(abs(eig_d))); % 必须小于1且小于0.95更稳妥如果最大模在0.995说明离散闭环已经很接近失稳边界无论连续域设计多完美上数字系统都会抖。这一步是离散域仿真的核心验收线。5.5 现象仿真曲线很漂亮答辩被问“控制量峰值是多少”直接卡住这条不属于软件报错而是设计与工程脱节。大作业评分常见追问就是“执行器限幅/-10V你的控制量达到多少”。很多人的K设计只盯响应曲线u-Kx的结果一眼没看结果u峰值超过100一看就是不可实现的。解决方法是固定一个饱和模块把输出限幅接到被控对象前再做仿真然后观察闭环是否还能稳住。% 控制峰值快速计算 x0 [1; 0; 0]; t 0:0.001:5; u_max 0; for i 1:length(t) ut -K * (expm(Acl*t(i)) * x0); % 假设无观测器动态 u_max max(u_max, abs(ut)); end disp([控制峰值约: , num2str(u_max)]);这个循环本质是在无噪声、初值1的情况下扫一遍控制量轨迹。若超限第一步想的不应该是增大限幅而是收回极点半径或调小Q中状态权重。反馈设计除了数学稳定还要在物理约束内稳定否则答辩一问就翻车。下面这个速查表可以打印出来贴屏幕边上现象首要怀疑优先排查方式rank显示缺秩数值容差svd看奇异值比值acker崩溃极点过近换place或拉开极点间距Simulink发散初值或K符号先跑step(A-B*K)隔离问题离散后抖动Ts或观测器极点查看离散极点最大模控制量超限极点过远计算u峰值并减极点半径6. 验证与扩展从跑通到能答辩的一次参数扫描跑通仿真只是拿到了及格线。大作业答辩时老师真正想听的是你的设计在参数变化下是否依然成立为什么这组极点合适。这就需要一个可复现的验证脚本我的习惯是做一次带增益扫描和极点扫描的批处理最后生成一张包络图。增益裕度检验是第一关把反馈增益乘以0.5倍和1.5倍观察闭环特征值和响应曲线是否仍稳定。这个操作等价于模拟系统模型误差15%甚至50%时控制器的强壮性。% 增益裕度扫描: 验证鲁棒性 alphas 0.5:0.1:1.5; % 增益缩放系数 max_real zeros(size(alphas)); for i 1:length(alphas) K_test alphas(i) * K; e eig(A - B*K_test); max_real(i) max(real(e)); end plot(alphas, max_real); grid on; xlabel(增益缩放系数); ylabel(闭环最大实部); title(闭环极点随增益变化);当缩放系数从0.5升到1.5时max_real始终小于0说明系统至少具备±50%的增益裕度。若某一段越过0线说明设计刚好压着稳定边界答辩时可能会被追问“模型不准确怎么办”。这个脚本是整个设计收口前必跑一遍的标准动作。参数扫描进阶玩法是对阻尼比和自然频率做二维网格搜索。对每组参数算K、跑一次sim将超调量、调节时间、控制峰值记录成表格最终选出满足所有约束的最优折中。% 极点参数扫描阻尼比 zeta 与自然频率 omega zeta_range [0.5 0.707 0.9]; omega_range [2 4 8]; result []; for zeta zeta_range for omega omega_range p1 -zeta*omega 1i*omega*sqrt(1-zeta^2); p [p1; conj(p1); -5]; K_t acker(A, B, p); % 直接解析计算第2个状态的零输入响应过冲 Acl_t A - B*K_t; % 用step获取单位阶跃响应信息 y step(ss(Acl_t, B, C, D), 0:0.01:5); overshoot (max(y) - y(end)) / y(end) * 100; settle_idx find(abs(y - y(end)) 0.02*y(end), 1, last); settle_time 0.01 * (length(y) - settle_idx 1); result [result; zeta, omega, overshoot, settle_time, max(abs(K_t))]; end end % 按超调量排序打印候选设计 disp(array2table(result, ... VariableNames, {zeta, omega, overshoot, settle, Kmax}));说明一下循环逻辑外循环换阻尼比内循环换自然频率每组重新用acker解一次增益再基于闭环状态空间算阶跃响应指标。结果矩阵每一行是一组候选设计最后一列是K的最大绝对值用于评估控制量负担。我通常把超调小于10%、调节时间小于3s、Kmax小于30的行勾出来那些行就是答辩PPT上的主推方案。有一次我赶进度把增益扫描跑了但没仔细看结果闭环保有极点刚好压在-0.05虚线一画看似稳定结果模型参数改动10%后直接翻到右半平面。从那以后每次交这类作业我强制自己把“增益扫描极点扫描控制量峰值”三张图全部生成并检查一遍确保所有曲线在最恶劣参数下仍然稳定。这套操作也能顺带回答答辩老师最常问的一个问题“你没有没做过鲁棒性分析”希望帮到你。本文还有配套的精品资源点击获取