受迫摆仿真中超调/振荡的原因及PID角度控制优化方案
我编写了一个离散时间的摆仿真程序，无外力作用时运行正常。但当尝试应用PI(D)控制使其移动并到达特定角度时，尽管调整增益后能到达目标角度，却无法稳定在该角度，无论如何调整参数都会出现异常。除了反复调整增益外，是否有更规范的PID控制方法来实现角度控制？
import numpy as np import pandas as pd import matplotlib.pyplot as plt m = 1.0 # mass (kg) k = 0.0 # spring constant k_{s} r = 1.0 # link length (m) g = 9.81 # gravity term (m/s^{2}) damp = 1.0 # damping term (for unforced pendulum) start_time = 0.0 # start time end_time = 50.0 # end time dt = 0.01 # time increment num_steps = int((end_time - start_time) / dt) t = np.linspace(start_time, end_time, num_steps) theta = np.zeros_like(t) thetadot = np.zeros_like(t) thetaddot = np.zeros_like(t) tau_u = np.zeros_like(t) error = np.zeros_like(t) # PID control gains kp = 1.2 ki = 0.0 kd = 0.0 # ICs theta[0] = -np.pi/2 thetadot[0] = 0.0 thetaddot[0] = 0.0 tau_u[0] = 0.0 target_angle = 0.0 error[0] = theta[0] - target_angle for i in range(1, num_steps): # mod term to keep angle between 0 and 2pi radians if theta[i] > 2*np.pi: theta[i] = theta[i] % 2*np.pi error[i] = theta[i - 1] - target_angle #target_angle - theta[i - 1] tau_u[i] = (kp * error[i - 1]) - (ki * np.sum(error) * dt) - (kd * thetadot[i - 1]) thetadot[i] = thetadot[i - 1] - ( ((3/(m*r**2))*tau_u[i]) + (((3*k)/(m*r**2))*theta[i - 1]) + ((3*g*np.cos(theta[i - 1]))/(2*r)) )* dt theta[i] = (theta[i - 1] + (thetadot[i] * dt) ) # plotting omitted for space
一、代码中的核心错误
你的仿真无法稳定的根本原因不是PID参数问题，而是动力学模型错误和PID实现逻辑错误，先修正这些基础问题，再谈参数整定：
动力学方程错误
单摆的重力矩由$\sin\theta$决定（而非$\cos\theta$），正确的转动动力学方程（质点摆，转动惯量$I=mr^2$）应为：
$$\ddot{\theta} = -\frac{g}{r}\sin\theta - \frac{damp}{I}\dot{\theta} + \frac{\tau_u}{I}$$
你代码中错误地使用了$\cos\theta$，导致重力项的方向完全错误，这是无法稳定的核心原因。另外，你定义了阻尼项damp但未加入方程，也会导致系统无阻尼振荡。角度误差的周期性处理错误
摆的角度是周期性的（周期$2\pi$），误差计算时需将角度差限制在$[-\pi, \pi]$范围内，否则会出现“绕圈”误差（比如从$-\pi/2$到0的误差是$\pi/2$，但如果错误计算成$3\pi/2$，PID会输出反向控制）。你当前的mod逻辑位置错误（theta[i]还未赋值就判断），完全起不到作用。PID误差符号与实现错误
- 误差应定义为目标角度 - 当前角度，你写反了，会导致正反馈，加剧振荡。
- 积分项直接用
np.sum(error)效率极低且容易累积错误，应维护一个单独的累积积分变量。 - 微分项的符号需与控制逻辑匹配，避免反向阻尼。
- 离散积分/微分的实现问题
你当前的欧拉积分顺序有问题，应先计算当前步的误差，再计算控制量，最后更新状态。
二、修正后的仿真代码
以下是修正后的代码，解决了上述所有问题：
import numpy as np import matplotlib.pyplot as plt m = 1.0 # mass (kg) r = 1.0 # link length (m) g = 9.81 # gravity term (m/s²) damp = 0.5 # damping coefficient I = m * r**2 # moment of inertia start_time = 0.0 end_time = 50.0 dt = 0.01 num_steps = int((end_time - start_time) / dt) t = np.linspace(start_time, end_time, num_steps) theta = np.zeros_like(t) thetadot = np.zeros_like(t) tau_u = np.zeros_like(t) error = np.zeros_like(t) # PID control gains kp = 15.0 ki = 0.5 kd = 2.0 # Initial conditions theta[0] = -np.pi / 2 thetadot[0] = 0.0 target_angle = 0.0 integral = 0.0 # 单独维护积分累积量 for i in range(1, num_steps): # 1. 计算当前误差并处理周期性 current_theta = theta[i-1] raw_error = target_angle - current_theta # 将误差限制在[-π, π] error[i] = np.arctan2(np.sin(raw_error), np.cos(raw_error)) # 2. 计算PID控制量 integral += error[i] * dt derivative = -thetadot[i-1] # 用角速度代替误差微分，符号匹配 tau_u[i] = kp * error[i] + ki * integral + kd * derivative # 3. 更新动力学状态（欧拉法） # 正确的摆动力学方程：θ'' = (τ - damp*θ' - mgr sinθ)/I thetaddot = (tau_u[i] - damp * thetadot[i-1] - m*g*r*np.sin(current_theta)) / I thetadot[i] = thetadot[i-1] + thetaddot * dt theta[i] = theta[i-1] + thetadot[i] * dt # 绘图 plt.figure(figsize=(12,6)) plt.subplot(211) plt.plot(t, theta, label='摆角') plt.axhline(target_angle, color='r', linestyle='--', label='目标角度') plt.legend() plt.subplot(212) plt.plot(t, tau_u, label='控制力矩') plt.legend() plt.show()
三、规范的PID参数整定方法
修正模型后，若仍需优化参数，可采用以下规范方法，而非盲目试凑：
基于系统模型的整定
对于线性化后的摆系统（小角度下$\sin\theta\approx\theta$），系统传递函数为：
$$G(s) = \frac{1}{I s^2 + damp s + mgr}$$
可通过极点配置法直接计算PID参数，使闭环系统的极点位于期望的位置（比如欠阻尼响应，阻尼比$\zeta=0.707$）。Ziegler-Nichols整定法
- 先设$K_i=0, K_d=0$，逐渐增大$K_p$直到系统出现持续等幅振荡，记录此时的临界增益$K_{cr}$和振荡周期$T_{cr}$。
- 按公式计算PID参数：
- PI控制：$K_p=0.45K_{cr}, K_i=0.54K_{cr}/T_{cr}$
- PID控制：$K_p=0.6K_{cr}, K_i=2K_p/T_{cr}, K_d=K_p T_{cr}/8$
- 手动整定步骤（经验法）
- 先调比例项$K_p$：从0开始增大，直到系统能快速趋近目标但略有振荡。
- 再调积分项$K_i$：从0开始增大，消除稳态误差，同时避免振荡加剧。
- 最后调微分项$K_d$：从0开始增大，抑制振荡，提高系统稳定性。
内容的提问来源于stack exchange，提问作者Charlie
