Simulink倒立摆仿真：从非线性建模到串级PID与LQR控制

简介这是一份面向控制理论与自动化初学者、以及正在学习MATLAB/Simulink仿真的学生的倒立摆仿真资源完整覆盖从非线性动力学建模、控制器设计到Simulink仿真验证的典型流程可用于课程设计、实验报告或科研入门。资源包内共9个文件压缩后仅372KB包含4个MATLAB脚本.m用于编写状态方程与控制律1个Simulink模型文件.slx及对应的r2014a版本、1个.mdl文件另有PPT和Word文档用于讲解原理与记录实验。平台上已有2356人学习浏览说明其可参考价值得到一定验证。内容不仅给出具体模型文件还详细阐述倒立摆的数学模型摆杆角度与角速度、PID/LQR等控制策略、Simulink子系统搭建以及仿真参数设置要点并附带参考文献便于对照学习与二次开发。对于希望在仿真环境中快速上手倒立摆控制的读者按文档步骤即可从搭模型、跑仿真到分析结果一气呵成是理论与实践结合的不错案例。1. 倒立摆仿真为什么总在第一步翻车建模约定比PID参数更先决倒立摆常被当成自动控制课的“hello world”但真要放到MATLAB/Simulink里仿真劝退率并不低许多人花一下午拖模块一按运行就发散最后只能把锅甩给“PID没调好”。我做过的项目里这类翻车多半不在参数而在建模这一层坐标系约定、单位、初始条件、反馈极性任何一个不对后面调三天增益都是白费。这篇文章讲一条不绕弯的路线在Simulink库下先用动力学方程搭出非线性倒立摆模型再用串级PID把杆立住、把小车拉回原点最后用状态反馈和动画回放做一次交叉验证。适合正在做课程设计的新手也适合需要快速验证控制算法、又不想被模型黑匣子坑到的从业者。读完你会得到一整套能复现、能排查、能说服人的仿真流程而不是一个“碰巧能站稳一次”的模型文件。2. 先把物理模型立住从非线性方程到Simulink的建模选型2.1 坐标系与自由度先定义θ0意味着“杆竖在天上”这里讨论的直线倒立摆由一节小车和一根摆杆组成小车在水平导轨上移动摆杆铰接在小车顶端控制输入是作用在小车上的水平力F。系统有两个自由度即小车位移x和摆杆角θ。第一步要做的不是拖模块而是把所有约定写死在注释里θ0必须定义为摆杆严格竖直向上、杆头朝天逆时针为正。很多教材把θ0定义在竖直向下两套约定下仿真结果看起来方向完全相反反馈极性、初始值、控制器符号都得跟着翻这是新手第一个隐蔽坑。单位也要在初始脚本里统一。角度全程用弧度长度用米质量用千克时间用秒力用牛顿。我见过有人在PID里角度用度、积分时间用分钟结果增益调了一个通宵都没找到规律。单位不统一时仿真里的数值就是无意义的数字任何“经验参数”都套不上去。符号含义初版取值单位M小车质量0.5kgm摆杆质量集中质心0.2kgl摆杆质心到铰链距离0.3mg重力加速度9.81m/s²b车体摩擦系数0首版N·s/md关节阻尼系数0首版N·m·s/radF控制力输入由控制器给出Nx小车位移初始0mθ杆角竖直向上为0初始0.1rad2.2 从牛顿方程到Simulink可用形式非线性方程的整理思路摆杆按集中质量处理质心距铰链l转动惯量写为ml²。对课程设计和快速原型这个近似足够。小车水平方向的受力方程为(Mm)ẍ ml·θ̈·cosθ − ml·θ̇²·sinθ F摆杆绕铰链的转动方程为暂不计摩擦阻尼ml²·θ̈ ml·ẍ·cosθ − mgl·sinθ 0第一式里ml·θ̈·cosθ是摆杆质心水平加速度带来的惯性耦合ml·θ̇²·sinθ是向心加速度项第二式三项分别对应角加速度、平动与转动耦合、重力项。首版把b和d都置0是想让模型只剩一个控制输入和两个二阶状态出问题时能一眼定位等控制器能立住杆再往两式右侧补−bẋ和−dθ̇。做仿真模型的优先序是先让结构对再谈参数复杂度。把两式联立消去中间变量得到适合在Simulink里直接计算的两条表达式θ̈ [−F·cosθ − ml·θ̇²·sinθ·cosθ (Mm)·g·sinθ] / [l·(M m·sin²θ)]ẍ [F − ml·θ̈·cosθ ml·θ̇²·sinθ] / (Mm)我习惯先算θ̈再把θ̈代回去求ẍ而不是做二次消元直接得ẍ的闭式。原因是Simulink里信号是“流”两级计算都有独立输出Scope一接就能确认哪一级算错。如果只用一个长表达式写完ẍ后期排查时根本分不清是符号错、单位错还是连错线。2.3 建模路线选型手算方程 vs Simscape物理建模Simulink库下搭倒立摆有两条常见路线。路线A是手推方程用Integrator、Mux、Gain、MATLAB Function这些基础模块搭出来路线B是用Simscape Multibody物理库直接拖一个小车、一个旋转关节、一根摆杆物理库自动解算。两条路线我都跑过倾向很明确第一版用A验证阶段用B。路线A的每个状态都看得见PID反馈接在哪、哪个信号进了控制器出任何现象都能一步步看Scope定位。路线B确实更接近真实机构可视化也漂亮但模型内部的坐标变换和解算对新手几乎是个黑匣子解算器更慢代数环排查也更难。两条路线的参数必须统一比如杆长Simscape里可能填整根杆的长度手算公式里l是质心到铰链的距离两边如果混用曲线就会硬差一个系数。建模路线主要优势主要代价适用场景基础库手算路线A信号完全可见、接线直观公式要自己推模型是理想化近似课程设计、PID/LQR设计验证Simscape Multibody路线B物理细节真实、自带三维可视化解算慢、内部状态不透明机构对标、方案演示、与A互验2.4 线性化与“模型够用”的判断在跑仿真之前建议把线性化模型也准备好。在θ0邻域取sinθ≈θ、cosθ≈1、θ̇²项直接舍去非线性方程变成状态空间形式。取状态向量X[x; ẋ; θ; θ̇]控制量为F得到A [0 1 0 0; 0 0 −m·g/M 0; 0 0 0 1; 0 0 (Mm)·g/(M·l) 0]B [0; 1/M; 0; −1/(M·l)]这个线性模型有两个用途。第一第6章的LQR增益就是由这组A、B算出来的第二它可以充当校验基准在同样的初始角度θ00.1下线性模型和非线性模型前一两秒曲线应几乎重合超过±30°后才明显分叉。如果一上来就分叉线性化大概率写错了。把A、B填进Simulink库的State-Space模块C取单位阵、D取零矩阵就能跑一个线性对照模型。这个“模型够用”的判断标准很实用控制器设计在线性模型上做最终可信度必须在非线性模型上验证。3. 在Simulink里搭出非线性倒立摆模块清单与初始化脚本3.1 最小可运行模型模块列表与信号路线新建一个空白模型建议按“从右到左”的风格搭最左边是控制器反馈中间是动力学计算最右边是积分器链和可视化。最小集合只需要这些模块模块数量作用MATLAB Function1计算θ̈和ẍ两条二阶导Integrator2把θ̈、ẍ分别积出速度和位置Mux1把F、ẋ、θ、θ̇打包成四路输入向量PID Controller2角度内环与位置外环Saturation1限制控制力输出首版设±20NScope2观察角度与位移曲线信号主干是F → Mux → 动力学块 → θ̈、ẍ → 各自积分链 → 反馈。位移链的中间速度ẋ要从第一个积分器输出端引出角度系统的θ̇同样从积分器链中间引出不要用位置信号差分去求速度差分会把数值噪声放大。Simulink里同一信号允许并联多个输出控制器反馈线和示波器线可以同源不会互相干扰。3.2 参数初始化脚本单位与初值都写清楚所有物理参数不要直接写在模块参数里而是放到基础工作区用脚本定义。模块参数一多漏改一个就会变成“不知道角度为什么总是发散的玄学问题”。把下面脚本存成prob_init.m每次改完参数先运行它再启动仿真。%% 倒立摆物理参数国际单位 M 0.5; % 小车质量 kg m 0.2; % 摆杆质量 kg l 0.3; % 摆杆质心到铰链距离 m g 9.81; % 重力加速度 m/s^2 b 0; % 车体摩擦首版先置0 d 0; % 关节阻尼首版先置0 %% 仿真与初始条件 theta0 0.1; % 初始杆角 rad竖直向上为00.1约5.7° x0 0.0; % 初始小车位置 m xd0 0.0; % 初始小车速度 m/s Tsim 5.0; % 仿真时长 stheta0取0.1rad而不是0是有意的初始完全静止在θ0、θ̇0时系统处在数学上的不动点控制器没有误差信号仿真再久也看不出控制器好坏。取0.1让反馈先“干活”又不会偏离线性区太远。b和d先置0的道理在前面讲过等基础结构调稳再加摩擦。3.3 用MATLAB Function封装动力学计算双击Simulink库里的MATLAB Function模块输入下面的函数体。模块对外显示四个输入端和两个输出端对应关系是输入1为控制力F输入2为当前车速xdot输入3为当前杆角theta输入4为当前角速度thetadot。function [xddot, thetaddot] pendulum_dyn(F, xdot, theta, thetadot) % 小车-摆杆非线性动力学 % 输入: F(N), xdot(m/s), theta(rad), thetadot(rad/s) % 输出: 位移二阶导xddot, 角度二阶导thetaddot M 0.5; m 0.2; l 0.3; g 9.81; % 与prob_init.m保持一致 den l * (M m * sin(theta)^2); thetaddot (-F*cos(theta) - m*l*thetadot^2*sin(theta)*cos(theta) ... (Mm)*g*sin(theta)) / den; xddot (F - m*l*thetaddot*cos(theta) m*l*thetadot^2*sin(theta)) / (Mm); end代码逻辑是把2.2的两条表达式逐段翻译。den先算分母避免同一个分母在长表达式里重复出现也方便检查量纲。特别提醒一句MATLAB Function里的物理参数最好不要依赖外部工作区直接把常量写死在函数体内与prob_init.m手动保持一致。有人会把M、m、l也做成输入端口那样更工程化但首版接口多了容易接错我宁可简洁优先。3.4 双积分器链与示波器把状态接进控制器接线顺序建议这样走动力学的第一个输出θ̈接第一节Integrator得到θ̇再接第二节得到θ第二个输出xddot同样接两条积分链得到xdot和x。角度环的反馈接到内环PID的测量端口位置环反馈接到外环PID的测量端口示波器接最终θ和x。一个特别值得养成的习惯跑第一次仿真前把PID增益全部置0让倒立摆在θ00.1rad下自然倒下去。如果θ在零点几秒内单调加速倒向水平说明动力学块、符号、积分链大概率正确如果一开始就乱跳或反向说明有根线接错或者坐标系约定没统一。这是动力学模块最便宜的验收方式不要跳过。提示开环让杆子倒下是正常物理不倒反而说明模型有隐藏约束或初始条件填错。4. 串级PID整定与求解器参数让杆立住、让车回中4.1 控制结构先定串级PD-PID双环为什么不能只控角度只把θ反馈回去、输出力F仿真三秒内θ可能稳住但小车位置会匀速向一个方向漂。原因是单角度环的控制目标是θ0系统对小车位移没有观测作用在小车上的力既要让杆不倒又会推动小车本身力里的任何直流分量都会让小车加速。物理上就是杆立住了、小车撞到导轨尽头照样失败。所以需要串级结构位置外环把x_ref与实际x的误差过PID_x输出θ_ref作为角度环的期望值角度内环把θ_ref与θ之差过PID_θ输出力F。内环必须比外环快一个数量级内环用来快速抵抗重力扰动外环只是慢慢把小车拉回原点。我一般让内环调节时间在0.3s内外环允许2到5s。如果外环输出的θ_ref超过±0.5rad说明位置环增益过大或饱和设置太松模型很容易翻。回路KpKiKd起始说明角度内环8006先从PD起积分留给后续稳态误差处理位置外环401.5外环要慢增益不要追着内环跑控制力饱和±20——防止极端误差下力过大导致步长塌陷4.2 增益起步值与整定顺序上表的数值不是最优解是“起步值”按2.1的物理参数量级估算。换一组参数要从系统特征频率重新估算角度方程线性化后θ̈≈(Mm)g/(Ml)·θ括号里的量就是开环不稳定极点量级Kp要大到能把这个极点拉到左半平面Kd再加阻尼。整定顺序按三步走先断外环把θ_ref直接给常数0只调角度环。Kp_θ从20开始每次×1.5往上加直到θ能在0.3s内回正且无持续振荡如果高频抖加大Kd_θ。再接外环位置环Kp_x从1起步逐步加到小车5s内回中。若发现θ_ref高频摆动优先减小Kd_x或降Kp_x。最后看控制力饱和Scope里观察F曲线如果频繁撞±20N限幅说明控制器用力过猛要降增益或检查模型参数。这条顺序的核心是“先让内环可信再让外环工作”。位置环没接时角度环其实是在一个不稳定的零极点结构上工作只要角度能回正就说明内环已经压住了主线动力学。4.3 求解器与步长设置不同表现下的选择Simulink默认的变步长ode45对大多数连续模型够用但倒立摆有个典型问题角速度可能到每秒数弧度位移却常以厘米每秒量级变化跨轴数值尺度差异大容易触发自适应步长剧烈收缩。求解器这块建议按表配置配置项推荐值适用场景求解器类型变步长连续模型 连续PID求解器名称ode45默认常规验证最大步长1e-3防止大步长跨过振荡峰值相对误差1e-4精度与速度的折中固定步长1e-3把PID改成离散模块之后使用如果仿真越跑越慢先从求解器入手把ode45换成ode15s或ode23t试试。另外当PID Controller里开启离散时间积分器或者直接把PID模块换成离散版本时模型就混入了固定采样此时建议把整个求解器改成固定步长ode4、步长1e-3否则会报采样时间冲突或者曲线出现奇怪的台阶。4.4 判定“立得住”的验收指标仿真跑完别只看“好像稳了”就收工。给一组我常用的验收指标角度稳态误差在±0.01rad以内位置稳态误差在±0.02m以内从θ00.1rad起步调节时间不超过5s控制力不长期停留在饱和限幅角度曲线没有高频毛刺。毛刺明显时优先怀疑仿真输出采样间隔太大或者Kd过大会放大测量噪声。把这组指标写成一个判断表每次改完参数对照打钩。倒立摆仿真最大的错觉就是“曲线没散架控制成功”没有量化指标后面换一组参数或换台机器跑结果可能完全不同。5. 五个高频倒立摆仿真翻车点与排查顺序5.1 代数环报警一运行就报“Algebraic loop”现象点击运行后Simulink提示检测到代数环模型要么卡在第一步要么仿真进度条几乎不动。原因控制器输出的力F经过动力学块和反馈又在某个没有积分器的闭合回路里影响同一个F。最常见的情况是动力学计算里直接引用了刚算出的加速度又用这个加速度参与计算力形成直通回路。解决按“F → 动力学 → 积分速度 → 积分位移 → 反馈 → 控制器”的顺序检查每一段确认所有影响力的路径都经过至少一个Integrator。确实需要直通的场景在反馈支路插入Memory或Unit Delay用一个小延迟打断闭环。代数环有时不报错也能跑但它会让变步长求解器每一步反复迭代表现就是“卡”。5.2 角度初始值设反θ0躺平不动一碰就爆现象把初始角度设为0仿真却“什么都没发生”稍给一个扰动模型立刻发散。θ0看起来就是竖直怎么都不该错。原因很多教材把θ0定义为竖直向下或Y轴方向取反。坐标系约定不一致导致初始状态、重力方向、控制器反馈极性全部错位。解决在prob_init.m注释里显式写清楚“θ0竖直向上逆时针为正”。跑前先画一个初始位形图把小车、摆杆按初值画出来亲眼确认杆子立着再启动仿真。想最快验证反馈极性把θ0设成0.1rad看前0.2s内角度是否回中如果角度被“推离”零点把PID输出取反即可。5.3 角度在±180°附近跳变控制器输出瞬间反向现象角度曲线整体正常偶尔竖着跳一下控制力同时反向模型随后发散。原因matlab里的atan2或wrapToPi返回值都落在(−π, π]角度跨过±π时θ从π直接跳到−π连续变化被显示成阶跃。控制器把这个阶跃当成真实误差输出瞬间反向。解决进PID之前对反馈角度做unwrap或把角度差用wrapToPi包回(−π, π]再计算。如果你的初始θ0超过π/2摆杆已经越过水平线普通PID很难救回来先把θ0约束在±30°以内。这个坑在“摆起控制”里尤其常见做摆动起摆时必须单独处理角度跳变。5.4 步长塌陷仿真越跑越慢像被卡住现象前1秒跑得飞快越往后进度条越慢CPU单核占用拉满仿真时长却几乎不前进。原因控制器饱和或角度接近水平时方程出现大梯度和高频振荡ode45自适应步长被压到极小也可能是5.1的代数环没清干净迭代开销随状态变化暴增。解决先在控制器输出串一个Saturation限幅比如±20N从物理上防止力爆炸再把求解器换成ode15s或ode23t并设置最大步长1e-3。如果仍然慢打开模型设置里的求解器统计信息看平均步长是否小到1e-6量级是就说明某个信号存在高频抖动优先查Kd和控制力饱和而不是继续调求解器。5.5 角度稳住但小车一直漂位置环没接现象杆子立住了小车却整体匀速或缓慢漂离原点直到撞上导轨边界。原因只做了角度反馈系统对x没有闭环控制x对应的极点是积分模式任何直流分量都会让位置持续累积。解决把x_ref0与x做差经过位置环PID后输出θ_ref。位置环增益从Kp_x1起步观察小车是否回中。如果位置环和角度环开始“互相打架”出现缓慢振荡把位置环Kp_x降到1以下让外环比内环慢一个数量级。这个坑在课程设计里出现频率极高但根因不是控制器不行而是控制目标没定义全。6. 用状态反馈与动画回放收尾仿真结果怎么才叫可信6.1 用lqr算增益状态反馈一次到位PID调稳只代表“这组增益在这个模型上有效”不代表你理解了系统能控性。把线性模型交给lqr是让结果可复现、可解释的兜底方案。在prob_init.m的参数基础上构造线性状态空间模型并直接求增益% 状态向量 X [x; xdot; theta; thetadot] A [0 1 0 0; 0 0 -m*g/M 0; 0 0 0 1; 0 0 (Mm)*g/(M*l) 0]; B [0; 1/M; 0; -1/(M*l)]; C eye(4); D zeros(4,1); Q diag([10 1 100 1]); % 角度权重远大于位移 R 0.01; % 控制力成本 K lqr(A,B,Q,R);Q和R的调法有规律可循Q里的角度对应分量给到100左右位移给10左右速度类状态给1让控制器优先保杆R越小控制越“用力”增益越大。求出的K是一个1×4行向量四个分量分别对应x、xdot、theta、thetadot。6.2 把K带入模型用Gain模块实现状态反馈把状态向量用Mux合成四路信号接一个Gain模块增益设为K输出取负后作为F。这样就不需要PID双环了结构更干净。建议在同一个Scope里画PID方案和LQR方案的θ曲线对比权重拉开后LQR的超调和调节时间会更有规律性。再把θ0改到0.3rad重跑一次看哪套方案还能稳住这是检验控制器鲁棒性最省钱的办法。6.3 动画验证与数据导出仿真结果可复现Scope截图不能算验证把数据导出成t、theta、x三个向量再用动画回放才算闭环。导出时注意把仿真数据设为离散输出否则变步长求解器给出的时间点不均匀动画会忽快忽慢。下面这个脚本按固定帧数回放figure; for k 1:20:length(t) xb x(k); th theta(k); xtip xb l*sin(th); ytip l*cos(th); clf; hold on; plot([xb-0.4 xb0.4], [0 0], k, LineWidth, 3); plot([xb xtip], [0 ytip], r, LineWidth, 2); plot(xtip, ytip, ro, MarkerSize, 6); axis equal; xlim([xb-1.5 xb1.5]); ylim([-0.2 1.2]); title(sprintf(t%.2fs, theta%.3f rad, t(k), th)); drawnow; endxlim跟随小车移动是为了防止杆子飞出视野每20个点取一帧减少刷新负担。我自己的习惯是无论PID调得多顺交出去之前一定用LQR加动画核一遍确认不是只能在一个初始角度下成立的脆弱模型。这个习惯救过我很多次让我少走不少冤枉路——仿真能复现、能解释才敢说这步做完了。希望帮到你。本文还有配套的精品资源点击获取