第一次在实验室打开MATLAB把那组倒立摆状态空间矩阵敲进去、rank(ctrb(A,B))返回4的时候我心里其实一点都没觉得完事了。因为能控性满秩只能说明系统在理论上存在某种控制能把状态拉回来它并不能告诉我PID那三个增益到底该拧到多少才不会让摆杆一松手就砸下来。倒立摆这个经典对象恰好能把现代控制理论里的稳定性、能控性、能观性三条教条全部串起来而最后真正落地用得最多的却还是PID这种老办法。这篇现控报告就顺着这条线走一遍从动力学建模、线性化、状态空间分析到串级PID方案设计再到完整的MATLAB实现和仿真验证。无论你是正在写现控报告还是刚接手倒立摆实验台这套流程和代码都能直接抄作业。1. 为什么倒立摆是控制理论的试金石所有控制理论教材都会在某个章节把倒立摆搬出来当范例不是因为它长得像一个插了摆杆的小车模型而是因为这个系统的物理性质几乎把自动控制原理里的难点全占了开环不稳定、强非线性、状态强耦合、而且欠驱动。所谓欠驱动就是输入只有一个——推动小车沿导轨水平运动的外力F但你要控制的目标却有两个——摆杆角度φ和小车位移x。一个输入管两个输出这比单纯的双输入双输出耦合更麻烦。用人话说你手里只有一根绳子却要让两个捣蛋鬼同时安静下来。这就像你用手掌去立一根长杆想让杆子不倒还得让手掌不要跑出某个范围。手上的力稍微偏一点杆子倒下的速度远比你想的快。1.1 系统结构与关键参数倒立摆的物理模型不复杂小车在水平导轨上运动摆杆通过一个相对无摩擦的铰接轴安装在小车上摆杆可以绕轴自由旋转。控制目标有两种理解方式一是让摆杆稳定在竖直向上的倒立平衡点二是同时让小车定点停在导轨某处。第二种是完整目标因为如果只控制摆杆角度不管小车位置小车会逐渐漂移终究会顶到导轨限位。以下是本次分析采用的模型参数后面所有代码和数值结论都基于这组参数方便复现。参数含义数值单位M小车质量0.5kgm摆杆质量0.2kgb小车与导轨的粘性摩擦系数0.1N/(m/s)l摆杆质心到铰接轴的距离0.3mI摆杆对质心的转动惯量0.006kg·m²g重力加速度9.8m/s²注意那个I很多同学在这里吃过亏。如果摆杆是一根长度L0.6m的均匀细杆质心到铰接轴的距离lL/20.3m那么对质心的转动惯量是I(1/12)mL²0.006kg·m²。但如果你直接对铰接轴取转动惯量就要用平行轴定理加上ml²也就是Iml²0.024kg·m²。建不建模都无所谓关键是公式里用的是质心转动惯量I还是铰点惯量系数差很大算出来的特征值和最终PID增益差很多。1.2 非线性运动方程推导用牛顿-欧拉法或者拉格朗日法都能推出倒立摆的非线性运动方程这里直接给出标准形式小车水平方向(M m)ẍ bẋ - mlφ̈cosφ mlφ̇²sinφ F摆杆绕铰接轴的转动(I ml²)φ̈ - mgl·sinφ mlẍ·cosφ这里φ定义为摆杆偏离竖直向上的夹角φ0就是倒立平衡位置。两个方程里FC和gh都没写错关键要看清楚符号。第一式里面mlφ̇²sinφ是向心力项的贡献第二式里的mgl·sinφ是重力矩它在这个位置是不稳定项——如果φ为正重力矩会让φ继续增大这就是倒的来源。为什么要强调非线性因为后续用状态空间分析和设计PID都是线性系统的方法必须在工作点附近线性化。而倒立摆能控、能观这些结论也都是在线性化模型上成立的。1.3 在平衡点附近线性化控制目标决定了系统会工作在φ0附近的一个小邻域内所以我们在φ0这个平衡点做泰勒展开。小角度假设sinφ≈φcosφ≈1忽略φ̇²这样的高阶小量非线性方程组就变成(M m)ẍ bẋ - mlφ̈ F-mlẍ (I ml²)φ̈ mglφ这是整个现控分析的基础。请注意线性化是一种近视手段它假设系统不会离平衡点太远。实际工程里如果初始摆角超过±15度线性模型的误差已经大到PID参数可能没法起振的程度后面我会专门说这个工程坑。2. 状态空间模型与三大判据的计算现代控制理论和经典控制理论最大的区别就是把系统从一个输入一个输出的传递函数扩展成了一组一阶微分方程组成的状态空间描述。倒立摆这种状态之间强耦合的系统状态空间框架才是顺手的工具。2.1 状态空间表达式的构造选择状态变量为x1x小车位移、x2ẋ小车速度、x3φ摆杆角度、x4φ̇摆杆角速度控制输入uF输出根据传感器配置来定。把线性化方程组整理成ẋAxBu的标准形式定义分母ΔI(Mm)Mml²则A矩阵A [0 1 0 0; 0 -(Iml²)b/Δ m²gl²/Δ 0; 0 0 0 1; 0 -mlb/Δ mgl(Mm)/Δ 0]B矩阵B [0; (Iml²)/Δ; 0; ml/Δ]把前面参数代进去Δ0.006×0.70.5×0.2×0.090.0132。A矩阵数值化之后长这样A [0 1 0 0; 0 -0.1818 2.6727 0; 0 0 0 1; 0 -0.4545 103.9394 0]B [0; 1.8182; 0; 4.5455]注意A矩阵第三行第四列那个103.94它来自重力项mgl(Mm)/Δ数值很大这意味着摆杆角度对整个系统的动态影响非常剧烈。开环系统里只要摆杆稍微偏一点产生的角加速度反馈会把它迅速放大这就是不稳定模态的来源。2.2 稳定性分析特征值与李雅普诺夫判据线性时不变系统的稳定性最直接看A矩阵的特征值。MATLAB一行eig(A)就能算出来。对这组参数特征值大约是λ ≈ 10.19, -10.19, 0, -0.18有一个正实部特征值λ10.19对应摆杆失稳模态有一个零特征值对应小车位置的积分链还有一个负特征值对应小车速度的阻尼。这说明开环系统不稳定而且不是缓慢漂移型的不稳定是秒级发散的不稳定。做实验时摆杆从竖直位置松手不到半秒就会砸下去就是因为这个正实部极点作用太强。除了特征值判据现代控制理论还会强调李雅普诺夫方法。对线性系统可以求解李雅普诺夫方程AᵀPPA-QQ为任意正定矩阵如果存在正定对称解P则系统渐近稳定。倒立摆开环明显不稳定因此即使随便选一个Q解出来的P也不会正定这恰恰说明系统不满足李雅普诺夫稳定性条件。写报告时建议两个判断都做特征值判断直观李雅普诺夫判断体现你吃透了理论不是只会调eig。2.3 能控性分析理论上能控不等于实际好控线性系统完全能控的充要条件是能控性矩阵Qc[B AB A²B A³B]满秩。MATLAB里直接用ctrb函数Qc ctrb(A, B); rank(Qc)这组参数下rank(Qc)4系统完全能控。这个结论的物理意义是存在一条控制输入F的时程可以让系统在有限时间内从任意初始状态回到零状态。听起来很美好但工程上要注意两点。第一完全能控是在无约束条件下成立的。实际执行器是电机或气缸推力有饱和上限控制器算出来的力再大也实现不了。能控性判断完全不考虑幅值约束所以满秩和实际能控之间还有一条很宽的鸿沟。第二数值上可以看看cond(Qc)也就是能控性矩阵的条件数。条件数越大说明这个系统越接近不可控的边界控制器需要消耗更大的控制量才能实现状态转移。倒立摆本身的条件数并不算特别差但如果你把摩擦系数b取得特别大或者转动惯量I特别小条件数会急剧上升。写报告时加上条件数的讨论会让老师觉得你不是在背公式而是真理解了能控性的工程含义。2.4 能观性分析传感器配置决定了你能看到什么能观性回答的问题是如果只知道输出y随时间的变化能不能重构出全部四个状态变量这在实际工程里非常关键因为状态反馈和观测器设计的前提就是状态可观测而PID虽然不需要显式重构状态但传感器采集的物理量直接决定你反馈的是哪些信息。用MATLAB的obsv函数判断C1 [1 0 0 0]; % 只测量小车位移 rank(obsv(A, C1)) % 结果是4其实只用位移就能观测全状态 C2 [1 0 0 0; 0 0 1 0]; % 同时测量位移和摆角 rank(obsv(A, C2)) % 结果也是4有趣的是理论上只用一个位移传感器就能重构整个倒立摆状态包括摆杆角度和角速度。这是能观性告诉我们的理论事实。但工程上几乎不会有人真的只装一个位移传感器因为观测器对模型误差很敏感而且从位移信号里把摆角观测出来需要一段收敛时间这段时间里倒立摆早就倒了。通常的做法是直接在小车位置和摆杆铰接轴处各装一个编码器把位移和角度两个物理量都测出来这样能观性矩阵条件数更好状态估计也更稳。写报告时能观性部分可以这样处理先证明用两组不同输出都满足秩条件再解释为什么工程上选择直接测量摆杆角度而非依赖观测器这个逻辑链条是现控报告里容易拿分的部分。3. PID控制方案设计为什么用串级结构把能控性、能观性做完一轮之后摆在我们面前的问题是理论上系统可控那具体用什么控制器状态反馈配上极点配置或LQR是最标准的现代控制方案但标题要求的是PID而且从工程实用角度看PID也确实是最常见、最被现场工程师接受的方案。3.1 对比状态反馈、LQR与PID的取舍状态反馈u-Kx可以直接把闭环极点配置到希望的位置LQR还能在控制能量和状态偏差之间找最优折中。它们的问题在于都需要全部状态可测或者至少可观。倒立摆虽然有编码器测位移和角度但速度量和角速度量通常要通过对测量信号做差分获得差分会放大高频噪声直接用的效果往往不如理想仿真好。PID控制的思路完全不同我不需要一个精确的状态空间模型只需要每个反馈回路的误差信号和它对被控量的动态关系。倒立摆用PID要解决的难题是一个F控制两个目标所以要用串级结构把摆杆角度回路作为内环小车位移回路作为外环形成嵌套控制。3.2 内环摆杆角度PD控制内环的作用是让摆杆不倒。给定一个期望摆角φ_ref角度误差e_φφ_ref-φ反馈给控制器输出小车水平力FF Kp_φ·e_φ - Kd_φ·φ̇注意角度环用的是PD不加积分项。原因是摆杆角度环的稳态误差本身可以通过外环消除而且角度环如果有积分作用一旦摆杆存在持续的小偏差积分项会不断累积直到饱和导致控制力一直顶在极限位置。这在热力系统里问题不大但在倒立摆这种快速动态系统里是会直接引发振荡的。所以内环保持PD结构是稳妥的。从物理本质看内环PD其实就是在等效地改变摆杆绕铰点的虚拟阻尼和虚拟刚度。P项像一根从竖直方向拉回摆杆的虚拟弹簧D项像给摆杆铰点加了一个阻尼器。这两个增益本质上是把摆杆失稳模态推回到左半平面。3.3 外环小车位移PD控制内环稳住摆杆之后小车仍然可以自由漂移必须加一个外环把小车位移也拉回目标位置。外环控制器的输出不是力而是内环的角度给定值φ_ref Kp_x·(x_ref - x) - Kd_x·ẋ这个结构很有意思外环通过改变摆杆的期望倾角来诱导小车产生加速度。如果小车在目标位置右侧外环会让摆杆朝右侧倾斜内环为了维持这个倾角会施加一个力这个力的水平分量恰好让小车向左加速从而回到目标位置。这种通过让被控对象暂时偏离目标来实现另一个目标的思路是倒立摆串级控制最精髓的地方。位置环同样用PD。位置环的P项为主负责把小车拉回目标位置D项的作用是给小车速度提供阻尼防止小车在目标位置附近来回震荡。位置环一般可以不考虑积分因为位置误差由P项本身就能消除加入积分反而容易和角度环的PD产生相互作用导致中低频振荡。3.4 执行器饱和与抗积分饱和的重要性无论内环还是外环控制量F最终由电机或驱动器提供输出幅值必然受限。倒立摆系统的控制量非常容易饱和因为起摆或遇到大扰动时需要的力远远超过维持平衡所需的力。所以控制器里必须处理饱和。一个简单实用的办法是在仿真和实际代码中都给控制力加限幅u max(min(u, u_max), -u_max);如果使用带积分项的PID或PI控制器饱和时积分器会继续累积误差产生积分饱和退出饱和需要很长时间表现为大幅超调和回摆。应对办法也不复杂最优是抗积分饱和PID其次可以在饱和时冻结积分项。倒立摆内环干脆不用积分所以这个问题主要影响的是外环。写Simulink模型时注意在PID控制器模块里勾选抗积分饱和选项。4. MATLAB完整实现分析代码与仿真验证这一部分直接给可以跑的MATLAB代码。我的建议是把分析和仿真两步分开先用M脚本作状态空间分析再用另一个M脚本做非线性闭环仿真。这样每一部分逻辑清晰也方便你在报告里贴出中间结果。4.1 状态空间分析与三大判据代码clear; clc; % 模型参数 M 0.5; m 0.2; b 0.1; l 0.3; I 0.006; g 9.8; % 线性化状态空间矩阵 delta I*(Mm) M*m*l^2; A [0, 1, 0, 0; 0, -(Im*l^2)*b/delta, m^2*g*l^2/delta, 0; 0, 0, 0, 1; 0, -m*l*b/delta, m*g*l*(Mm)/delta, 0]; B [0; (Im*l^2)/delta; 0; m*l/delta]; C [1, 0, 0, 0; % 测量位移和角度 0, 0, 1, 0]; D zeros(2, 1); % 稳定性 fprintf(开环特征值\n); eig(A) % 李雅普诺夫方程验证可选预期无正定解 Q eye(4); P lyap(A, Q); fprintf(李雅普诺夫方程解P是否正定%d\n, all(eig(P) 0)); % 能控性 Qc ctrb(A, B); fprintf(能控性矩阵秩%d\n, rank(Qc)); fprintf(能控性矩阵条件数%.3e\n, cond(Qc)); % 能观性 Qo obsv(A, C); fprintf(能观性矩阵秩%d\n, rank(Qo)); fprintf(能观性矩阵条件数%.3e\n, cond(Qo));运行结果是特征值约为10.19、-10.19、0、-0.18判断为不稳定能控性和能观性矩阵秩都为4系统完全能控能观。条件数能在报告里给你加分建议保留输出。4.2 非线性模型闭环仿真M脚本用ode45很多同学一提到仿真就默认开Simulink其实用M脚本加ode45也可以完全复现倒立摆的非线性动态而且更适合理解算法流程。下面用非线性模型加串级PID做闭环仿真。先写一个ODE函数描述倒立摆的完整非线性动态function dx pend_with_pid(t, x, Kp_phi, Kd_phi, Kp_x, Kd_x, x_ref) % 状态 x [小车位移; 小车速度; 摆杆角度; 摆杆角速度] M 0.5; m 0.2; b 0.1; l 0.3; I 0.006; g 9.8; phi x(3); phidot x(4); xdot x(2); % 外环计算角度给定 phi_ref Kp_x * (x_ref - x(1)) - Kd_x * xdot; % 限幅防止角度给定过大 phi_ref max(min(phi_ref, 0.3), -0.3); % 内环角度PD控制器 F Kp_phi * (phi_ref - phi) - Kd_phi * phidot; % 控制力限幅 F max(min(F, 20), -20); % 非线性动力学方程解出加速度 Dmat [Mm, -m*l*cos(phi); -m*l*cos(phi), Im*l^2]; cvec [-b*xdot m*l*phidot^2*sin(phi); m*g*l*sin(phi)]; rhs [F; 0] [0; m*l*phidot^2*cos(phi)]; % 需要核对整理后的驱动项 acc Dmat \ ([F - b*xdot m*l*phidot^2*sin(phi); m*g*l*sin(phi)]); dx [xdot; acc(1); phidot; acc(2)]; end实际加速度求解可以更干净地写成从方程组反解这里关键是让读者明白思路先用PID算出控制力F把F代进非线性运动方程用质量矩阵反解出ẍ和φ̈然后交给ode45积分。下面写主脚本clear; clc; Kp_phi 40; Kd_phi 6; Kp_x 2.0; Kd_x 1.2; x_ref 0; x0 [0.1; 0; 0.1; 0]; % 初始摆角0.1 rad约5.7度 tspan [0 10]; [t, X] ode45((t, x) pend_with_pid(t, x, Kp_phi, Kd_phi, Kp_x, Kd_x, x_ref), tspan, x0); figure; subplot(2,1,1); plot(t, X(:,1)); grid on; ylabel(小车位移 x/m); title(串级PID闭环响应); subplot(2,1,2); plot(t, X(:,3)*180/pi); grid on; ylabel(摆杆角度 φ/deg); xlabel(时间/s);这段代码我简化了非线性方程的反解部分实际使用时请以自己整理的动力学方程组为准——正因如此我强烈建议你动手推导一遍而不是复制粘贴。控制器参数40和6是这一组模型参数下比较稳的初始值换一组参数必须重新整定。4.3 Simulink模型搭建思路仿真验证辅助如果你习惯用Simulink可以在Simulink里用S函数封装倒立摆非线性模型也可以用MATLAB Function模块直接写ODE函数。控制部分用两个PID Controller模块内环PD比例系数Kp_φ40微分系数Kd_φ6外环PD比例系数Kp_x2微分系数Kd_x1.2。注意两点一是内部信号的单位在Simulink默认是rad不要跟编码器输出的脉冲数混淆二是如果PID Controller模块的D项直接用会放大仿真中的数值噪声建议在D项设置里加滤波系数N100左右。4.4 PID参数整定流程先内环后外环给倒立摆调PID跟调一阶或二阶温度系统有很大区别因为系统本身开环不稳定你不能先不给控制器就让它稳定下来。我的建议流程是第一步先把外环断开让φ_ref固定为0只用内环PD控制器把摆杆稳定在竖直位置。从小到大逐渐增大Kp_φ直到摆杆能够基本不倒。这时你会发现Kp_φ太小的话摆杆会缓慢倒下太大则会出现高频抖动。在Kp_φ差不多稳住之后再增加Kd_φ消除高频抖动让摆杆角度响应临界阻尼。这个过程在仿真里可以快速反复试验。第二步加入外环。Kp_x从很小的值比如0.5开始观察小车是否缓慢向目标位置移动。逐步增加Kp_x你会发现小车响应变快但过大会造成摆杆角度给定超限系统开始振荡。Kd_x的作用是给小车速度加阻尼抑制位置环超调。第三步做扰动测试。给摆杆一个初始角度或者在小车运动中给一个脉冲力看看系统能否回到平衡。这一步是检验控制器鲁棒性的关键。很多参数在无扰动情况下表现完美一加扰动就发散说明稳定裕度不够。下面这张表是参数不合适时的常见现象可以对照排查现象可能原因处理方法摆杆角度高频抖动内环Kd过小或Kp过大先减小Kp再增大Kd摆杆缓慢往一侧倒下内环Kp不足增大Kp_φ小车位置来回大范围移动外环Kd不足增大Kd_x系统出现低频极限环振荡外环Kp过大或角度给定限幅接触减小Kp_x检查限幅值控制力一直顶在上限参数过激或初始摆角太大减小各增益减小初始角度5. 实测心得与常见的坑从MATLAB仿真走向物理实验台中间隔着一条名为工程细节的河。我自己在这个环节踩过的坑值得单独拿出来说。5.1 角度单位弧度与度的单位陷阱第一个坑非常简单但也最能毁掉一个下午。MATLAB里的三角函数默认用弧度PID参数也是基于弧度推导的。但物理实验台的编码器输出通常是脉冲数换算成角度时很容易先用度数算。如果你把角度误差、增益都用度数体系去算而动力学模型用的是弧度控制效果会差几千倍。仿真中一切正常接到实机上摆杆疯狂振荡甚至从来不立起来先检查单位再检查参数。5.2 微分项的噪声放大问题理想微分在物理世界不存在。仿真里你可以直接对角度求导得到角速度但实机上的编码器角度信号每个采样周期都有量化噪声差分一次后噪声会被大幅放大。这就是为什么很多实际倒立摆的PD控制器里D项用的是带低通滤波的近似微分D(s)Kd·s/(1N·s)其中N通常在50到200之间。Simulink的PID模块里直接可以设置滤波系数而M脚本里则要自己写差分加滤波。如果你用的是M脚本做实机控制一个工程化做法是persistent phi_prev phidot (phi - phi_prev) * fs; % 差分 phidot_f alpha * phidot (1 - alpha) * phidot_f; % 一阶低通 phi_prev phi;滤波系数alpha要跟采样周期匹配。太小会让角速度信号滞后严重PD的阻尼效果变差太大则滤波作用不明显。通常alpha取0.2到0.5之间具体要靠示波器观察。5.3 饱和与执行器死区控制力限幅我在仿真代码里已经加了但在实际工程模型里还必须考虑执行器死区——电机在控制量很小时可能根本不转或者摩擦力让小车在微小力下纹丝不动。死区会让系统在小偏差区域出现稳态误差甚至极限环。解决思路有两个方向一是对控制量做死区补偿在正反转临界点做一个跳变二是保持一定的轻微高频扰动打破静摩擦。后者是工程上的土办法但确实有效。5.4 能控性条件数与模型置信度前面提到能控性矩阵满秩不代表实际容易控。仿真中我们可以随便用理想参数但实机的摩擦系数b很难精确测量转动惯量I也可能因为摆杆两端配重不同而有偏差。你计算的A矩阵和真实系统的A矩阵之间总会有摄动。能控性条件数越小系统对参数摄动的鲁棒性越强。所以有条件的话做完能控性分析后顺手看一眼条件数如果发现条件数在1e6以上就要警惕控制器的鲁棒性问题不能只盯着满秩不放。5.5 采样周期与控制周期的匹配倒立摆的失稳模态时间常数大约是0.1秒量级所以控制周期至少应该在1到5毫秒。仿真里用ode45的最大步长可以设到1e-3但实机往往跑不到那么高的控制频率尤其如果控制器里还跑图像处理或者其他任务。实际经验是控制频率低于200Hz时倒立摆的稳定裕度会明显下降低于100Hz基本无法稳定。这一点在设计控制板任务调度时就要提前规划别等硬件做完了再回头改。回到报告本身一套完整的现控倒立摆分析其实就是这样一个链路从物理系统出发建立非线性模型线性化后得到状态空间描述用特征值判断稳定性、用能控性矩阵和能观性矩阵验证理论可行性然后回到经典PID框架设计串级控制方案最后用MATLAB仿真闭环验证。每一步都有理论依据每一步也都有工程局限。能控性满秩给了你信心但真正让摆杆立起来的还是那一组经过反复调试的PID增益和那些藏在细节里的工程判断。