简介面向自动控制原理、现代控制理论课程考试及课程设计围绕单倒置摆这一高阶次、多变量、严重不稳定非线性系统完整给出状态空间建模与MATLAB仿真流程。文档从物理模型出发忽略次要因素后推导运动方程选取小车位移、速度及摆角、角速度作为状态变量建立状态空间表达式继而用秩判据证明系统能控通过特征值判断系统不稳定并分别展示全维状态观测器与降维观测器方案采用状态反馈极点配置方法设计控制器附有可直接运行的MATLAB代码与阶跃响应截图。资源为1个PDF文件压缩包仅581KB体量轻但内容紧凑。已有191人学习浏览适合备考、实验仿真或课程设计参考可帮助读者理顺建模推导细节并快速复现仿真结果。1. 单倒置摆控制系统的状态空间建模与MATLAB仿真先解决“能不能控”再谈“怎么控”单倒置摆惯性倒立摆是控制理论里最经典的验证对象它本身是一个不稳定的非线性系统摆杆不施加控制时轻轻一碰就会倒向一侧。所谓“单倒置摆控制系统的状态空间建模”本质是把摆杆角度、角速度、小车位移和速度四个物理量写成一组一阶微分方程组再在平衡点附近做线性化得到一个形如ẋ Ax Bu的连续时间状态空间模型。MATLAB仿真在这个场景里承担两件事一是验证模型是否正确二是验证设计的控制器能否在模型上把摆稳定住。这个标题之所以值得认真拆解是因为状态空间建模的质量直接决定后续LQR、极点配置或MPC控制器的效果。模型里一个摩擦系数取错、一个线性化假设被忽略到仿真阶段就会出现“明明控制器参数很合理但曲线还是发散”的情况。适合读这篇文章的人包括正在做课程设计或毕业设计的自动化/机械专业学生、刚接触现代控制理论想用MATLAB做验证的工程师以及需要快速搭建一个倒立摆仿真demo来验证算法的研究人员。2. 从物理方程到状态空间表达式单倒置摆建模的关键推导与参数选取2.1 单倒置摆的非线性微分方程推导标准的单倒置摆模型由一辆小车和一根固定在车上的摆杆组成。小车在水平方向受外力F驱动摆杆在重力作用下有自然下垂的趋势但控制目标是把摆杆维持在竖直向上的倒置位置。建模时通常取以下物理量小车质量M单位kg、摆杆质量m单位kg、摆杆长度l单位m指质心到转轴的距离、小车位移x单位m、摆杆角度θ单位rad与竖直向上方向的夹角、施加在小车上的水平力F单位N。提示摆杆角度θ的定义方向不同会直接改变后续线性化模型中A矩阵的符号建议在开始写代码之前先在纸上画出坐标系避免MATLAB里矩阵符号对不上。用牛顿第二定律分别对小车和摆杆做受力分析。小车水平方向的力平衡方程是Mẍ F - N其中N是摆杆通过转轴作用在小车上的水平反作用力。摆杆的水平方向运动为m(d²/dt²)(x l·sinθ) N把摆杆的x坐标展开求二阶导数后整理再加上摆杆绕质心的转动方程涉及摆杆转动惯量I最终可以得到两个耦合的非线性微分方程。这个推导过程本身并不复杂但很容易在医院里算错符号。常见做法是引入拉格朗日方程来避免受力分析时方向判断出错拉格朗日量L 动能 - 势能对广义坐标x和θ分别代入欧拉-拉格朗日方程得到的结果与牛顿法完全一致。完整非线性方程组的标准形式为忽略摩擦和摆杆阻尼时(M m)ẍ ml·θ̈·cosθ - ml·θ̇²·sinθ F ml·ẍ·cosθ (I ml²)θ̈ - mgl·sinθ 0其中I是摆杆绕质心的转动惯量对于均质细杆取I (1/12)ml²。这套方程就是后面所有工作的起点也是整个单倒置摆控制系统状态空间建模的核心依据。大部分教材为了让书写简洁会把摆杆当成集中质量点处理此时直接令I 0也可以得到一组可用的简化模型。2.2 在倒置平衡点线性化得到四阶状态空间模型单倒置摆的平衡位置是摆杆竖直向上也就是θ 0或θ π取决于坐标系定义。在竖直向上这个平衡点附近θ和其导数都是小量可以做线性化近似sinθ ≈ θcosθ ≈ 1θ̇²·sinθ ≈ 0把这三条代入非线性方程组得到(M m)ẍ ml·θ̈ F ml·ẍ (I ml²)θ̈ - mgl·θ 0从这两个方程中消去ẍ和θ̈得到分别关于ẍ和θ̈的显式表达式。然后定义状态变量为x₁ x小车位移x₂ ẋ小车速度x₃ θ摆杆角度x₄ θ̇摆杆角速度以u F作为控制输入输出取小车位移x和摆杆角度θ可以写成标准的状态空间形式ẋ₁ x₂ ẋ₂ (-m²l²g / Δ)·x₃ ((I ml²) / Δ)·u ẋ₃ x₄ ẋ₄ (mgl(M m) / Δ)·x₃ (-ml / Δ)·u其中Δ (M m)(I ml²) - m²l²。写成矩阵就是A [0 1 0 0; 0 0 (-m²l²g/Δ) 0; 0 0 0 1; 0 0 (mgl(Mm)/Δ) 0] B [0; (Iml²)/Δ; 0; -ml/Δ]注意A矩阵中与θ相关的项不带任何阻尼项这正是倒立摆系统“自不稳定”的数学体现在没有控制力u时状态方程的特征值会出现在右半平面。这一点在后续MATLAB仿真里会直接表现为阶跃响应发散。2.2.1 状态空间模型参数典型取值设计仿真时参数取值的合理性比数学推导更容易被忽略。常见做法是选择一组能够代表真实实验台架的量级并保持一致。以下是一组在课程设计和仿真验证中非常常用的参数参数符号数值单位小车质量M0.5kg摆杆质量m0.2kg摆杆长度质心距l0.3m重力加速度g9.8m/s²摆杆转动惯量I0.006kg·m²将上述参数代入矩阵算得Δ (0.50.2)×(0.0060.2×0.09) - 0.2²×0.09 0.7×0.024 - 0.0036 0.0132。以后在MATLAB里构造A和B矩阵时直接用这些数值即可不需要每次手算。需要说明的是转动惯量的取值对仿真结果有明显影响。如果手里没有实验数据用均质杆公式(1/12)ml²计算得到0.0015比上表中的0.006小很多此时系统的自然不稳定性会更强控制器需要施加更大的力矩才能稳定。建议做仿真时把I作为一个可调参数写进脚本里方便对比。2.3 验证模型的标准方法代码、特征值与阶跃响应构造好状态空间模型后第一步不是急着设计控制器而是先在MATLAB里验证A矩阵的特征值是否真的落在右半平面。这个验证步骤只需要几个简单命令% 单倒置摆系统参数定义 M 0.5; % 小车质量 kg m 0.2; % 摆杆质量 kg l 0.3; % 摆杆质心到转轴距离 m g 9.8; % 重力加速度 m/s^2 I 0.006; % 摆杆转动惯量 kg*m^2 % 计算分母 Delta Delta (M m)*(I m*l^2) - m^2*l^2; % 状态矩阵 A 和输入矩阵 B A [0 1 0 0; 0 0 -m^2*l^2*g/Delta 0; 0 0 0 1; 0 0 m*g*l*(Mm)/Delta 0]; B [0; (I m*l^2)/Delta; 0; -m*l/Delta]; C eye(4); % 输出全部四个状态变量 D zeros(4,1); % 无直馈 % 检查系统特征值 eig_A eig(A); disp(系统开环特征值); disp(eig_A); % 观察开环阶跃响应这里只观察角度状态 sys_open ss(A, B, C, D); figure(1); step(sys_open, 0.5); grid on; title(单倒置摆开环阶跃响应发散验证);这段代码逻辑非常直接先定义参数计算Δ再按推导出的代数式填入A和B最后用eig函数观察特征值分布。运行后你会看到特征值中至少有一个为正实数大约是5.5左右这对应摆杆角度θ呈指数发散。用step命令画出的曲线中角度状态对应的子图会迅速冲向无穷大这就验证了模型与理论一致。若在特征值结果中看到虚数部分很大或出现正负共轭复根而角度曲线呈现振荡增大而不是单调发散说明线性化时参数可能用错建议核对θ定义方向。3. 能控性分析与LQR控制器设计MATLAB必须做的两组计算3.1 用rank(ctrb(A,B))判定单倒置摆能控性状态空间模型建好之后设计控制器的第一个先决条件是确认系统是能控的。能控性的含义是是否存在一个控制输入u(t)能在有限时间内把任意初始状态驱动到原点。对线性时不变系统能控性判据是能控性矩阵Qc [B AB A²B A³B]满秩即秩等于状态维数n。在MATLAB里这条判断只需要一行代码% 计算能控性矩阵并判断秩 Qc ctrb(A, B); rank_Qc rank(Qc); fprintf(能控性矩阵的秩%d系统维度为4\n, rank_Qc); if rank_Qc 4 disp(系统完全能控可以设计状态反馈控制器); else disp(系统不完全能控需要检查模型参数或执行器位置); end运行结果必然是秩为4因为倒立摆的标准模型在数学上是完全能控的。但有一点值得说明能控性矩阵的条件数condition number可能很大尤其当M和m的比值悬殊时说明系统“接近”不能控状态。这在实际物理台架上有意义——如果小车质量非常大而摆杆质量很小那么靠推动小车来控制摆杆的效率会急剧下降。遇到这种情况可以在MATLAB中用cond(Qc)查看条件数如果超过10⁴建议调整执行器增益或重新考虑建模假设。3.2 用lqr函数整定状态反馈增益K能控性验证通过后最常用的控制方案是线性二次型调节器LQR。LQR的核心思想是寻找一个状态反馈控制律u -Kx使如下性能指标取最小值J ∫(xᵀQx uᵀRu)dt其中Q是状态加权矩阵半正定R是控制输入加权矩阵正定。Q中某个对角线元素越大说明越看重对应状态的快速收敛R越大说明越限制控制力的大小。对单倒置摆而言权重调整的直观逻辑是摆杆角度θ的偏差比小车位移x的偏差更应该被“及时纠正”因为角度误差稍大系统就会倒下。% LQR 控制器设计 Q diag([100, 10, 1000, 100]); % 状态加权矩阵位移、速度、角度、角速度 R 1; % 控制输入加权矩阵 K lqr(A, B, Q, R); disp(LQR 状态反馈增益矩阵 K); disp(K); % 计算闭环系统矩阵 Acl Acl A - B*K; eig_cl eig(Acl); disp(闭环系统特征值); disp(eig_cl);Q矩阵的第四个权重取100而不是更小的值目的是限制角速度过大避免在仿真开始时出现瞬间的剧烈摆动。R取1意味着允许控制力有一个相对宽裕的范围。如果仿真中控制量u饱和或过大把R从1增大到10或100即可K矩阵的整体数值会下降控制更温和但调节时间会变长。这里有一个在MATLAB仿真中容易踩的细节lqr函数默认要求A、B矩阵是连续时间系统的状态空间描述它使用的Riccati方程解法与离散时间的dlqr完全不同不要拿离散化后的矩阵去调用lqr否则算出的K是完全错误的值。提示判断这组Q、R是否合理的经验法则是——闭环特征值的实部绝对值应与系统开环不稳定极点的实部保持在同一数量级。若闭环极点实部过大控制力会超出执行器能力仿真中体现为曲线高频振荡。3.3 闭环系统仿真从阶跃响应曲线判断控制效果设计出K之后下一步是搭建闭环系统并观察状态响应。常见做法是把闭环状态空间模型写成sys_cl ss(Acl, B, C, D)然后使用step或initial命令来观察。区别在于step对应从零状态出发对参考输入的响应而initial对应零输入下对非零初值的响应。对于倒立摆来说initial仿真更有物理意义——模拟用手把摆杆拨开一个小角度然后松开观察控制器能否把它拉回竖直。% 闭环系统阶跃响应 sys_cl ss(Acl, B, C, D); t 0:0.001:5; % 仿真时长5秒步长1ms % 假设初始角度偏离0.1弧度其他状态为零 x0 [0; 0; 0.1; 0]; [y, t, x] initial(sys_cl, x0, t); % 绘制摆杆角度与小车位移响应 figure(2); subplot(2,1,1); plot(t, x(:,3)*180/pi, LineWidth, 1.5); xlabel(时间 (s)); ylabel(摆杆角度 (deg)); title(摆杆角度响应初始偏移0.1rad); grid on; subplot(2,1,2); plot(t, x(:,1), LineWidth, 1.5); xlabel(时间 (s)); ylabel(小车位移 (m)); title(小车位移响应); grid on;运行后理想的响应应该是摆杆角度从5.7°左右快速衰减到0调节时间在2秒以内超调量不超过30%小车位移在初始阶段会先向一侧移动一段距离这是为了“接住”倒下的摆杆然后慢慢回到原点附近。如果角度响应出现持续振荡甚至发散优先检查Q和R的数量级是否差距过大或者从initial换到step观察一下是否有数值异常。另外注意MATLAB中角度单位默认为弧度绘图前必须乘以180/pi换算成度否则曲线看起来会比预期小很多或大很多容易产生误判。4. MATLAB仿真全流程脚本构建、Simulink框图与参数扫描4.1 用脚本文件组织从模型到仿真的完整流程把上面各部分代码整合成一个完整脚本是最多人采用的做法。一个规范的m文件应该包含清晰的区段注释这样后续修改参数时不需要在满是数字的命令行窗口里找历史记录。下面是一个推荐的脚本结构其中包含了从参数定义到结果输出的全过程并增加了一个简单的控制力记录模块% 单倒置摆状态空间建模与LQR控制仿真 % 作者基于经典倒立摆模型整理 % 日期2024-06 clear; clc; close all; %% 1. 参数定义 M 0.5; m 0.2; l 0.3; g 9.8; I 0.006; Delta (M m)*(I m*l^2) - m^2*l^2; %% 2. 连续时间状态空间模型 A [0 1 0 0; 0 0 -m^2*l^2*g/Delta 0; 0 0 0 1; 0 0 m*g*l*(Mm)/Delta 0]; B [0; (I m*l^2)/Delta; 0; -m*l/Delta]; C eye(4); D zeros(4,1); %% 3. 能控性验证 assert(rank(ctrb(A,B)) 4, 系统不完全能控); %% 4. LQR控制器设计 Q diag([100, 10, 1000, 100]); R 1; K lqr(A, B, Q, R); Acl A - B*K; %% 5. 闭环仿真 sys_cl ss(Acl, B, C, D); t 0:0.001:5; x0 [0; 0; 0.1; 0]; % 初始角度偏移 0.1 rad [y, t, x] initial(sys_cl, x0, t); %% 6. 控制力计算与绘图 u -K * x; % 注意这里的维度处理K为1x4x为4xN figure(3); subplot(3,1,1); plot(t, x(:,3)*180/pi, b, LineWidth, 1.5); grid on; ylabel(角度 (deg)); title(单倒置摆LQR控制响应); subplot(3,1,2); plot(t, x(:,1), r, LineWidth, 1.5); grid on; ylabel(位移 (m)); subplot(3,1,3); plot(t, u, k, LineWidth, 1.2); grid on; ylabel(控制力 (N)); xlabel(时间 (s));逻辑说明第3步用assert在脚本开头强制检查能控性如果不满足条件脚本直接报错避免后续计算出无意义的结果。第6步计算控制力时用了u -K * x这里x是一个4行N列的矩阵K是1行4列两者相乘得到1行N列的控制力序列正好对应每个时刻的施加力。在MATLAB中矩阵乘法要求内维匹配K乘x的含义恰是把每个时刻的四个状态分别乘以对应增益并求和这与控制律u -Kx的定义完全一致。4.2 用Simulink搭建单倒置摆仿真框图脚本仿真速度快、便于批量改参数但可视化程度不如Simulink。如果想把状态空间模型和控制器用框图的形式呈现Simulink搭建是更直观的选择。标准做法是使用State-Space模块封装系统模型用Gain模块实现状态反馈K用Sum模块完成u -Kx的运算。打开MATLAB在命令行输入simulink新建空白模型。添加以下模块State-Space模块Simulink/Continuous库设置A为上面计算出的A矩阵B为B矩阵C为eye(4)D为zeros(4,1)初始条件填写[0;0;0.1;0]Gain模块Simulink/Math Operations库增益值填K这里需要以[K]的形式从工作区引入Sum模块设置为“-”即正输入减负输入Scope模块用于观察角度和位移曲线Mux模块把四个状态合成一个向量以便在Scope中显示连线关系是State-Space的输出4维向量分支两路一路进Mux和Scope另一路进Gain模块乘以KGain输出进入Sum的负输入端口Sum的输出接回State-Space的输入端口形成闭环。注意State-Space模块的输入是单路标量控制力u因此Sum模块必须是标量输出。Simulink里这个闭环模型有一个容易犯的错误State-Space模块默认输入端是向量如果你把u定义成1×1标量直接连接没问题但一旦Gain模块的输出维度变成4×1因为K乘以4维状态会得到标量但Gain模块如果设置不当可能输出4×1向量Sum模块就无法正常求差。解决方法是双击Gain模块在“Multiplication”选项中选择Matrix(K*u)或者直接用Fcn模块写表达式。如果不确定建议在连线之间插入Display模块观察维度确认u的维度是1×1。4.3 参数扫描观察Q矩阵对控制效果的影响仿真模型搭好后用一个循环脚本做参数扫描比手动改参数再重跑效率高得多。这里以Q矩阵中角度权重q₃为扫描变量观察超调量和调节时间的变化规律% 参数扫描q3 从 100 扫描到 5000观察响应变化 q3_list [100, 500, 1000, 2000, 5000]; settle_time zeros(size(q3_list)); max_angle zeros(size(q3_list)); for i 1:length(q3_list) Qi diag([100, 10, q3_list(i), 100]); Ki lqr(A, B, Qi, R); Acli A - B*Ki; sys_i ss(Acli, B, C, D); [y_i, t_i, x_i] initial(sys_i, x0, t); % 计算调节时间角度稳定在±2%以内 angle_deg x_i(:,3)*180/pi; idx find(abs(angle_deg) 0.2, 1, first); % 0.2度即约2%的0.1rad if ~isempty(idx) settle_time(i) t_i(idx); else settle_time(i) NaN; end max_angle(i) max(abs(angle_deg)); fprintf(q3%d: 调节时间%.2fs, 最大角度%.2f°\n, ... q3_list(i), settle_time(i), max_angle(i)); end % 绘制结果 figure(4); yyaxis left; plot(q3_list, settle_time, bo-, LineWidth, 1.5); ylabel(调节时间 (s)); yyaxis right; plot(q3_list, max_angle, rs--, LineWidth, 1.5); ylabel(最大摆角 (deg)); xlabel(Q矩阵角度权重 q3); grid on;这段代码用find函数找到角度第一次衰减到0.2°以内的时间点作为调节时间。0.2°是0.1rad初始偏移的约2%符合自动控制原理中调节时间的定义习惯。扫描结果的规律通常是q₃增大时控制器对角度偏差更“敏感”调节时间缩短最大摆角减小但控制力峰值也会同步增大。扫描曲线能从整体上帮助你理解Q矩阵权重与控制性能之间的权衡关系这也是把状态空间建模和MATLAB仿真结合起来最实际的收益。5. 仿真发散排查与模型验证技巧三个必查项和一组验证流程仿真发散、曲线异常、结果与理论不符这些问题在倒立摆仿真中出现频率非常高。本节不列大而全的通查清单只给三个从实践中提炼出的必查项以及一套能在5分钟内完成验证的流程。5.1 必查项一A矩阵中角度项符号与开环极点方向新建一个m文件只做一件事用eig(A)查看开环特征值。理论上单倒置摆模型的四个特征值应当为一个正实根、一个负实根、一对纯虚根或近似如此。正实根对应倒立摆自不稳定的倒下运动负实根对应另一个方向的稳定模态纯虚根来自系统无阻尼自由运动。如果特征值中出现两个正实根或全部为左半平面极点说明A矩阵里与θ相关的项正负号写反了。检查方法很简单手动计算A(2,3)和A(4,3)这两个元素前者应为负数-m²l²g/Δ后者应为正数mgl(Mm)/Δ。符号一正一反恰恰是倒立摆与普通质量-弹簧系统最大的差别也是最容易抄错的地方。5.2 必查项二矩阵乘法维度与initial/step命令的使用initial和step的用法相似但物理含义完全不同。initial模拟的是零输入u0时系统对初始状态的响应因此闭环系统矩阵Acl直接决定响应曲线形态而step模拟的是从零初始状态出发在输入端加一个阶跃信号后的响应。对倒立摆而言如果使用step相当于在平衡状态下突然给小车一个恒定的外力控制器需要同时克服摆杆重力和阶跃外力响应曲线中摆杆角度会经历较大偏移后才回到零。初学者容易交叉使用这两个命令导致曲线不符合预期。另外在矩阵乘法中K*x和x*K得到的结果完全一致前者是1×N向量后者转置后形状不同但数值相同不要因为维度问题在循环里写出u(i) -K*x(:,i)以外的高复杂度写法直接向量化计算更稳妥。提示若仿真输出中角度曲线立刻发散且控制力u始终为0检查Simulink模型中Sum模块的符号和Gain模块的输出维度。这个情况在脚本仿真里很少出现但在Simulink框图里属于高频故障。5.3 必查项三仿真步长与数值稳定性直接使用MATLAB的ode45或Simulink的变步长求解器时倒立摆系统的快速不稳定极点会迫使求解器自动减小步长一般不会因数值原因发散。但如果你把仿真配置改为固定步长比如固定0.01秒就非常容易出现数值发散——因为这个系统的开环正极点决定了特征时间尺度在0.2秒左右0.01秒的步长刚刚能分辨这个动态但遇到LQR高增益控制时闭环系统特征值虚部变大步长0.01秒可能低于快速振荡模态的采样要求。实际处理方式是把固定步长缩小到0.001秒或者改用ode45自适应求解器并在仿真时间轴上均匀插值。5.4 模型验证的标准流程从开环验收到闭环对比做完以上三个必查项推荐按下面的流程把整个模型全面验证一遍这个过程对交付课程报告或项目文档非常有帮助。第一步运行开环eig验证正极点存在第二步用rank(ctrb(A,B))验证能控性第三步用小Q值小R值先做一个保守控制器观察闭环是否稳定第四步增大Q权重观察调节时间是否缩短如果曲线反而震荡加剧说明权重过高、控制力饱和需要减小Q或增大R第五步绘制控制力曲线确认u没有超过执行器限幅。这套流程做完之后如果你还需要做更深入的演示可以在MATLAB中写一个简单的动画用plot(x, y)绘制一个小车与摆杆的几何图形再用drawnow在循环中更新位置与角度。动画代码不过几十行但对直观展示控制效果、让评审或同学一眼看懂你在做什么作用比任何曲线图都大。实现思路是仿真的时间向量t和状态矩阵x已经生成在每个时间点上根据x(i,1)画小车矩形、根据x(i,3)计算摆杆端点坐标画摆杆线段然后用xlim固定坐标轴范围就能看到摆杆从初始偏角被拉回竖直位置的全过程。动画与定量曲线结合呈现在汇报中非常加分这一步只是把已经算好的仿真结果可视化不会改变模型本身的任何行为。本文还有配套的精品资源点击获取