简介本资源是一套面向控制工程专业本科生与初学者的球杆系统建模与稳定性分析实践材料聚焦线性系统理论在典型机电装置中的建模、仿真与稳定性判据应用。资源包含3个核心文件2个MATLAB源码文件极点配置.m、判断系统能控能观.m用于实现状态空间建模、能控能观性验证及闭环极点配置1份Word文档球杆系统建模分析.docx系统梳理了牛顿-欧拉法与拉格朗日法建模思路、特征值判稳原理及Simulink仿真要点。压缩包为RAR格式总大小605KB轻量易下载结构紧凑便于快速上手。已有1442人学习下载适合课程设计、控制原理实验及考研复试项目复现。读者可直接运行代码验证理论结论结合文档理解从物理建模→数学推导→MATLAB实现→稳定性分析的完整技术链路掌握控制系统建模与分析的关键实践能力。1. 球杆系统不是玩具而是线性系统理论的“压力测试仪”球杆系统常被误认为是教学演示用的简化模型——一根杆绕固定点摆动、小球沿杆滑动看似结构简单。但实际在控制工程实践中它是一个典型的非最小相位、欠驱动、强耦合二自由度系统小球位置与杆角速度相互牵制输入如电机扭矩无法独立控制两个状态变量且存在自然不稳定平衡点竖直向上。这种特性使它成为检验控制器鲁棒性、验证能控/能观性判据、测试极点配置效果的黄金标尺。MATLAB 不是仅用来画图或跑仿真而是通过 symbolic math toolbox 推导解析动力学方程、用 Control System Toolbox 构建状态空间模型、借助 Robust Control Toolbox 分析摄动影响。适合正在啃《现代控制理论》教材的本科生、准备数学建模国赛中控制类题目的参赛队以及需要快速验证 PID 或 LQR 设计流程的现场工程师——你调参时看到的阶跃响应超调、稳态误差、甚至仿真发散背后全是拉格朗日方程里漏掉的科氏力项或能观性矩阵秩不足的真实反馈。2. 从物理约束出发推导状态空间模型拉格朗日法 MATLAB 符号计算2.1 为什么必须用拉格朗日方程而非牛顿第二定律球杆系统含两个广义坐标杆转角 θrad和小球沿杆位置 xm二者运动相互耦合。若强行用牛顿法需分解约束反力杆对球的法向力、摩擦力、引入额外未知量方程数膨胀且易出错。而拉格朗日方程直接基于能量$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right) - \frac{\partial L}{\partial q_i} Q_i$$其中 $L T - V$ 为拉格朗日函数$T$ 是系统总动能$V$ 是势能$Q_i$ 是非保守力广义力。该方法自动消去理想约束力仅保留输入扭矩 τ 作为广义力天然适配多自由度机械系统。MATLAB 的 Symbolic Math Toolbox 可符号化推导全部偏导项避免手算微分错误——这是本项目球杆系统建模分析.docx中明确强调却常被跳过的前提。2.2 在 MATLAB 中完成符号建模的四步实操2.2.1 定义符号变量与物理参数syms theta(t) x(t) tau(t) real syms m M g l J b_k b_x % 小球质量、杆质量、重力加速度、杆长、转动惯量、杆阻尼系数、球滑动阻尼系数 % 注意J 是杆绕支点的转动惯量若为均质细杆则 J (1/3)*M*l^2 Dtheta diff(theta, t); Dx diff(x, t); D2theta diff(theta, t, 2); D2x diff(x, t, 2);提示real属性确保后续simplify不引入复数共轭b_k和b_x必须显式定义否则默认无阻尼导致能控性矩阵奇异——这是.rar包中判断系统能空能观.m运行失败的常见原因。2.2.2 构建动能 T 与势能 V 表达式% 杆动能转动 小球动能平动转动此处按质点处理忽略球自转 T_rod (1/2)*J*Dtheta^2; T_ball (1/2)*m*(Dx^2 (x*Dtheta)^2 2*x*Dx*Dtheta*cos(theta)); % 含科氏力项 T T_rod T_ball; % 势能杆重心高度 小球高度以支点为零势能点 V_rod M*g*(l/2)*cos(theta); V_ball m*g*x*cos(theta); V V_rod V_ball; L T - V;注意T_ball中2*x*Dx*Dtheta*cos(theta)是科氏力对应的动能交叉项漏掉此项将导致线性化后 A 矩阵缺失关键耦合元素后续极点配置失效。球杆系统建模分析.docx第 3.2 节明确指出此为高频错误点。2.2.3 符号求解运动微分方程% 拉格朗日方程对 theta 和 x 分别求导 eq1 diff(diff(L, Dtheta), t) - diff(L, theta) tau - b_k*Dtheta; % 杆方程含输入τ和阻尼 eq2 diff(diff(L, Dx), t) - diff(L, x) -b_x*Dx; % 小球方程无直接输入仅受摩擦 % 解出二阶导数 D2theta, D2x sol solve([eq1, eq2], [D2theta, D2x]); D2theta_eq simplify(sol.D2theta); D2x_eq simplify(sol.D2x);此时D2theta_eq和D2x_eq是含theta,x,Dtheta,Dx,tau的非线性表达式即系统原始动力学方程。2.2.4 线性化并生成状态空间模型% 在平衡点 [theta0, x0, Dtheta0, Dx0] 处线性化小角度近似 A_sym jacobian([Dtheta; Dx; D2theta_eq; D2x_eq], [theta; x; Dtheta; Dx]); B_sym jacobian([Dtheta; Dx; D2theta_eq; D2x_eq], tau); A0 double(subs(A_sym, {theta,x,Dtheta,Dx,tau}, {0,0,0,0,0})); B0 double(subs(B_sym, {theta,x,Dtheta,Dx,tau}, {0,0,0,0,0})); % 构建 ss 对象注意状态顺序[theta; x; Dtheta; Dx] sys_lin ss(A0, B0, [1 0 0 0], 0); % 输出为 theta杆角关键参数说明A0是 4×4 系统矩阵其特征值决定开环稳定性B0是 4×1 输入矩阵反映扭矩对各状态的影响权重。.rar中极点配置.m直接读取此sys_lin进行设计若此处线性化点选错如选 π 而非 0后续所有控制律将失效。3. 能控性与能观性验证不只是秩判据更是控制器部署的准入门槛3.1 能控性矩阵的物理意义与秩缺陷诊断球杆系统的能控性本质是能否通过单一输入扭矩 τ在有限时间内将系统从任意初始状态驱动到原点θ0, x0, ω0, v0数学上由能控性矩阵 $ \mathcal{C} [B\ AB\ A^2B\ A^3B] $ 的秩判定。但单纯rank(C)4不够——需检查条件数cond(C)。若cond(C) 1e6说明矩阵接近奇异数值计算中微小扰动会导致控制增益剧烈震荡实际硬件执行时电机易饱和。3.1.1 在 MATLAB 中执行完整能控性分析C_mat ctrb(A0, B0); fprintf(能控性矩阵秩: %d\n, rank(C_mat)); fprintf(能控性矩阵条件数: %.2e\n, cond(C_mat)); % 若条件数过大检查是否漏掉阻尼项 b_k 或 b_x if cond(C_mat) 1e5 warning(能控性矩阵病态请检查 b_k, b_x 是否为0或过小); end % 可视化能控性 Gramian李雅普诺夫方程解 Qc lyap(A0, -C_mat*C_mat); % 近似能控性 Gramian eig_Qc eig(Qc); fprintf(能控性 Gramian 特征值: [%.3f, %.3f, %.3f, %.3f]\n, eig_Qc);注意lyap(A0, -C_mat*C_mat)求解 $A_0 P P A_0^T -C C^T$其特征值反映各状态方向上的能控能量。若某特征值接近零如 1e-8对应状态几乎不可控——这在球杆系统中常表现为小球位置 x 的能控性远弱于杆角 θ需在控制器中降低 x 通道权重。3.2 能观性矩阵与传感器布置的强关联能观性回答仅测量杆角 θ输出 yθ能否唯一重构全部状态 [θ, x, ω, v]能观性矩阵 $ \mathcal{O} [C; CA; CA^2; CA^3]^T $ 的秩必须为 4。但现实中若只装一个编码器测 θC[1 0 0 0]则rank(O)通常为 3 —— 小球位置 x 成为不可观状态。.rar中判断系统能空能观.m的核心逻辑即在此。3.2.1 用 MATLAB 验证不同传感器配置的效果% 方案1仅测杆角 C1 [1 0 0 0]; O1 obsv(A0, C1); fprintf(仅测θ时能观性矩阵秩: %d\n, rank(O1)); % 通常为3 % 方案2增加小球位置传感器如直线电位器 C2 [1 1 0 0]; % 同时测θ和x O2 obsv(A0, C2); fprintf(测θx时能观性矩阵秩: %d\n, rank(O2)); % 应为4 % 方案3测θ和杆角速度ω编码器带速反馈 C3 [1 0 1 0]; O3 obsv(A0, C3); fprintf(测θω时能观性矩阵秩: %d\n, rank(O3)); % 验证是否足够提示球杆系统建模分析.docx第 4.1 节指出方案2虽能观但x传感器噪声大方案3更实用因ω可由θ微分获得需加低通滤波。.rar中未提供滤波代码实际部署时需在极点配置.m前插入y_filtered filter([1 0.9], [1 -0.9], y_raw);。4. 极点配置实现LQR控制器从理论公式到可部署的离散化代码4.1 为什么极点配置比PID更适合球杆系统PID 控制器在球杆系统中面临根本局限它本质是单输入单输出SISO设计而球杆系统是多输入多输出MIMO强耦合对象。当小球偏离杆中心时仅调节杆角无法快速抑制 x 振荡需同时协调 θ 与 x 的动态响应。极点配置通过状态反馈 $ u -Kx $直接将闭环极点置于期望位置实现多变量协同控制。MATLAB 的place()函数可精确配置 4 个极点但需满足能控性前提——这正是第 3 章验证的必要性。4.1.1 手动选择极点的工程准则极点类型期望位置物理意义球杆系统典型值主主导极点$-5 \pm 3j$决定整体响应速度与超调σ5 对应调节时间 ~0.8s快速衰减极点$-20$抑制高频振荡提升鲁棒性避免电机带宽限制零极点对消$-10$抵消系统右半平面零点若存在球杆系统通常无RHP零点注意.rar中极点配置.m默认使用p [-53j, -5-3j, -20, -10]但若你的杆长 l 增大需同比例减小实部如 l 加倍则 σ 减半否则控制器过度激进导致电机饱和。4.1.2 生成状态反馈增益 K 并验证闭环性能p [-53j, -5-3j, -20, -10]; % 期望极点 K place(A0, B0, p); % 计算反馈增益 fprintf(状态反馈增益 K [%.3f, %.3f, %.3f, %.3f]\n, K); % 构建闭环系统 A_cl A0 - B0*K; sys_cl ss(A_cl, B0, [1 0 0 0], 0); figure; step(sys_cl); title(闭环阶跃响应); % 检查闭环极点是否匹配 eig_cl eig(A_cl); fprintf(实际闭环极点: \n); disp(eig_cl);关键参数说明K是 1×4 行向量对应[k_theta, k_x, k_omega, k_v]。若eig_cl与p偏差 0.1说明place()数值不稳定应改用acker()或手动构造K (place(A0,B0,p))。4.2 从连续到离散嵌入式部署前的采样周期选择MATLAB 默认设计连续控制器但实际 DSP 或 STM32 需离散化。采样周期 $T_s$ 选择不当会导致性能恶化$T_s$ 过大10ms离散化引入相位滞后控制器响应迟钝$T_s$ 过小0.1msCPU 负载过高且 ADC 采样噪声放大。4.2.1 使用 c2d() 进行零阶保持离散化Ts 0.005; % 200Hz 采样率兼顾响应与负载 sys_d c2d(sys_cl, Ts, zoh); % 零阶保持离散化 % 提取离散状态空间矩阵 Ad sys_d.A; Bd sys_d.B; Cd sys_d.C; % 生成可用于 C 语言移植的差分方程 % x(k1) Ad*x(k) Bd*u(k) % y(k) Cd*x(k) fprintf(离散化后 Ad \n); disp(Ad); fprintf(离散化后 Bd \n); disp(Bd);提示.rar中未提供离散化代码但极点配置.m输出的K是连续域增益。若直接用于离散系统需同步离散化KKd K * (eye(4) - Ad)\Bd见《Digital Control of Dynamic Systems》第 4.3 节否则实际控制效果严重劣化。5. 稳定性边界验证用 Lyapunov 函数量化鲁棒裕度5.1 为什么特征值判据在非线性系统中不够用球杆系统的原始模型是非线性的含 sinθ, cosθ, x·ω² 项线性化仅在平衡点附近有效。当小球大幅滑动或杆角超过 ±15°线性控制器可能失稳。Lyapunov 稳定性理论提供全局/局部稳定域估计若存在正定函数 $V(x)0$ 且 $\dot{V}(x)0$则系统在该区域内渐近稳定。MATLAB 的lyap()可求解李雅普诺夫方程 $A^TP PA -Q$其中 $P$ 定义椭球形稳定域 $x^TPx c$。5.1.1 计算最大不变椭球稳定域Q eye(4); % 权重矩阵通常取单位阵 P lyap(A_cl, -Q); % 解 A_clP P*A_cl -Q % 计算稳定域半径 c需满足 x^TPx c 时闭环稳定 % 通过仿真找到最大 c 使得所有轨迹收敛 c_candidates logspace(-2, 1, 50); c_max 0; for c c_candidates % 初始化状态在椭球边界 x0 sqrt(c)*inv(chol(P))*randn(4,1) x0 sqrt(c) * inv(chol(P)) * randn(4,1); [~, y] ode45((t,x) (A_cl - B0*K)*x, [0 5], x0); if max(abs(y(:,1))) 0.1 max(abs(y(:,2))) 0.1 % θ和x均收敛 c_max c; end end fprintf(估计最大稳定域半径 c %.3f\n, c_max); % 可视化稳定域投影到θ-x平面 theta_grid linspace(-0.5, 0.5, 100); x_grid linspace(-0.3, 0.3, 100); [THETA, X] meshgrid(theta_grid, x_grid); V_vals zeros(size(THETA)); for i 1:length(theta_grid) for j 1:length(x_grid) x_vec [THETA(j,i); X(j,i); 0; 0]; % 初始速度为0 V_vals(j,i) x_vec * P * x_vec; end end contour(THETA, X, V_vals, [c_max c_max], LineWidth, 2, Color, r); title(Lyapunov 稳定域θ-x平面投影); xlabel(\theta (rad)); ylabel(x (m));注意此代码计算的是线性闭环系统的稳定域但实际非线性系统稳定域更小。.rar中未包含此验证而球杆系统建模分析.docx第 5.3 节强调若要求小球初始位置 |x₀| 0.15m则需将c_max设为 0.02 并重新设计 K。5.2 用 μ-分析评估参数摄动鲁棒性真实系统中杆长 l、小球质量 m 存在制造公差±5%阻尼系数 b_k 会随温度变化±30%。μ-分析结构奇异值量化系统在这些摄动下保持稳定的最大允许不确定性。MATLAB Robust Control Toolbox 提供musyn()和mussv()函数。5.2.1 构建摄动模型并计算鲁棒稳定裕度% 定义摄动块delta_l (l 变化), delta_m (m 变化), delta_bk (b_k 变化) delta_l ultidyn(delta_l, [1 1], Bound, 0.05); delta_m ultidyn(delta_m, [1 1], Bound, 0.05); delta_bk ultidyn(delta_bk, [1 1], Bound, 0.3); % 将摄动注入 A0 矩阵示例l 影响 J 和重力项 A_perturbed A0 delta_l*Adl delta_m*Adm delta_bk*Adbk; % Adl, Adm, Adbk 为各参数的灵敏度矩阵需从符号模型导出 % 构建不确定系统 sys_unc ss(A_perturbed, B0, [1 0 0 0], 0); % 计算结构奇异值 [mu_val, mu_frequencies] mussv(sys_unc, [], m); fprintf(最大结构奇异值 μ %.3f\n, max(mu_val)); if max(mu_val) 1 fprintf(系统对 ±5%% 参数摄动鲁棒稳定\n); else fprintf(系统在摄动下可能失稳请加强控制器鲁棒性\n); end提示.rar中未实现 μ-分析但球杆系统建模分析.docx第 6.2 节指出当 μ 0.8 时建议在 LQR 中增大 Q 矩阵中状态权重如Q diag([10, 5, 1, 1])以提升对参数变化的容忍度。本文还有配套的精品资源点击获取 球杆系统建模与控制：MATLAB符号推导与鲁棒极点配置
简介本资源是一套面向控制工程专业本科生与初学者的球杆系统建模与稳定性分析实践材料聚焦线性系统理论在典型机电装置中的建模、仿真与稳定性判据应用。资源包含3个核心文件2个MATLAB源码文件极点配置.m、判断系统能控能观.m用于实现状态空间建模、能控能观性验证及闭环极点配置1份Word文档球杆系统建模分析.docx系统梳理了牛顿-欧拉法与拉格朗日法建模思路、特征值判稳原理及Simulink仿真要点。压缩包为RAR格式总大小605KB轻量易下载结构紧凑便于快速上手。已有1442人学习下载适合课程设计、控制原理实验及考研复试项目复现。读者可直接运行代码验证理论结论结合文档理解从物理建模→数学推导→MATLAB实现→稳定性分析的完整技术链路掌握控制系统建模与分析的关键实践能力。1. 球杆系统不是玩具而是线性系统理论的“压力测试仪”球杆系统常被误认为是教学演示用的简化模型——一根杆绕固定点摆动、小球沿杆滑动看似结构简单。但实际在控制工程实践中它是一个典型的非最小相位、欠驱动、强耦合二自由度系统小球位置与杆角速度相互牵制输入如电机扭矩无法独立控制两个状态变量且存在自然不稳定平衡点竖直向上。这种特性使它成为检验控制器鲁棒性、验证能控/能观性判据、测试极点配置效果的黄金标尺。MATLAB 不是仅用来画图或跑仿真而是通过 symbolic math toolbox 推导解析动力学方程、用 Control System Toolbox 构建状态空间模型、借助 Robust Control Toolbox 分析摄动影响。适合正在啃《现代控制理论》教材的本科生、准备数学建模国赛中控制类题目的参赛队以及需要快速验证 PID 或 LQR 设计流程的现场工程师——你调参时看到的阶跃响应超调、稳态误差、甚至仿真发散背后全是拉格朗日方程里漏掉的科氏力项或能观性矩阵秩不足的真实反馈。2. 从物理约束出发推导状态空间模型拉格朗日法 MATLAB 符号计算2.1 为什么必须用拉格朗日方程而非牛顿第二定律球杆系统含两个广义坐标杆转角 θrad和小球沿杆位置 xm二者运动相互耦合。若强行用牛顿法需分解约束反力杆对球的法向力、摩擦力、引入额外未知量方程数膨胀且易出错。而拉格朗日方程直接基于能量$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right) - \frac{\partial L}{\partial q_i} Q_i$$其中 $L T - V$ 为拉格朗日函数$T$ 是系统总动能$V$ 是势能$Q_i$ 是非保守力广义力。该方法自动消去理想约束力仅保留输入扭矩 τ 作为广义力天然适配多自由度机械系统。MATLAB 的 Symbolic Math Toolbox 可符号化推导全部偏导项避免手算微分错误——这是本项目球杆系统建模分析.docx中明确强调却常被跳过的前提。2.2 在 MATLAB 中完成符号建模的四步实操2.2.1 定义符号变量与物理参数syms theta(t) x(t) tau(t) real syms m M g l J b_k b_x % 小球质量、杆质量、重力加速度、杆长、转动惯量、杆阻尼系数、球滑动阻尼系数 % 注意J 是杆绕支点的转动惯量若为均质细杆则 J (1/3)*M*l^2 Dtheta diff(theta, t); Dx diff(x, t); D2theta diff(theta, t, 2); D2x diff(x, t, 2);提示real属性确保后续simplify不引入复数共轭b_k和b_x必须显式定义否则默认无阻尼导致能控性矩阵奇异——这是.rar包中判断系统能空能观.m运行失败的常见原因。2.2.2 构建动能 T 与势能 V 表达式% 杆动能转动 小球动能平动转动此处按质点处理忽略球自转 T_rod (1/2)*J*Dtheta^2; T_ball (1/2)*m*(Dx^2 (x*Dtheta)^2 2*x*Dx*Dtheta*cos(theta)); % 含科氏力项 T T_rod T_ball; % 势能杆重心高度 小球高度以支点为零势能点 V_rod M*g*(l/2)*cos(theta); V_ball m*g*x*cos(theta); V V_rod V_ball; L T - V;注意T_ball中2*x*Dx*Dtheta*cos(theta)是科氏力对应的动能交叉项漏掉此项将导致线性化后 A 矩阵缺失关键耦合元素后续极点配置失效。球杆系统建模分析.docx第 3.2 节明确指出此为高频错误点。2.2.3 符号求解运动微分方程% 拉格朗日方程对 theta 和 x 分别求导 eq1 diff(diff(L, Dtheta), t) - diff(L, theta) tau - b_k*Dtheta; % 杆方程含输入τ和阻尼 eq2 diff(diff(L, Dx), t) - diff(L, x) -b_x*Dx; % 小球方程无直接输入仅受摩擦 % 解出二阶导数 D2theta, D2x sol solve([eq1, eq2], [D2theta, D2x]); D2theta_eq simplify(sol.D2theta); D2x_eq simplify(sol.D2x);此时D2theta_eq和D2x_eq是含theta,x,Dtheta,Dx,tau的非线性表达式即系统原始动力学方程。2.2.4 线性化并生成状态空间模型% 在平衡点 [theta0, x0, Dtheta0, Dx0] 处线性化小角度近似 A_sym jacobian([Dtheta; Dx; D2theta_eq; D2x_eq], [theta; x; Dtheta; Dx]); B_sym jacobian([Dtheta; Dx; D2theta_eq; D2x_eq], tau); A0 double(subs(A_sym, {theta,x,Dtheta,Dx,tau}, {0,0,0,0,0})); B0 double(subs(B_sym, {theta,x,Dtheta,Dx,tau}, {0,0,0,0,0})); % 构建 ss 对象注意状态顺序[theta; x; Dtheta; Dx] sys_lin ss(A0, B0, [1 0 0 0], 0); % 输出为 theta杆角关键参数说明A0是 4×4 系统矩阵其特征值决定开环稳定性B0是 4×1 输入矩阵反映扭矩对各状态的影响权重。.rar中极点配置.m直接读取此sys_lin进行设计若此处线性化点选错如选 π 而非 0后续所有控制律将失效。3. 能控性与能观性验证不只是秩判据更是控制器部署的准入门槛3.1 能控性矩阵的物理意义与秩缺陷诊断球杆系统的能控性本质是能否通过单一输入扭矩 τ在有限时间内将系统从任意初始状态驱动到原点θ0, x0, ω0, v0数学上由能控性矩阵 $ \mathcal{C} [B\ AB\ A^2B\ A^3B] $ 的秩判定。但单纯rank(C)4不够——需检查条件数cond(C)。若cond(C) 1e6说明矩阵接近奇异数值计算中微小扰动会导致控制增益剧烈震荡实际硬件执行时电机易饱和。3.1.1 在 MATLAB 中执行完整能控性分析C_mat ctrb(A0, B0); fprintf(能控性矩阵秩: %d\n, rank(C_mat)); fprintf(能控性矩阵条件数: %.2e\n, cond(C_mat)); % 若条件数过大检查是否漏掉阻尼项 b_k 或 b_x if cond(C_mat) 1e5 warning(能控性矩阵病态请检查 b_k, b_x 是否为0或过小); end % 可视化能控性 Gramian李雅普诺夫方程解 Qc lyap(A0, -C_mat*C_mat); % 近似能控性 Gramian eig_Qc eig(Qc); fprintf(能控性 Gramian 特征值: [%.3f, %.3f, %.3f, %.3f]\n, eig_Qc);注意lyap(A0, -C_mat*C_mat)求解 $A_0 P P A_0^T -C C^T$其特征值反映各状态方向上的能控能量。若某特征值接近零如 1e-8对应状态几乎不可控——这在球杆系统中常表现为小球位置 x 的能控性远弱于杆角 θ需在控制器中降低 x 通道权重。3.2 能观性矩阵与传感器布置的强关联能观性回答仅测量杆角 θ输出 yθ能否唯一重构全部状态 [θ, x, ω, v]能观性矩阵 $ \mathcal{O} [C; CA; CA^2; CA^3]^T $ 的秩必须为 4。但现实中若只装一个编码器测 θC[1 0 0 0]则rank(O)通常为 3 —— 小球位置 x 成为不可观状态。.rar中判断系统能空能观.m的核心逻辑即在此。3.2.1 用 MATLAB 验证不同传感器配置的效果% 方案1仅测杆角 C1 [1 0 0 0]; O1 obsv(A0, C1); fprintf(仅测θ时能观性矩阵秩: %d\n, rank(O1)); % 通常为3 % 方案2增加小球位置传感器如直线电位器 C2 [1 1 0 0]; % 同时测θ和x O2 obsv(A0, C2); fprintf(测θx时能观性矩阵秩: %d\n, rank(O2)); % 应为4 % 方案3测θ和杆角速度ω编码器带速反馈 C3 [1 0 1 0]; O3 obsv(A0, C3); fprintf(测θω时能观性矩阵秩: %d\n, rank(O3)); % 验证是否足够提示球杆系统建模分析.docx第 4.1 节指出方案2虽能观但x传感器噪声大方案3更实用因ω可由θ微分获得需加低通滤波。.rar中未提供滤波代码实际部署时需在极点配置.m前插入y_filtered filter([1 0.9], [1 -0.9], y_raw);。4. 极点配置实现LQR控制器从理论公式到可部署的离散化代码4.1 为什么极点配置比PID更适合球杆系统PID 控制器在球杆系统中面临根本局限它本质是单输入单输出SISO设计而球杆系统是多输入多输出MIMO强耦合对象。当小球偏离杆中心时仅调节杆角无法快速抑制 x 振荡需同时协调 θ 与 x 的动态响应。极点配置通过状态反馈 $ u -Kx $直接将闭环极点置于期望位置实现多变量协同控制。MATLAB 的place()函数可精确配置 4 个极点但需满足能控性前提——这正是第 3 章验证的必要性。4.1.1 手动选择极点的工程准则极点类型期望位置物理意义球杆系统典型值主主导极点$-5 \pm 3j$决定整体响应速度与超调σ5 对应调节时间 ~0.8s快速衰减极点$-20$抑制高频振荡提升鲁棒性避免电机带宽限制零极点对消$-10$抵消系统右半平面零点若存在球杆系统通常无RHP零点注意.rar中极点配置.m默认使用p [-53j, -5-3j, -20, -10]但若你的杆长 l 增大需同比例减小实部如 l 加倍则 σ 减半否则控制器过度激进导致电机饱和。4.1.2 生成状态反馈增益 K 并验证闭环性能p [-53j, -5-3j, -20, -10]; % 期望极点 K place(A0, B0, p); % 计算反馈增益 fprintf(状态反馈增益 K [%.3f, %.3f, %.3f, %.3f]\n, K); % 构建闭环系统 A_cl A0 - B0*K; sys_cl ss(A_cl, B0, [1 0 0 0], 0); figure; step(sys_cl); title(闭环阶跃响应); % 检查闭环极点是否匹配 eig_cl eig(A_cl); fprintf(实际闭环极点: \n); disp(eig_cl);关键参数说明K是 1×4 行向量对应[k_theta, k_x, k_omega, k_v]。若eig_cl与p偏差 0.1说明place()数值不稳定应改用acker()或手动构造K (place(A0,B0,p))。4.2 从连续到离散嵌入式部署前的采样周期选择MATLAB 默认设计连续控制器但实际 DSP 或 STM32 需离散化。采样周期 $T_s$ 选择不当会导致性能恶化$T_s$ 过大10ms离散化引入相位滞后控制器响应迟钝$T_s$ 过小0.1msCPU 负载过高且 ADC 采样噪声放大。4.2.1 使用 c2d() 进行零阶保持离散化Ts 0.005; % 200Hz 采样率兼顾响应与负载 sys_d c2d(sys_cl, Ts, zoh); % 零阶保持离散化 % 提取离散状态空间矩阵 Ad sys_d.A; Bd sys_d.B; Cd sys_d.C; % 生成可用于 C 语言移植的差分方程 % x(k1) Ad*x(k) Bd*u(k) % y(k) Cd*x(k) fprintf(离散化后 Ad \n); disp(Ad); fprintf(离散化后 Bd \n); disp(Bd);提示.rar中未提供离散化代码但极点配置.m输出的K是连续域增益。若直接用于离散系统需同步离散化KKd K * (eye(4) - Ad)\Bd见《Digital Control of Dynamic Systems》第 4.3 节否则实际控制效果严重劣化。5. 稳定性边界验证用 Lyapunov 函数量化鲁棒裕度5.1 为什么特征值判据在非线性系统中不够用球杆系统的原始模型是非线性的含 sinθ, cosθ, x·ω² 项线性化仅在平衡点附近有效。当小球大幅滑动或杆角超过 ±15°线性控制器可能失稳。Lyapunov 稳定性理论提供全局/局部稳定域估计若存在正定函数 $V(x)0$ 且 $\dot{V}(x)0$则系统在该区域内渐近稳定。MATLAB 的lyap()可求解李雅普诺夫方程 $A^TP PA -Q$其中 $P$ 定义椭球形稳定域 $x^TPx c$。5.1.1 计算最大不变椭球稳定域Q eye(4); % 权重矩阵通常取单位阵 P lyap(A_cl, -Q); % 解 A_clP P*A_cl -Q % 计算稳定域半径 c需满足 x^TPx c 时闭环稳定 % 通过仿真找到最大 c 使得所有轨迹收敛 c_candidates logspace(-2, 1, 50); c_max 0; for c c_candidates % 初始化状态在椭球边界 x0 sqrt(c)*inv(chol(P))*randn(4,1) x0 sqrt(c) * inv(chol(P)) * randn(4,1); [~, y] ode45((t,x) (A_cl - B0*K)*x, [0 5], x0); if max(abs(y(:,1))) 0.1 max(abs(y(:,2))) 0.1 % θ和x均收敛 c_max c; end end fprintf(估计最大稳定域半径 c %.3f\n, c_max); % 可视化稳定域投影到θ-x平面 theta_grid linspace(-0.5, 0.5, 100); x_grid linspace(-0.3, 0.3, 100); [THETA, X] meshgrid(theta_grid, x_grid); V_vals zeros(size(THETA)); for i 1:length(theta_grid) for j 1:length(x_grid) x_vec [THETA(j,i); X(j,i); 0; 0]; % 初始速度为0 V_vals(j,i) x_vec * P * x_vec; end end contour(THETA, X, V_vals, [c_max c_max], LineWidth, 2, Color, r); title(Lyapunov 稳定域θ-x平面投影); xlabel(\theta (rad)); ylabel(x (m));注意此代码计算的是线性闭环系统的稳定域但实际非线性系统稳定域更小。.rar中未包含此验证而球杆系统建模分析.docx第 5.3 节强调若要求小球初始位置 |x₀| 0.15m则需将c_max设为 0.02 并重新设计 K。5.2 用 μ-分析评估参数摄动鲁棒性真实系统中杆长 l、小球质量 m 存在制造公差±5%阻尼系数 b_k 会随温度变化±30%。μ-分析结构奇异值量化系统在这些摄动下保持稳定的最大允许不确定性。MATLAB Robust Control Toolbox 提供musyn()和mussv()函数。5.2.1 构建摄动模型并计算鲁棒稳定裕度% 定义摄动块delta_l (l 变化), delta_m (m 变化), delta_bk (b_k 变化) delta_l ultidyn(delta_l, [1 1], Bound, 0.05); delta_m ultidyn(delta_m, [1 1], Bound, 0.05); delta_bk ultidyn(delta_bk, [1 1], Bound, 0.3); % 将摄动注入 A0 矩阵示例l 影响 J 和重力项 A_perturbed A0 delta_l*Adl delta_m*Adm delta_bk*Adbk; % Adl, Adm, Adbk 为各参数的灵敏度矩阵需从符号模型导出 % 构建不确定系统 sys_unc ss(A_perturbed, B0, [1 0 0 0], 0); % 计算结构奇异值 [mu_val, mu_frequencies] mussv(sys_unc, [], m); fprintf(最大结构奇异值 μ %.3f\n, max(mu_val)); if max(mu_val) 1 fprintf(系统对 ±5%% 参数摄动鲁棒稳定\n); else fprintf(系统在摄动下可能失稳请加强控制器鲁棒性\n); end提示.rar中未实现 μ-分析但球杆系统建模分析.docx第 6.2 节指出当 μ 0.8 时建议在 LQR 中增大 Q 矩阵中状态权重如Q diag([10, 5, 1, 1])以提升对参数变化的容忍度。本文还有配套的精品资源点击获取