简介上海交大仪器科学与技术研究生课程《线性系统分析与设计》的大作业资源面向需要完成线性系统建模、稳定性分析与控制器设计任务的研究生适合配合MATLAB控制系统工具箱和Simulink仿真学习使用。压缩包为rar格式共3个文件包含2个MATLAB脚本和1份Word文档脚本负责系统分析与设计计算文档整理完整设计过程与结论整体仅193KB内容精炼。目前已有680人学习可作为课程作业的完整参考。作业覆盖系统建模、稳定性判据应用、频率响应分析、PID校正以及阶跃/脉冲响应对比等关键环节读者可借此理解线性系统理论从推导到代码实现的完整路径并借鉴其报告结构、仿真设置与分析思路。1. 某高校线性系统大作业一个把现代控制理论从课本拽进代码的完整闭环很多人在上完线性系统课程之后依然说不清状态反馈、观测器、能控性这些概念到底能拿来干什么——直到他需要交一份大作业要求把某个真实物理对象建模成状态空间表达式再设计控制器让它达到指定的动态指标。这篇笔记要讲的就是某高校线性系统大作业这类任务最常见的完整拆解路径从选角度、建模型、判能控能观到极点配置、设计状态观测器再到仿真验证和参数整定。它面向的不是想混过验收的人而是想知道“用什么思路做出来、每一步卡点在哪、怎么调参不翻车”的从业者和高年级学生。这套流程其实和工程里做一款运动控制系统的原型验证是同一条路唯一区别是作业里没有成本上限和量产压力。读完你会得到一份能直接照着改的模板工程以及几条花钱买不来的踩坑记录。2. 从题目到方案选哪个被控对象决定了你后面会不会舒服2.1 为什么说选对象是整份作业的地基大作业通常不限定具体被控对象只要求“一个物理系统 完整的现代控制设计流程”。最常见的翻车点是选了传递函数只有二阶但带严重非线性、或者模型参数根本查不到的系统导致后续建模和仿真处处受制。我一般建议选一个“结构简单、能写出解释清楚的微分方程、参数可在合理范围内自定”的对象比如直流电机位置伺服系统、弹簧-质量-阻尼系统或者单级倒立摆的线性化模型。直流电机最合适它的电气时间常数远小于机械时间常数可以自然简化成二阶模型方便验证能控性而且状态变量——电枢电流和转速——都有明确的物理意义画状态反馈框图时也能给评委讲清楚每个增益的含义。2.2 数学建模的常见做法怎么从物理方程转成标准状态空间以他励直流电机带动惯性负载为例忽略电感动态即认为电枢回路电感 L 很小可近似为 0机械方程为 ( J\frac{d\omega}{dt} K_t i - b\omega )电气方程为 ( u R i K_e \omega )。取状态 ( x_1 \omega )转速、( x_2 i )电枢电流输入为电枢电压 ( u )输出为转速 ( x_1 )可以得到下面这个标准形式。import numpy as np # 电机参数转动惯量 J粘性摩擦系数 b转矩常数 Kt反电动势常数 Ke电枢电阻 R J, b, Kt, Ke, R 0.01, 0.0002, 0.05, 0.05, 1.0 A np.array([[-b/J, Kt/J], [-Ke/R, -1.0/R]]) # 注意这里的第二行是推导后的简化形式 B np.array([[0.0], [1.0/R]]) C np.array([[1.0, 0.0]]) D np.array([[0.0]]) print(A , A) print(B , B)这段代码有几处需要特别说明。A矩阵第二行的写法背后有推导细节把电气方程离散掉电感后电流的微分方程实际上约束退化为代数关系这里为了保持状态空间形式一种常见处理是引入足够小的电感值后取极限或者直接使用包含电感的完整四阶模型再降阶。如果评审老师不接受简化模型就改成两个状态分别为电机的电枢电流和转速、并保留电感 L 的完整模型代码里的A[1,0]和A[1,1]会多出-Ke/L与-R/L项形态对应的是“电流的导数等于电压减去反电动势减去电阻压降再除以电感”。另一个大坑是单位制如果转速用 rpm、转矩用 kg·cmKt 和 Ke 的换算差一个系数后面算出的极点位置会完全偏离预期。我习惯全部转成国际单位rad/s、N·m、V、A。写完A,B,C,D后执行一次np.linalg.matrix_rank(ctrb(A,B))确认能控性矩阵满秩这是整个后续设计能成立的前提。检查能观性同样简单。2.3 能控性和能观性为什么这一步会在答辩时被追问能控性回答的问题是“系统状态能不能在有限时间内被输入驱动到任意目标”能观性回答“能不能通过输出测量把全部状态估计出来”。对直流电机模型能控性矩阵满秩几乎是必然的——毕竟电压能直接影响电流电流又作用于转速。但能观性这里就有一个典型陷阱如果输出只测转速而电流状态对转速的作用被模型近似削弱到几乎不可观那么设计观测器时会发现极点的实数部分虚得厉害表现为估计值迟迟不收敛。遇到这种情况要么改输出为转速加电流联测要么在观测器设计里接受较慢的观测动态。答辩时老师通常会追问“你的输出矩阵为什么不选单位阵”这时如果你能直接说出“物理可测性优先于数学完备性”并且补充一句“我在观测器里用到了输出注入来补偿不可测状态的影响”这一问就算过关了。3. 状态反馈与极点配置把系统的动态特性按指标“焊死”3.1 性能指标怎么转成期望极点大作业要求通常写“超调量小于 10%调节时间小于 0.5 s”这两个指标需要转成二阶系统的期望极点位置。假设主导极点为一对共轭复根阻尼比 ( \zeta ) 由超调量公式 ( \sigma e^{-\pi\zeta/\sqrt{1-\zeta^2}} ) 反推调节时间按 ( t_s 4 / (\zeta\omega_n) ) 估算。对 10% 超调量( \zeta \approx 0.59 )对 0.5 s 调节时间并留一点余量取 ( \zeta\omega_n \approx 10 )那么 ( \omega_n \approx 17 )对应极点为 ( -10 \pm j 13.8 )。这个计算过程要用代码固化下来因为评审可能会要求你解释每个数值的来历而不是拍脑袋给一个极点。import numpy as np # 根据超调量求阻尼比 sigma_desired 0.10 zeta -np.log(sigma_desired) / np.sqrt(np.pi**2 np.log(sigma_desired)**2) print(阻尼比 zeta , round(zeta, 3)) # 根据调节时间求自然频率 ts_desired 0.5 wn 4.0 / (zeta * ts_desired) print(自然频率 wn , round(wn, 3)) # 期望闭环极点一对主导极点 real_part -zeta * wn imag_part wn * np.sqrt(1 - zeta**2) poles_desired [complex(real_part, imag_part), complex(real_part, -imag_part)] print(期望极点 , poles_desired)上面这段代码的作用是核心换算逻辑上先由时域指标反推阻尼比和自然频率再构造闭环极点候选。参数说明因为调节时间与极点的实部绝对值直接相关所以算出来的real_part必须严格落在设计带宽之内如果调试时发现稳态误差太大可以增加一个积分环节而不是把极点越拉越左否则控制量会饱和。这个取舍在答辩时也是加分项说明“我把主导极点放在预期位置同时保留系统零点的影响评估而不是一味提高增益”。算出期望极点后下一步就是用place函数反解状态反馈增益。3.2 用python-control的place求状态反馈增益命令、参数与多解风险place是python-control里做极点配置的标准入口底层使用鲁棒特征结构配置算法。先对前面的电机模型封装成状态空间对象再调用它import control as ct # 从上一节的 A, B, C, D 构造连续时间状态空间系统 sys_plant ct.ss(A, B, C, D) print(系统能控性矩阵的秩:, np.linalg.matrix_rank(ct.ctrb(A, B))) # 期望极点需要是共轭成对出现的复数place 会自动处理 K ct.place(A, B, poles_desired) print(状态反馈增益 K , K) A_cl A - B K print(闭环特征值:, np.linalg.eigvals(A_cl))这里要重点说清楚place的参数行为第一个参数是系统矩阵第二个是输入矩阵第三个是期望极点列表其中复数极点必须成对给出否则函数直接报错。另一个容易踩的坑是期望极点数必须等于系统状态维数直流电机的二阶模型就需要两个极点如果你在 3.1 里只给了一对主导极点还不打算保留其余极点可能要从系统零极点里找到第三个极点并安排到更左侧但要注意place会同时移动所有极点不是只处理你指定的那对。代码最后打印的闭环特征值用来校验极点配置确实生效这是最容易自查的一步如果算出来和期望极点不一致优先检查A - B K这一步是不是用了np.dot导致矩阵乘法方向反了。多解问题上place会在满足极点条件的解里选择有较好条件数的那一个所以不需要你自己做线性方程组的拼装。3.3 状态反馈的物理限制控制量饱和不是 MATLAB 能帮你解决的问题理论设计完成后仿真里最容易出现的问题就是控制电压超出执行器上限。你在极点配置里把极点放得越远需要的控制能量越大而在参数表里写上“电枢电压最大 24 V”之后控制器输出却经常跑出 200 V。这背后的原因是状态反馈增益是根据模型求出的但模型的线性范围假设了一个不受限的输入而真实的物理系统有输入饱和。仿真里加入限幅块或者代码里做np.clip(u, -24, 24)后你通常会看到系统变得松弛甚至出现极限环。解决的常见做法有两个方向其一设计极点时留足带宽余量不要让wn超过物理可实现范围其二在控制器输出端加 anti-windup 结构。作业层面我更推荐第一种因为它更简单而且能在答辩里讲清楚“我用极点位置反映了执行器能力”不需要引入额外的积分器。4. 状态不可测怎么办设计龙伯格观测器并做仿真联调4.1 全维观测器还是降维观测器按可测状态数量和噪声情况选状态反馈的前提是全状态可测但直流电机里电枢电流往往是通过采样电阻估算出来的噪声大成本高。如果你的输出只能测转速那两个状态里有一个要依赖观测器来重构。全维观测器的含义是用系统模型加输出误差修正项来估计全部状态实现简单降维观测器只估计不可测的那部分收敛性能更好但推导复杂、容易出符号错误。对大作业来说除非老师明确要求设计降维观测器否则全维观测器是性价比最高的选择。它在工程上的对应物是各种基于模型的状态估计方案区别只是工作点在稳态还是动态过程。要注意的是观测器的极点一般取期望闭环反馈极点的 2 到 5 倍远但也不能无限远因为观测器带宽过高会把测量噪声放大进状态估计进而污染控制量这在仿真的噪声测试里会直接看到“锯齿状”的转速曲线。4.2 用observer函数完成观测器增益设计的完整代码python-control没有单独的全维观测器函数但可以借用极点配置方法因为观测器增益 L 的本质是把对偶系统的极点放到目标位置。写代码时先把系统转置再调用place最后转置回来。import control as ct import numpy as np # 设定观测器期望极点是反馈极点实部的 3 倍左右且共轭成对 obs_poles [complex(-30, 30), complex(-30, -30)] # 全维观测器增益 L 的计算对偶系统极点配置 A_T np.transpose(A) C_T np.transpose(C) L np.transpose(ct.place(A_T, C_T, obs_poles)) print(观测器增益 L , L) # 构造观测器系统矩阵 A_obs A - L C print(观测器特征值:, np.linalg.eigvals(A_obs))这段代码的注意点在于place的入参形态对偶系统里“输入矩阵”是 ( C^T )输出矩阵是 ( B^T )期望极点的个数依然是系统阶数。在代码逻辑上ct.place返回的增益矩形成型是(1, n)转置之后L正好是(n, 1)的列向量才能正确乘在C后面。有同学在这里翻车把转置漏掉然后A - L C的维度直接不匹配合报错信息是“形状不一致”会误导人往矩阵乘法方向排查。如果极点是复数place会在内部自动处理共轭但你在传参时仍然需要显式地构造成对复数列表。打印的A_obs特征值是用来自检的如果它们和obs_poles对不上多半是C取错了行而不是算法有问题。4.3 状态反馈和观测器的组合仿真闭环矩阵的装配套路控制器实际使用的是估计状态所以闭环系统写为 ((A - BK)) 在真实状态上作用而观测器误差动态 (e) 是独立的。仿真时有两种接法一种是把真实模型和观测器拼成一个更大的状态向量用scipy.integrate.solve_ivp联立积分另一种是把输出直接接回控制律形成带反馈的闭环系统再用ct.forced_response仿真。第一种更直观也更适合检查“估计状态是否收敛到真实状态”。import numpy as np from scipy.integrate import solve_ivp # 闭环系统组装真实状态 x估计状态 x_hat def closed_loop_dynamics(t, x_all): x x_all[:2].reshape(2, 1) x_hat x_all[2:].reshape(2, 1) u -K x_hat # 控制器用的是估计状态 dx A x B u dx_hat A x_hat B u L (C x - C x_hat) return np.concatenate([dx.ravel(), dx_hat.ravel()]) x0 np.array([0.0, 0.0, 0.2, 0.0]) # 真实转速初始为0观测器初始转速0.2 rad/s T np.linspace(0, 1.0, 1000) sol solve_ivp(closed_loop_dynamics, [0, 1.0], x0, t_evalT, methodRK45)这里的参数有讲究x0里故意让观测器的转速初值和真实值不一致目的就是看收敛过程。如果观测器开始就设成与真实值一致仿真曲线会完美重合就验证不了观测器的纠错能力。u的表达式里的负号来自反馈方向的定义写成u -K x_hat时状态反馈回路是负反馈对应闭环矩阵A - BK一旦方向搞反系统立刻发散仿真结果会出现指数增长曲线而不是收敛。如果仿真步长太大会看到估计值震荡但真实值平稳这时缩小t_eval的密度或者改用methodLSODA解决刚性数值问题。5. 线性系统大作业避坑清单四条最值得记住的翻车现场5.1 能控性矩阵满秩但place报错先查是不是单位制埋了雷现象ctrb算出来秩等于状态数可ct.place直接抛出“system is not controllable”的异常信息代码逻辑看起来没有疏漏参数和教材一致。原因几乎总是单位制混用转速用 rpm、电压用伏特、转矩用 kg·cm但模型推导里 Kt 和 Ke 却假设为国际单位导致 A 和 B 的量纲不一致数值上伤害了条件数极点配置矩阵变成一个弱可控矩阵。解决方式是把所有物理量先统一到 SI 再重新推导一遍状态空间这比检查代码更高效因为代码只是把错模型喂给算法。这个坑在答辩中问到的概率极高而且当你主动说出“我统一成国际单位后发现算法正常了”评审一般不会继续在这个点上深挖。5.2 仿真发散但极点明明位于左半平面观测器初值才是元凶现象单独仿真状态反馈时一切正常一旦把观测器接进去转速曲线直接飞掉打印极点位置又全部为负实部逻辑上不应该发散。原因观测器初值与真实状态的偏差很大观测器误差在初始阶段衰减过程中产生巨大的瞬态控制量如果系统模型里还有零点或未建模动态这个瞬态会激励出高频分量形成数值发散。解决不要为了展示“观测器收敛能力”而故意设置离谱的初值比如真实转速 100 rad/s 而观测器初值 0控制电压会在第一个采样周期里冲出上百伏把初值设为真实值的 10% 到 20% 内既能看到收敛过程又不会让仿真爆掉。这是仿真和真实系统最大的差别真实系统有物理限幅保护而纯数值积分没有任何保护。5.3 加了积分环节后系统反而超调变大极点重新配置的必要性现象为了消除稳态误差在控制器里加了积分器结果阶跃响应的超调量从 8% 变成 20%完全超出指标。原因积分器改变了系统阶数但你没有把期望极点重新分配等效地把积分器极点放在了一个高增益区域改变了整个闭环的阻尼特性。解决把积分器状态纳入状态向量重新计算新的状态空间系统矩阵并重新用place配置期望极点保证主导极点位置和阻尼比不变积分器极点额外放到左侧不影响主导动态的区域比如取主导极点实部的 4 倍。这个坑背后其实是“动态增阶后必须重新做完整设计流程”而不是在原增益上加一个自由积分项。5.4 仿真结果在答辩图和正文里不一致参数版本控制缺失的教训现象README 里写电机惯量 0.01 kg·m²仿真图用的却是 0.1 的参数结果答辩时被现场操作还原不出来。原因大作业往往在最后一周反复调参不同参数版本的代码被随意覆盖到交稿时根本不知道哪张图对应哪份代码。解决从第一天开始就按参数组维护仿真配置把电机惯量、阻尼系数、期望极点、观测器极点全部集中放在一个config.py或params.yaml里每次跑仿真前先打印参数核验值并给仿真图标注参数版本号。这在作业里看着有点过度工程但实际工作里这是基本习惯——没有参数版本控制的仿真结果本质上都是不可复现的。6. 让大作业从“能交”变“能讲清楚”一个值得做的小验证技巧如果你时间还够我最建议做的一件进阶小事是把状态反馈用的极点与观测器用的极点之间的关系变成一个可以交互调整的参数做一组极点相对位置扫描实验。具体方法是用脚本循环修改观测器极点为反馈极点实部的 2、3、4、6 倍每次跑同一段阶跃响应输出转速曲线、控制电压曲线和观测器误差的 2-范数。你会发现三个结论倍数越大估计状态收敛越快超过 6 倍后控制量的噪声放大非常明显2 倍时观测器动态与控制器动态交织响应曲线会产生一个明显的瞬态低谷。把这组图放进附录答辩时主动说“我扫描了观测器带宽对闭环的影响选 3 倍是兼顾收敛速度与控制噪声”评委不会再质疑你只会背公式。这个技巧背后的工程逻辑是参数扫描——不单看一个工作点而是看设计参数在合理区间内变动的敏感性反映的是你对控制系统设计的理解深度。我自己的习惯是把这类扫描结果和最终选定的参数组合保存在一个独立文件里方便回溯哪个结论是基于什么前提得到的。交完作业之后不妨再做一个小延伸把模型里的粘性摩擦系数改成不确定区间重复极点配置和仿真看闭环稳定性是否还能保住——这项工作正是从“大作业”走向“鲁棒控制的初步认知”最近的一步。希望帮到你。本文还有配套的精品资源点击获取