简介：这是一份面向控制工程、机器人及农业装备方向研究者的学术论文资源，聚焦伸缩臂在作业过程中的抖动抑制难题，采用微分平坦理论与自抗扰控制（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₁ | 一阶模态质量 | ∝ L |
| c₁ | 模态耦合系数 | ∝ 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, zeta1=0.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(J=J, m1=m1, c1=c1, w1=w1, mu=mu, zeta1=zeta1, b0=1.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, input=0) z, p, _ = tf2zpk(num[0], den) return rank, z, p A, B, p = boom_model(L=16.0, EI=2.0e7, rhoA=120.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 m，b₀ 会掉一个数量级以上。下面这组数按 √(EI/ρA) = 408 m²/s、ρA = 120 kg/m 折算，实际项目里用自己臂架的参数重算一遍。
| 臂长 L (m) | ω₁ (rad/s) | b₀ 标称值 | 建议 ω_c | 建议 ω_o |
|---|---|---|---|---|
| 8 | 22.4 | 4.9e-5 | 4.0 | 16.0 |
| 12 | 10.0 | 1.5e-5 | 2.0 | 8.0 |
| 16 | 5.6 | 6.1e-6 | 1.2 | 5.0 |
| 20 | 3.6 | 3.1e-6 | 0.8 | 3.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, L=16.0, EI=2.0e7, rhoA=120.0, zeta1=0.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 估计 theta，z2 估计 dtheta，z3 估计总扰动 self.u_prev = 0.0 def reset(self, y, dy=0.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_adrc=True, T=2.0, dt=0.001, U_MAX=4.0e4): boom = Boom(L=16.0) ctrl = LADRC(b0=boom.b0, wc=1.2, wo=5.0, dt=dt) 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, band=5e-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_adrc=True) r_adrc, s_adrc = metrics(t, tip, T=2.0) boom2, t2, th2, tip2, u2 = run(use_adrc=False) r_pd, s_pd = metrics(t2, tip2, T=2.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₀ 最可靠的标定方法是在每个臂长段做一次小角度阶跃：记录稳态角加速度 θ̈ 和控制量增量 Δu，b₀ = θ̈/Δu。注意剔除阶跃开始的前 50 ms，那一段还在克服阀的死区。把各段的 b₀ 和 ω₁ 列成调度表烧进控制器，中间臂长按线性插值取。
几个典型的误调特征可以直接对照判断：控制量低频摆动是 ω_o 太低，z₃ 估计跟不上柔性模态；臂端出现高频颤振是 ω_o 太高，噪声进了控制量；动作结束瞬间有一段持续等幅振荡，通常是 b₀ 偏小导致闭环实际增益偏高，把 b₀ 乘以 0.8 再试。这套组合真正的价值不在某个单一参数有多高，而在于前馈负责把已知的、可解析的部分一次性做对，自抗扰只处理剩下的那一小部分，两者的分工越清楚，需要现场拿捏的参数就越少。
本文还有配套的精品资源，点击获取
