高空作业车伸缩臂抖动抑制：微分平坦前馈与自抗扰控制

简介这是一份面向控制工程、机器人及农业装备方向研究者的学术论文资源聚焦伸缩臂在作业过程中的抖动抑制难题采用微分平坦理论与自抗扰控制ADRC相结合的思路展开研究。论文面向具备一定自动控制与动力学基础的中高级读者可解决臂体末端变形量不易测量、变臂长工况下控制器整定困难等实际问题。压缩包内仅含1个PDF文件大小约2.46MB为《农业机械学报》2020年第51卷第3期的正式排版论文含中英文摘要、正文、图表与参考文献便于直接阅读与引用。已有149人学习下载。文中将伸缩臂等效为扭簧连接的两刚性杆系统以变幅力矩为输入、两杆仰角为输出建立拉格朗日动力学模型并合成微分平坦输出构建单输入单输出二阶系统配合线性自抗扰控制器按杆长实时更新扩张状态观测器参数仿真与实验表明不同臂长下均能在2秒内消除抖动并保持仰角稳定对农机、机器人臂及医疗设备等场景的抑振控制设计具有直接参考价值。1. 伸缩臂抖动从哪里来为什么只靠 PID 压不住一台 40 米级的高空作业车臂架完全伸出后做一次 30° 的变幅动作末端往往要晃上三四秒才收得住。操作手最直接的感受是“车在摇”松手柄那一刻臂端还在余振。现场最常见的处理是把变幅回路的比例增益往下压、积分时间往上拉代价是动作变肉、末端定位靠人眼瞄。这种“调软换稳定”的妥协本质原因是抖动来自结构柔性而 PID 并不区分“刚体位置误差”和“柔性模态被激发”这两件性质完全不同的事——前者该用高增益压后者只能靠增益回避。微分平坦解决的是前一半问题把伸缩臂的状态和控制输入写成少数几个平坦输出及其导数的函数于是可以按解析式反推出走完一条给定轨迹所需的标称力矩顺手把容易激起振荡的高频成分从参考轨迹里剔掉。自抗扰控制解决后一半液压力滞、负载变化、未建模的二三阶模态、风载全部打包成一个“总扰动”交给扩张状态观测器估计并实时补偿。下面按“建模 → 平坦性验证 → 自抗扰整定 → 闭环仿真 → 实车取舍”的顺序把这条路线落成能跑的代码和能调的表。2. 用假设模态法建伸缩臂模型并验证微分平坦性伸缩臂是变截面箱型结构直接把偏微分方程塞进控制器不现实。工程上普遍的做法是假设模态法把臂架挠度展开成若干阶振型的线性组合只保留对抖动贡献最大的前一到两阶剩下的高阶模态当成噪声交给自抗扰去处理。这一步决定了后面所有参数的物理量纲值得先把模型压到最简。2.1 把连续伸缩臂压成 4 阶状态空间对根部铰接、末端自由的伸缩臂第 i 阶弯曲振型取悬臂梁形式φ_i(s) cosh(β_i s) − cos(β_i s) − σ_i [sinh(β_i s) − sin(β_i s)]其中 s ∈ [0,1] 是沿臂长的归一化坐标β_i 是频率方程 cos β · cosh β −1 的第 i 个根β₁1.8751β₂4.6941β₃7.8548σ_i (cosh β_i cos β_i)/(sinh β_i sin β_i)。用拉格朗日方程处理基座转角 θ 与一阶模态坐标 q₁ 的耦合得到两行方程J θ̈ c₁ q̈₁ u q̈₁ 2ζ₁ω₁ q̇₁ ω₁² q₁ −(c₁/m₁) θ̈其中 J ρA L³/3 是臂架绕根部的转动惯量m₁ ρA L ∫φ₁² ds 是一阶模态质量c₁ ρA L² ∫ s φ₁ ds 是模态耦合系数ω₁ β₁²√(EI/ρA)/L² 是一阶弯曲频率。把这些派生量对臂长的依赖列成表后面整定参数时全靠它。符号含义对臂长 L 的依赖J绕根部转动惯量L³m₁一阶模态质量∝ Lc₁模态耦合系数∝ L²ω₁一阶弯曲频率∝ L⁻²μ模态耦合因子 c₁²/(J m₁)近似与 L 无关b₀标称控制增益 1/[J(1−μ)]∝ L⁻³消去 q̈₁ 后以 x [θ, θ̇, q₁, q̇₁] 为状态、u 为输入得到一个四阶线性模型(1−μ) θ̈ u/J (c₁/J)(2ζ₁ω₁ q̇₁ ω₁² q₁) q̈₁ −(c₁/m₁) θ̈ − 2ζ₁ω₁ q̇₁ − ω₁² q₁μ c₁²/(J m₁) 是模态耦合因子对实际箱型伸缩臂一般小于 0.15所以 (1−μ) 接近 1但不能直接当 1 用——正是这一项把 b₀ 从 1/J 拉开。2.2 平坦输出选关节角还是臂端位移平坦输出的选择直接决定前馈逆是否稳定。如果平坦输出选得不好前馈里会出现非最小相位的零点逆运算就会把高频误差放大成发散的控制量。把上一节的模型写成传递函数取 y θ 时G(s) θ/u (1/J) · (s² 2ζ₁ω₁s ω₁²) / [ s² ( (1−μ)s² 2ζ₁ω₁s ω₁² ) ]分子是一阶模态的二阶多项式零点即 s² 2ζ₁ω₁s ω₁² 0 的根实部为 −ζ₁ω₁。只要 ζ₁ 0零点就在左半平面系统是最小相位的内部动态稳定。换成臂端位移输出 y θ γ q₁分子变成 (1 − γc₁/m₁)s² 2ζ₁ω₁s ω₁²一旦 γc₁/m₁ 1零点实部翻正前馈逆立刻不稳。这就解释了一个反直觉的现象想抑制臂端抖动平坦输出反而常常取关节角而不是臂端位移——臂端位移是控制目标关节角才是平坦输出。提示这一步一定要用数值真的算一遍零点实部别只看传递函数的形状。不同伸缩臂的 γ 差异很大凭直觉判断容易翻车。2.3 数值验证平坦性并生成前馈项下面这段脚本把建模、可控性检查、零点计算和前馈公式一次性串起来。import numpy as np from scipy.signal import ss2tf, tf2zpk def boom_model(L, EI, rhoA, zeta10.02): 单柔性模态伸缩臂返回 A、B 及派生物理量。 beta1 1.8751 s np.linspace(0, 1, 4001) sig (np.cosh(beta1) np.cos(beta1)) / (np.sinh(beta1) np.sin(beta1)) phi (np.cosh(beta1*s) - np.cos(beta1*s) - sig * (np.sinh(beta1*s) - np.sin(beta1*s))) J rhoA * L**3 / 3.0 # 绕根部转动惯量 m1 rhoA * L * np.trapz(phi**2, s) # 一阶模态质量 c1 rhoA * L**2 * np.trapz(s * phi, s) # 模态耦合系数 w1 beta1**2 * np.sqrt(EI / rhoA) / L**2 # 一阶弯曲频率 mu c1**2 / (J * m1) # 耦合因子 A np.zeros((4, 4)); B np.zeros((4, 1)) A[0, 1] 1.0 A[1, 2] c1 / (J * (1 - mu)) * w1**2 A[1, 3] c1 / (J * (1 - mu)) * 2 * zeta1 * w1 A[2, 3] 1.0 A[3, 2] -w1**2 / (1 - mu) A[3, 3] -2 * zeta1 * w1 / (1 - mu) B[1, 0] 1.0 / (J * (1 - mu)) B[3, 0] -c1 / (m1 * J * (1 - mu)) return A, B, dict(JJ, m1m1, c1c1, w1w1, mumu, zeta1zeta1, b01.0/(J*(1-mu))) def check_flatness(A, B): 可控性秩与传递函数零点实部判断平坦输出是否可用。 ctrb np.hstack([B, A B, A A B, A A A B]) rank np.linalg.matrix_rank(ctrb) C np.array([[1.0, 0.0, 0.0, 0.0]]) # 输出取关节角 theta num, den ss2tf(A, B, C, 0.0, input0) z, p, _ tf2zpk(num[0], den) return rank, z, p A, B, p boom_model(L16.0, EI2.0e7, rhoA120.0) rank, z, poles check_flatness(A, B) print(可控性秩 , rank, 状态维数 , A.shape[0]) print(零点实部 , np.round(z.real, 4)) print(极点实部 , np.round(poles.real, 4)) print(b0 标称值 , p[b0], 一阶频率 , p[w1], rad/s)以 L 16 m、EI 2.0e7 N·m²、ρA 120 kg/m 为例可控性矩阵满秩 4说明四个状态全部可控平坦输出存在零点实部约为 −0.113极点实部约为 −0.128两者都在左半平面前馈逆稳定。b₀ 的标称值落在 1e-5 量级这个数量级后面整定时要特别注意数值精度。boom_model里的EI和rhoA可以从臂架图纸估算也可以用一次实测的自由衰减振荡反推——测出臂端余振频率 ω₁ 和衰减包络按 ω₁ β₁²√(EI/ρA)/L² 反解 √(EI/ρA) 就行。check_flatness里输出矩阵C取关节角如果想试别的平坦输出只改C即可零点实部会跟着变。有了平坦性前馈力矩可以直接写出u_ff J(1−μ) θ̈_ref − c₁ (2ζ₁ω₁ q̇₁_ref ω₁² q₁_ref)其中 q₁_ref 由参考轨迹经模态方程滤波得到把 θ̈_ref 作为输入驱动 q̈₁ 2ζ₁ω₁q̇₁ ω₁²q₁ −(c₁/m₁)θ̈_ref再代入上式。3. 自抗扰控制器的关键参数ESO 带宽与 b0 调度3.1 把柔性耦合与摩擦滞回并进总扰动把第 2 章的模型改写成输入输出形式θ̈ b₀ u f其中 b₀ 1/[J(1−μ)]f 包含柔性模态反作用项 (c₁/[J(1−μ)])(2ζ₁ω₁q̇₁ ω₁²q₁)、未建模的高阶模态、液压阀的死区与滞回、以及风载和负载变化。把 f 扩张成第三个状态 x₃得到ẋ₁ x₂ ẋ₂ x₃ bu ẋ ḟ这里有个容易被忽略的取舍b₀ 越准z₃ 需要承担的就越少。如果直接把 b₀ 当成 1/J 用柔性模态那一项就会被划进总扰动ESO 的负担加倍。反过来 b₀ 偏小会让控制量整体偏大闭环等效增益偏移全靠 z₃ 去补观测器带宽不够时就会出现低频摆动。3.2 三阶 ESO 增益与控制器带宽的带宽参数化线性扩张状态观测器写成ż₁ z₂ − β₁(z₁ − θ) ż₂ z₃ − β₂(z₁ − θ) b₀u ż −β(z₁ − θ)把误差 e z₁ − θ 的特征多项式配成 (s ω_o)³直接得到三组增益β₁ 3ω_o, β₂ 3ω_o², β ω_o³控制律同样按带宽参数化u [ θ̈_ref k_p(θ_ref − z₁) k_d(θ̇_ref − z₂) − z₃ ] / b₀ k_p ω_c², k_d 2ω_c这里的 θ̈_ref 就是微分平坦给出的前馈项z₃ 负责把柔性耦合和扰动补上。两件事各管一段这也是这套组合比单独用其中任何一个更稳的原因。参数物理含义推荐取值主要约束ω_o观测器带宽(3~8)·ω₁且 ≤ 1/(3τ)编码器噪声、执行器时滞 τω_c控制器带宽ω_o/(3~5)与 ω₁ 至少留 2 倍裕度b₀控制增益标称值1/[J(1−μ)]随 L 在线更新用阶跃响应标定ζ₁模态阻尼比0.01~0.03钢结构中实测值三个参数里最需要拿捏的是 ω_o 和 ω_c 的比值这条经验规则来自 ESO 的估计精度——观测器比控制器快 3 到 5 倍时z₃ 对扰动的估计滞后已经不影响内环稳定性。继续拉大比值收益递减反而先把噪声放大。3.3 b0 随臂长 L 的调度表ω₁ ∝ L⁻²J ∝ L³所以臂架伸出的过程同时也是 b₀ 掉得飞快的过程。从 8 m 伸到 20 mb₀ 会掉一个数量级以上。下面这组数按 √(EI/ρA) 408 m²/s、ρA 120 kg/m 折算实际项目里用自己臂架的参数重算一遍。臂长 L (m)ω₁ (rad/s)b₀ 标称值建议 ω_c建议 ω_o822.44.9e-54.016.01210.01.5e-52.08.0165.66.1e-61.25.0203.63.1e-60.83.2调度方式有两种硬件资源紧张就按臂长分段查表段内取保守值控制周期宽裕就在线插值用实测臂长实时算 ω₁ β₁²√(EI/ρA)/L² 和 b₀ 1/[J(1−μ)]。伸缩动作本身是个慢过程插值带来的额外开销通常可以接受。4. 微分平坦前馈 自抗扰反馈的 Python 闭环实现4.1 参考轨迹五次多项式与 jerk 约束要让前馈真正起到抑制抖动的作用参考轨迹本身不能含高频。最省事的选择是五次多项式两端的速度、加速度、jerk 全为零θ_ref(t) θ₀ Δθ · (10τ³ − 15τ⁴ 6τ⁵), τ t/T对 τ 求导得到 θ̇_ref Δθ · (30τ² − 60τ³ 30τ⁴)/Tθ̈_ref Δθ · (60τ − 180τ² 120τ³)/T²。jerk 的峰值出现在 τ (3±√3)/6 处把动作时间 T 稍微拉长一点jerk 峰值下降的幅度比加速度峰值更快这对柔性臂特别划算。4.2 完整仿真回路import numpy as np class Boom: 单柔性模态伸缩臂被控对象。 def __init__(self, L16.0, EI2.0e7, rhoA120.0, zeta10.02): beta1 1.8751 s np.linspace(0, 1, 4001) sig (np.cosh(beta1)np.cos(beta1))/(np.sinh(beta1)np.sin(beta1)) phi (np.cosh(beta1*s) - np.cos(beta1*s) - sig*(np.sinh(beta1*s) - np.sin(beta1*s))) self.L, self.phi_tip L, 2.0 self.J rhoA*L**3/3.0 self.m1 rhoA*L*np.trapz(phi**2, s) self.c1 rhoA*L**2*np.trapz(s*phi, s) self.w1 beta1**2*np.sqrt(EI/rhoA)/L**2 self.zeta1 zeta1 self.mu self.c1**2/(self.J*self.m1) self.b0 1.0/(self.J*(1-self.mu)) def deriv(self, x, u): th, dth, q, dq x J, m1, c1, w1, z self.J, self.m1, self.c1, self.w1, self.zeta1 d2th (u/J (c1/J)*(2*z*w1*dq w1**2*q)) / (1 - self.mu) d2q -(c1/m1)*d2th - 2*z*w1*dq - w1**2*q return np.array([dth, d2th, dq, d2q]) def tip_angle(self, x): 臂端等效角度 刚体转角 一阶模态在末端的贡献。 return x[0] self.phi_tip * x[2] / self.L class LADRC: 三阶线性自抗扰控制器带宽参数化。 def __init__(self, b0, wc, wo, dt): self.b0, self.dt b0, dt self.kp, self.kd wc**2, 2.0*wc self.b1, self.b2, self.b3 3*wo, 3*wo**2, wo**3 self.z np.zeros(3) # z1 估计 thetaz2 估计 dthetaz3 估计总扰动 self.u_prev 0.0 def reset(self, y, dy0.0): self.z[:] [y, dy, 0.0] self.u_prev 0.0 def step(self, y, ref, dref, ddref): e self.z[0] - y # 欧拉离散z2 更新用上一拍控制量避免代数环 self.z[0] self.dt * (self.z[1] - self.b1*e) self.z[1] self.dt * (self.z[2] - self.b2*e self.b0*self.u_prev) self.z[2] self.dt * (-self.b3*e) u0 ddref self.kp*(ref - self.z[0]) self.kd*(dref - self.z[1]) u (u0 - self.z[2]) / self.b0 self.u_prev u return u def s_curve(t, T, th0, th1): 五次多项式参考轨迹返回 (角度, 角速度, 角加速度)。 if t 0.0: return th0, 0.0, 0.0 if t T: return th1, 0.0, 0.0 tau t / T s 10*tau**3 - 15*tau**4 6*tau**5 ds 30*tau**2 - 60*tau**3 30*tau**4 dds 60*tau - 180*tau**2 120*tau**3 dth th1 - th0 return th0 dth*s, dth*ds/T, dth*dds/T**2 def run(use_adrcTrue, T2.0, dt0.001, U_MAX4.0e4): boom Boom(L16.0) ctrl LADRC(b0boom.b0, wc1.2, wo5.0, dtdt) x np.zeros(4); ctrl.reset(x[0]) ts, th, tip, us [], [], [], [] for k in range(int(T/dt) 1): t k * dt r, dr, ddr s_curve(t, T, 0.0, np.deg2rad(30.0)) if use_adrc: u ctrl.step(x[0], r, dr, ddr) else: # 对照组同样的带宽参数但不做 ESO 扰动补偿 u (ddr ctrl.kp*(r - x[0]) ctrl.kd*(dr - x[1])) / boom.b0 u float(np.clip(u, -U_MAX, U_MAX)) k1 boom.deriv(x, u); k2 boom.deriv(x 0.5*dt*k1, u) k3 boom.deriv(x 0.5*dt*k2, u); k4 boom.deriv(x dt*k3, u) x x dt/6.0*(k1 2*k2 2*k3 k4) ts.append(t); th.append(x[0]); tip.append(boom.tip_angle(x)); us.append(u) return boom, np.array(ts), np.array(th), np.array(tip), np.array(us) def metrics(t, tip, T, band5e-4): 残余振幅动作结束后的最大偏离与稳定时间。 mask t T resid np.max(np.abs(tip[mask] - tip[mask][-1])) idx np.where(np.abs(tip - tip[-1]) band)[0] settle t[idx[-1]] if len(idx) else 0.0 return resid, settle boom, t, th, tip, u run(use_adrcTrue) r_adrc, s_adrc metrics(t, tip, T2.0) boom2, t2, th2, tip2, u2 run(use_adrcFalse) r_pd, s_pd metrics(t2, tip2, T2.0) print(f带 ESO 补偿 残余振幅{r_adrc:.3e} rad 稳定时间{s_adrc:.3f} s) print(f纯 PD 前馈 残余振幅{r_pd:.3e} rad 稳定时间{s_pd:.3f} s)代码里有几个点值得单独说明。LADRC.step中 z₂ 的更新用的是self.u_prev而不是当前拍刚算出的 u这是为了避免同一拍里 u 既依赖 z₃ 又参与 z₂ 更新形成代数环dt 1 ms 下这个滞后完全可以忽略但如果控制周期拉长到 20 ms 以上就应改用双线性离散或直接推离散化增益。np.clip那行是必须的ESO 在启动瞬间 z₃ 尚未收敛未经限幅的控制量会打出远超液压系统能力的指令。Boom.tip_angle把一阶模态贡献按 φ₁(1)·q₁/L 折算成等效角度φ₁(1) ≈ 2.0 是归一化振型在自由端的取值。这个换算关系决定了后面残余振幅的量纲——直接拿它和角度传感器读数比较就能判断模型保真度够不够。4.3 抖动抑制效果的量化对比跑完之后重点看两组数。第一组是残余振幅动作结束后臂端角度围绕终值的最大偏离。如果这个值在 ω₁ 附近有明显频率成分说明总扰动中还有 ESO 没跟上的部分把 ω_o 往 5ω₁ 以上抬一档再试。第二组是稳定时间臂端进入 ±0.5 mrad 带后不再出去的时刻。稳定时间降不下来但残余振幅已经很小问题通常不在控制器而在参考轨迹本身的 jerk 还不够小把动作时间 T 拉长 30% 再看。纯 PD 对照组的意义在于给出下界——如果两者差距不到 2 倍多半是 ω_c 设得太低或者 b₀ 用错了量纲。5. 实车落地观测器带宽、噪声与时滞的取舍5.1 ωo 提升到哪一步就该停ESO 的代价是把测量噪声一阶差分放大z₃ 里会叠进编码器量化噪声和倾角传感器的抖动再通过 u (u₀ − z₃)/b₀ 直接打进阀控电流。经验规律是 ω_o 每提高一倍控制量中的噪声幅值大约也翻一倍。有两条硬线不能碰一是 1/(3τ)τ 是实测的阀芯响应加管路滞后二是编码器分辨率换算出的等效角速度噪声水平一旦噪声幅度接近 z₃ 的真实扰动幅度继续加带宽就只剩坏处。实操上从 3ω₁ 起调每次加 20%盯控制电流的纹波和臂端高频颤振出现两者之一就退回来。5.2 电液执行器时滞的补偿比例阀的响应、管路容腔和油液可压缩性等效成一阶或纯延时。如果实测 τ 超过控制周期的 3 倍单纯提高 ω_o 已经无效此时按 ω_o ≤ 1/(3τ) 把带宽封顶再用 Smith 预估器把时滞从观测器回路里摘出去预估器输出的是无延时位置用这个位置去驱动 ESO实际采回的 θ 只用于修正预估器偏差。5.3 b0 分段整定与现场验证b₀ 最可靠的标定方法是在每个臂长段做一次小角度阶跃记录稳态角加速度 θ̈ 和控制量增量 Δub₀ θ̈/Δu。注意剔除阶跃开始的前 50 ms那一段还在克服阀的死区。把各段的 b₀ 和 ω₁ 列成调度表烧进控制器中间臂长按线性插值取。几个典型的误调特征可以直接对照判断控制量低频摆动是 ω_o 太低z₃ 估计跟不上柔性模态臂端出现高频颤振是 ω_o 太高噪声进了控制量动作结束瞬间有一段持续等幅振荡通常是 b₀ 偏小导致闭环实际增益偏高把 b₀ 乘以 0.8 再试。这套组合真正的价值不在某个单一参数有多高而在于前馈负责把已知的、可解析的部分一次性做对自抗扰只处理剩下的那一小部分两者的分工越清楚需要现场拿捏的参数就越少。本文还有配套的精品资源点击获取