MATLAB控制系统仿真源码包:建模、PID校正与Simulink批量仿真
简介:《MATLAB与控制系统仿真实践》第3版配套源码包(非电子书),面向学习控制工程与系统仿真的高校学生、工程师及自学者,提供与教材章节对应的可运行实例,帮助把根轨迹、频率响应、状态空间等理论落到建模、分析与控制器设计实践中。包内共390个文件,以255个.m脚本和函数为主体,辅以Simulink模型(.mdl、.slx)、.mat数据、.fig图形、说明文档及少量Java、MEX等支撑文件,压缩包约1.41MB,结构便于按知识点检索。已有394人学习。内容覆盖MATLAB数值与符号计算、Simulink动态系统建模、PID/LQR与状态反馈设计、系统辨识、优化调参、实时仿真和代码生成等环节,读者可对照教材复现实验、观察输出响应、理解函数调用与模型搭建思路,也可作为课程实验与自学练习的参考。
1. 从一条阶跃曲线说起:MATLAB 控制系统仿真源码包到底能干什么
课设答辩前一晚,导师一句「把阶跃响应和 Bode 图调出来看看」,比任何算法题都让人清醒。手边这份 MATLAB 与控制系统仿真实践 第 3 版源码,解决的就是这个场景:它不提供电子书正文,给的是能直接在命令行跑起来的 .m 脚本、Simulink 模型和配套数据,从传递函数建模、时域指标计算、频域判稳一路走到 PID 校正,基本覆盖了控制系统课设会用到的那一整套流程。适合三类人:正在做课设、卡在「不知道第一行写什么」的学生;想把零散脚本整理成可复用模板的初级控制工程师;以及需要一套现成例子给新人练手的人。下面按「包里有什么、怎么跑通、怎么改」的顺序拆开讲。
2. 资源构成与技术栈:解压后先认清 .m、.slx 和那几个必须装的工具箱
2.1 目录结构:脚本、模型、数据文件各占哪一块
解开压缩包,第一眼看到的通常不是一个大工程,而是按内容分好的一堆文件夹,每个文件夹对应一个主题——建模、时域分析、频域分析、根轨迹、校正设计、离散化、Simulink 仿真等等。每个文件夹里基本是同一套组合:一到多个 .m 脚本、若干 .mdl 或 .slx 模型、一些 .mat 数据、以及脚本跑出来的 .fig 或 .png 图片。
这种「一个知识点一个文件夹」的组织方式非常适合自学和改造,因为每个脚本都是自洽的,不用先理解整个工程。我的习惯是先把目录列一遍,标出哪些是纯脚本、哪些带模型,因为带模型的脚本对版本和环境最敏感,出问题概率最高。
| 文件类型 | 典型文件名 | 作用 | 注意点 |
|---|---|---|---|
| .m 脚本 | ch03_step.m、bode_demo.m | 建模、分析、绘图的主入口 | 直接双击运行,但依赖当前目录 |
| .m 函数 | calc_spec.m、pid_tune.m | 被脚本调用的工具函数 | 必须在路径上,否则报未定义 |
| .mdl / .slx | motor_ctrl.mdl | Simulink 图形化模型 | 老 .mdl 在新版本里需要另存 |
| .mat | data_sweep.mat | 预置的仿真结果或参数 | 变量名要跟脚本对上 |
| .fig / .png | step_resp.fig | 输出图,便于对照 | 不是运行必需,删掉不影响 |
提示:解压路径里不要出现空格和中文字符。Simulink 的
load_system和部分老脚本对路径编码很敏感,放在D:\work\ctrl_sim\这类纯英文短路径下最省事。
判断一个脚本能不能直接跑,最省时的办法是看它开头有没有clear; close all; clc;和addpath。带addpath(genpath(...))的脚本说明作者考虑过路径问题,可以直接跑;缺少这两个语句的,一般要在它所在目录下运行,否则calc_spec这类本地函数会找不到。
2.2 工具箱依赖:tf、bode、sim 分别归哪个 Toolbox 管
这套代码里出现频率最高的函数,基本都集中在 Control System Toolbox 和 Simulink 两块。tf、ss、zpk、step、impulse、lsim、bode、nyquist、margin、stepinfo、c2d、feedback、series、parallel这些属于 Control System Toolbox;sim、set_param、load_system、Simulink.SimulationInput属于 Simulink;linearize、slLinearizer属于 Simulink Control Design;ilaplace、laplace属于 Symbolic Math Toolbox,只有部分推导脚本会用到,没有也能跑,只是那些脚本会报错。
装没装、有没有授权,跑一段检查代码最快。别靠「我记得装了」,尤其是学校机房,装是装了但 license 过期的情况并不少见。
% 一次性核对这份源码包依赖的工具箱是否可用 need = {'Control System Toolbox', 'Simulink', 'Simulink Control Design'}; info = ver; % ver 返回已安装且已授权的工具箱列表 names = {info.Name}; % 提取名称,方便做字符串比较 for k = 1:numel(need) if any(strcmp(names, need{k})) % 精确匹配,避免子串误判 fprintf('[OK] %s\n', need{k}); else fprintf('[MISS] %s\n', need{k}); end end % 再确认关键函数在路径上,比看工具箱名字更直接 fcn = {'tf', 'stepinfo', 'margin', 'sim'}; for k = 1:numel(fcn) fprintf('%-12s -> %s\n', fcn{k}, which(fcn{k})); endver只列出已安装且授权有效的工具箱,所以它比「文件夹里有没有这个目录」可靠。第二段用which查函数实际所在路径,如果返回空字符串,说明函数根本不在搜索路径上,这时要么是工具箱没装,要么是路径被改乱了。我一般会把这段存成check_env.m放在包根目录,换机器第一件事就是跑它。
还有一个容易忽略的点:pidtune、tunePID这类自动整定函数在较新的 Control System Toolbox 里才有,早期版本没有。如果脚本里用了它们而你的版本偏老,报错是「Undefined function」,而不是提示版本不够,很容易误判成工具箱缺失。
2.3 版本兼容:老 .mdl 与新 .slx 的打开顺序
较老版本保存的模型是 .mdl,新版 MATLAB 打开时会提示「该模型使用旧格式保存」,可以打开,但一旦保存就变成新格式,旧版本再也读不回来。如果这份源码包要来回在机房老版本和自己电脑新版本之间倒腾,正确顺序是:先在新版本里用「另存为」存一份 .slx 副本,原始 .mdl 保留不动,之后所有改动只动副本。
脚本层面的版本差异主要出现在仿真调用上。早期写法是simset配sim,现在推荐Simulink.SimulationInput,后者能一次性把参数、初值、求解器配置打包,还支持parsim并行。两种写法在新版本里都能跑,但simset已被标记为不推荐,长期维护的模板应该往 SimulationInput 靠。
% 判断当前 MATLAB 版本是否支持新版仿真输入对象 if ~verLessThan('matlab', '9.7') % R2019b 起 SimulationInput 比较完善 in = Simulink.SimulationInput('motor_ctrl'); in = in.setVariable('Kp', 3.2); % 覆盖模型工作区里的同名变量 in = in.setModelParameter('StopTime', '8'); out = sim(in); else % 老版本回退到 set_param + sim 的组合 load_system('motor_ctrl'); set_param('motor_ctrl', 'StopTime', '8'); set_param('motor_ctrl', 'Kp', '3.2'); % 直接改模型参数 out = sim('motor_ctrl'); endverLessThan的第一个参数是工具箱名,第二个是版本号字符串,返回真表示当前版本低于给定值。这段分支的意义在于让同一份脚本在机房老版本和本地新版本上都能跑,不用每次手动改。setVariable改的是模型工作区变量,setModelParameter改的是求解器、停止时间这类模型级配置,两者作用域不一样,混用会出问题——这一点在第 4 章会展开。
3. 建模与三大分析:一个二阶系统在代码里被拆了四遍
3.1 三种建模写法:tf、ss、zpk 之间的等价换算
同一个二阶系统,tf、ss、zpk三种写法描述的是同一个对象,只是表达角度不同。传递函数形式直观,适合手算对照;状态空间形式方便做多输入多输出和初值响应;零极点形式一眼能看出阻尼和自然频率。源码包里这三种写法基本都出现了,读懂它们之间的换算,后面看任何脚本都不会卡。
% 同一个二阶系统:s^2 + 0.6s + 1 的分母 num = 1; % 分子系数,按 s 降幂排列 den = [1 0.6 1]; % 分母系数:s^2 系数、s 系数、常数项 sys_tf = tf(num, den); % 传递函数形式 A = [0 1; -1 -0.6]; % 状态矩阵,2x2 对应两个状态变量 B = [0; 1]; % 输入矩阵,列数 = 输入个数 C = [1 0]; % 输出矩阵,行数 = 输出个数 D = 0; % 直接传递项,严格真分式时为 0 sys_ss = ss(A, B, C, D); % 状态空间形式 sys_zpk = zpk(sys_tf); % 从传递函数转到零极点增益形式 disp(sys_zpk) % 直接打印,能看到极点位置和增益 % 反向验证:三种写法转换回来应该完全一致 err = norm(sys_tf - tf(sys_ss)); % 理论上应为 0(数值误差内) fprintf('两种形式的偏差 = %g\n', err);den里的系数是从最高次到常数项依次排列,少写一项结果就完全变了,比如写成[1 0.6]会被当成一阶系统。D项非零意味着系统是「非严格真」的,会出现输入直接跳到输出的情况,阶跃响应在 t=0 那一刻就有值,很多教材例子默认D=0,看到曲线起点不为零先别怀疑代码写错。
状态空间矩阵不是唯一解,同一系统可以有无数种相似变换,但tf(sys_ss)换回来必须和原来的传递函数一致。上面最后一个norm就是干这个用的,我习惯在改造任何模型之前先跑一次这个校验,确认转换链路没问题再往下做。
还有一个坑:用s = tf('s')直接写G = 1/(s^2+0.6*s+1)确实更接近手写公式,但变量名s会被占用成 tf 对象,如果脚本里后面又拿s当循环变量或者拉普拉斯变量,就会翻车。看到别人的脚本里先s = tf('s')然后for s = 1:10,基本可以断定这段跑不通。
3.2 时域指标:stepinfo 一行拿到超调量和调节时间
时域分析部分的核心是step、impulse、lsim三个函数加stepinfo。step算单位阶跃响应,impulse算冲激响应,lsim算任意输入下的响应。真正省事的是stepinfo,把上升时间、调节时间、超调量、峰值、稳态值一次性全算出来,不用自己写公式找最大值。
zeta = 0.3; wn = 2; % 阻尼比与自然频率 sys = tf(wn^2, [1 2*zeta*wn wn^2]); % 标准二阶系统 t = 0:0.01:10; % 时间向量,步长决定曲线平滑度 [y, t] = step(sys, t); % 指定时间向量,避免自动选点疏密不均 S = stepinfo(y, t); % 从响应数据反算指标 fprintf('超调量 = %.2f %%\n', S.Overshoot); fprintf('上升时间 = %.3f s\n', S.RiseTime); fprintf('调节时间 = %.3f s\n', S.SettlingTime); fprintf('峰值时间 = %.3f s\n', S.PeakTime); % 对照理论值,验证参数没写错 sigma = exp(-pi*zeta/sqrt(1-zeta^2)) * 100; fprintf('理论超调 = %.2f %%\n', sigma);| 字段 | 含义 | 常见误读 |
|---|---|---|
| Overshoot | 超调量,百分比 | 不是绝对峰值,是相对稳态值的比例 |
| RiseTime | 上升时间,10%→90% | 老版本默认定义不同,看文档确认 |
| SettlingTime | 调节时间,进入 ±2% 误差带 | 误差带可配置,默认 2% |
| PeakTime | 达到峰值的时间 | 单调系统可能为空或等于末端 |
| SteadyStateValue | 稳态值 | 由最后一个采样点估计,时间太短会偏 |
时间向量t一定要自己指定。不指定时step会自动选点,系统时间常数很小时点数可能只有几十个,画出来的曲线看着像折线;系统时间常数很大时又会生成上万点,stepinfo算得慢。步长取0.01是个比较通用的折中,真正做精度对比时我会再减半跑一次,看指标变化是否超过 1%,超过就说明步长还不够细。
stepinfo还有个重载形式可以直接传系统对象stepinfo(sys),不用自己算t和y。两种写法结果基本一致,但传数据的形式更可控——比如你想对 Simulink 里导出的响应数据算指标,就只能用数据形式,这也为后面第 4 章的模型对拍埋了伏笔。
3.3 频域判稳:margin 的四个返回值怎么读
频域部分主要看bode、nyquist和margin。bode画幅频和相频,nyquist画奈奎斯特图用于判稳,margin直接给出幅值裕度和相位裕度,是校正设计里最常用的一个函数。
s = tf('s'); G = 10 / (s*(s+1)*(s+5)); % 开环传递函数,含一个积分环节 [Gm, Pm, Wcg, Wcp] = margin(G); % 四个输出一次拿全 Gm_dB = 20*log10(Gm); % 幅值裕度换算成 dB 更直观 fprintf('幅值裕度 = %.2f dB @ %.2f rad/s\n', Gm_dB, Wcg); fprintf('相位裕度 = %.2f deg @ %.2f rad/s\n', Pm, Wcp); figure; bode(G); grid on; % 幅频相频一张图,方便看穿越点 title('开环 Bode 图');Gm是绝对倍数不是分贝,很多人第一次打印出来看到Gm = 3.16会以为是 dB 值,其实换算过来才 10 dB。Wcg是相位穿越频率(相角到 -180° 处的频率),Wcp是幅值穿越频率(幅值到 0 dB 处的频率),两个频率经常被搞混,记住Gm对应Wcg、Pm对应Wcp就不会错。
判稳的经验门槛:相位裕度 45° 到 60°、幅值裕度 6 dB 以上,闭环一般既有足够的阻尼又不至于太迟钝。源码包里校正设计的例子基本都是围绕把Pm从十几度拉到 50° 左右展开的。如果margin返回Inf或NaN,通常是没有穿越点——系统要么全程幅值都低于 0 dB,要么是纯积分环节没有相位穿越,不是代码写错了。
nyquist图在实际调试里用得比bode少,但对含右半平面极点的系统,它是唯一能直接判稳的工具。看奈奎斯特图时要数包围 (-1, j0) 点的圈数,nyquist(sys)画出来的曲线方向和包围关系容易看反,稳妥做法是配合margin和闭环阶跃响应一起看,三条证据一致才算过。
3.4 PID 校正:从 s 域拼接到参数扫描
校正设计部分最实用的还是 PID。源码包里的写法有两种:一种是在 s 域直接拼C = Kp + Ki/s + Kd*s,另一种是调pid对象。前者直观,后者方便接pidtune和离散化。
s = tf('s'); G = 10 / (s*(s+1)*(s+5)); Kp = 2; Ki = 1; Kd = 0.3; C = Kp + Ki/s + Kd*s; % 理想 PID 的 s 域写法 L = C * G; % 开环 T = feedback(L, 1); % 单位负反馈闭环 figure; subplot(2,1,1); step(T, 5); grid on; title('闭环阶跃响应'); subplot(2,1,2); bode(L); grid on; title('开环 Bode'); % 参数扫描:看不同 Kp 下超调和调节时间怎么变 Kp_list = 1:0.5:5; for k = 1:numel(Kp_list) Ck = Kp_list(k) + Ki/s + Kd*s; Tk = feedback(Ck*G, 1); Sk = stepinfo(Tk); fprintf('Kp=%.1f 超调=%5.1f%% 调节=%.2fs\n', ... Kp_list(k), Sk.Overshoot, Sk.SettlingTime); endfeedback(L, 1)的第二个参数是反馈通道的传递函数,写1表示单位负反馈;写成feedback(L, 1, +1)就变成正反馈,曲线会一路发散。这个符号错误是新手最常见的翻车点之一,因为函数不报错,只是结果不对。
比例项加大加快响应但抬高超调,积分项消除稳态误差但拖慢响应甚至引入振荡,微分项压制超调但对噪声敏感——这三句话背下来容易,真调起来还是得靠扫描。上面这段循环跑完,基本能看出Kp从 1 加到 5 时超调从十几冲到五十以上的趋势,再回头选一个折中点。
如果版本支持,pidtune(G, 'PID')能直接给一组满足相位裕度要求的初值,比自己从零试凑快很多。但自动整定的结果要拿stepinfo复核一遍,我遇到过整定出来的参数相位裕度很漂亮、调节时间却长得离谱的情况,原因是算法优化的是裕度不是响应速度。
4. Simulink 与脚本联动:把图形模型变成可批量跑的仿真引擎
4.1 参数注入:set_param 改的是哪个层级的参数
Simulink 模型里的参数分几层:模型级(停止时间、求解器类型、最大步长)、模块级(增益值、PID 的 P/I/D)、还有模型工作区变量。set_param能改前两层,格式是set_param('模型名/模块路径', '参数名', '值')。模块路径从模型名开始,中间用斜杠分隔,多层嵌套就写成'motor_ctrl/Controller/PID'。
mdl = 'motor_ctrl'; load_system(mdl); % 加载进内存但不弹出界面 % 模型级:改停止时间和求解器 set_param(mdl, 'StopTime', '10', 'Solver', 'ode45'); % 模块级:改某个 Gain 模块的值 set_param([mdl '/Plant/Gain'], 'Gain', '2.5'); % 模块级:改 PID 模块的三个系数 set_param([mdl '/Controller/PID'], 'P', '2.0'); set_param([mdl '/Controller/PID'], 'I', '1.0'); set_param([mdl '/Controller/PID'], 'D', '0.3'); save_system(mdl); % 需要保留改动时才存盘模块路径写错是最常见的报错来源,提示是「block does not exist」或者「invalid block path」。定位办法是先用find_system(mdl, 'Type', 'Block')把所有模块路径打印出来,复制粘贴过去,别手敲。模型里如果是 Gain 模块,参数名是Gain;如果是 PID 模块,参数名是P、I、D;不同类型模块的参数名不一样,用get_param(路径, 'DialogParameters')能列出全部可改项。
load_system和open_system的区别值得记一下:前者只加载到内存,不出界面,适合脚本批量跑;后者会弹出模型窗口。批量仿真一定用前者,不然跑一百次弹出一百个窗口,机器直接卡死。
4.2 批量仿真:for 循环 + sim + 结果落盘
单次仿真用sim就够了,真正的价值在批量:扫一组参数,把每次的超调量、调节时间存成表,最后画趋势图。这是做课设、写报告、做参数敏感性分析最实用的一段代码。
mdl = 'motor_ctrl'; load_system(mdl); set_param(mdl, 'StopTime', '8', 'SaveOutput', 'on', 'OutputSaveName', 'yout'); Kp_list = 0.5:0.5:5; n = numel(Kp_list); res = zeros(n, 3); % 三列:Kp、超调、调节时间 for k = 1:n set_param([mdl '/Controller/PID'], 'P', num2str(Kp_list(k))); out = sim(mdl); % 返回 SimulationOutput 对象 y = out.get('yout'); % 兼容不同版本,比 out.yout 稳妥 S = stepinfo(y, out.tout); res(k, :) = [Kp_list(k), S.Overshoot, S.SettlingTime]; fprintf('Kp=%.1f -> 超调 %.1f%%, 调节 %.2fs\n', ... Kp_list(k), S.Overshoot, S.SettlingTime); end save('kp_sweep.mat', 'res'); % 落盘,后面画图不用重跑 plot(res(:,1), res(:,2), 'o-'); grid on; xlabel('Kp'); ylabel('超调量 (%)');num2str不能省,set_param的值必须是字符串,直接传数字会报类型错误。out.get('yout')比点语法稳,因为 R2019a 前后SimulationOutput的属性访问方式变过,get方法一直兼容。out.tout是输出时间向量,只有勾选了「保存输出」并且输出端口有接线才有值,模型没接输出端口的话yout是空的,stepinfo会报错说输入维度不匹配。
循环里每次都set_param会有一点点开销,参数组合上千时可以把sim换成parsim走并行。但对几十组参数这种量级,串行完全够用,先保证结果对,再谈快。
4.3 与解析结果对拍:怎么确认模型没接错线
Simulink 模型最大的问题是「看起来对」——线都连上了,仿真也能跑,但接成正反馈、增益接反、或者单位不对,曲线还是会出来一条,只是数值不对。我的习惯是每个模型都用同一组参数,分别跑脚本和 Simulink,把两条曲线叠在一张图上,肉眼能看出来就直接定位问题。
% 解析解:用 Control System Toolbox 算同一个系统 s = tf('s'); G = 2.5 / (s*(s+1)); % 和模型里的 Plant 保持一致 C = 2 + 1/s + 0.3*s; % 和模型里的 PID 保持一致 T = feedback(C*G, 1); [y_ref, t_ref] = step(T, 8); % 仿真解:从 Simulink 里取 load_system('motor_ctrl'); out = sim('motor_ctrl'); y_sim = out.get('yout'); t_sim = out.tout; % 叠图 + 定量误差 figure; plot(t_ref, y_ref, 'b-', 'LineWidth', 1.5); hold on; plot(t_sim, y_sim, 'r--', 'LineWidth', 1.2); legend('解析解', 'Simulink'); grid on; y_sim_i = interp1(t_sim, y_sim, t_ref, 'linear', 'extrap'); fprintf('最大绝对误差 = %.4f\n', max(abs(y_sim_i - y_ref)));interp1那一步是必须的,两个时间向量步长一般不一样,直接相减会维度不匹配。误差在 1e-3 量级以内基本可以认为模型接线正确;差到 0.1 以上,先查反馈符号,再查增益和积分环节的初值,最后查求解器步长。我吃过一次亏:模型里 PID 的 I 项用了离散积分器,脚本里用的是连续1/s,两条曲线差了 20%,折腾了一晚上才定位到是离散化方式不同。
5. 踩坑排查:五个让仿真结果对不上的常见问题
5.1 现象:load_system报「文件不存在」,路径肉眼看着没问题
原因基本是路径里带了中文、空格或特殊符号。Simulink 底层对文件路径的编码处理比 MATLAB 脚本更严格,脚本能cd进去跑,模型加载却会失败。解决方式是把整个源码包挪到纯英文短路径下,比如D:\ctrl_sim\,然后在脚本开头统一用cd(fileparts(mfilename('fullpath')))切到脚本所在目录,不依赖当前工作目录。
5.2 现象:Undefined function 'tf',但明明装过 Control System Toolbox
先跑ver确认工具箱在列表里,再跑which tf看函数路径。两种情况最常见:一是装了但 license 没覆盖,ver里不会出现;二是路径被自己写的同名tf.m覆盖了——很多人练习时在某个文件夹里建过tf.m,之后每次addpath都把它加在前面。用which tf -all能看到所有同名文件,把非工具箱路径的那个删掉或改名即可。
5.3 现象:out.yout取不到数据,报「无此属性」
SimulationOutput的访问方式在不同版本间变过,老版本out.yout直接可用,新版本有时要用out.get('yout'),还有的配置下输出变量名被改成了logsout。稳妥写法是out.get('yout')加一层空值判断:取到空就先检查模型是否勾选了「保存输出」、输出端口有没有接线、OutputSaveName是不是被改过。别急着怀疑代码,先怀疑配置。
5.4 现象:闭环阶跃响应一路发散,越跑越大
feedback的符号接反了,把负反馈写成了正反馈。feedback(L, 1)默认是负反馈,加第三个参数+1就变成正反馈,函数不报错,结果却是发散的。另一种情况是 Simulink 里 Sum 模块的符号列表写成了++或+-接反。排查办法是把开环L的 Bode 图画出来,看相位裕度是否为正——如果margin显示Pm < 0,系统本身就是不稳定的,先调参数再谈曲线。
5.5 现象:仿真跑得极慢,或者曲线有高频抖动
求解器选错。默认ode45适合非刚性系统,遇到含快速动态或刚性环节的模型(比如电机电气时间常数和机械时间常数差两个数量级),步长会被压到极小,仿真时间成倍增长。换成ode15s或ode23tb通常能解决。曲线抖动则多半是最大步长没限制,加一句set_param(mdl, 'MaxStep', '0.01')就好。另外代数环也会导致报错或结果异常,在反馈回路上串一个单位延迟模块能打断代数环。
6. 二次开发:把示例脚本改造成自己的课设模板
源码包里的脚本是「一个知识点一个文件」,跑通没问题,但真要做课设或者工程验证,需要的是「一套参数进去、一张表出来」的批处理能力。我的做法是抽一个统一入口函数,把建模、参数扫描、指标统计、结果落盘全封进去,脚本只负责传参。
function res = sweep_zeta(zeta_list, wn) %SWEEP_ZETA 批量扫描阻尼比对二阶系统阶跃指标的影响 % zeta_list 阻尼比向量,wn 自然频率(标量) t = 0:0.01:15; % 统一时间向量,保证各次结果可比 n = numel(zeta_list); M = zeros(n, 4); % 四列:zeta、超调、调节时间、峰值时间 for k = 1:n sys = tf(wn^2, [1 2*zeta_list(k)*wn wn^2]); y = step(sys, t); % 只取响应值,时间向量已知 S = stepinfo(y, t); M(k, :) = [zeta_list(k), S.Overshoot, S.SettlingTime, S.PeakTime]; end res = array2table(M, 'VariableNames', ... {'zeta', 'Overshoot', 'SettlingTime', 'PeakTime'}); end调用侧就一行res = sweep_zeta(0.1:0.1:0.9, 2);,返回的 table 可以直接writetable成 Excel 交作业,也可以拿res.Overshoot画趋势图。这套结构有三个好处:时间向量统一,各次结果能直接比较;指标计算集中在一处,改公式不用翻十个脚本;返回 table 而不是散落的变量,后续做报告、做图都省事。
改造时的验证方法,我一般分两步走。第一步拿原始脚本里的一组已知参数跑一遍新函数,指标必须完全一致,差一点点都说明循环里的时间向量或参数传递出了问题;第二步故意把wn改成一个极端值(比如 0.1 和 50),看会不会出现时间向量不够长导致SettlingTime截断在末端的情况,出现了就把t的上限改成15/wn自适应。这两步走完,函数才算能进模板库。
还有一个从项目里带出来的习惯:任何批量脚本都先把结果save成.mat,再画图。图可以重画,仿真跑一遍可能是几十分钟。我就遇到过一次,参数扫描跑了四十分钟,画图时手滑关了窗口,中间变量全没了,只能重跑。从那以后我每次写完循环,保存语句一定写在画图语句前面,跑之前先强制走一遍「小样本试跑、确认落盘、再放开全量」的流程。
这份源码包的价值不在于它写得多漂亮,而在于每个知识点都给了能跑的例子,改起来有参照。希望帮到你。
本文还有配套的精品资源,点击获取