【自动控制入门2B】从零搭建全连续控制系统：基于抗积分饱和PID的输入限制直线运动物体位置控制仿真程序

引言

本文分享【自动控制入门2A】从零搭建全连续控制系统：基于抗积分饱和PID的输入限制直线运动物体位置控制文章里仿真使用的程序，包括从零搭建的python仿真程序和用于对比的matlab/simulink仿真模型，其中python程序仅用到常用数据处理库和绘图库，控制系统的核心逻辑均为从零实现。编程不难，贵在逻辑梳理。

1 python 仿真程序

程序主要包括被控对象类、控制器类、传感器类、参数设置、仿真循环、控制性能计算、仿真数据导出、仿真结果绘图等几个模块，均有详细注释。

程序中设有打开输入限幅、打开抗积分饱和、仿真结果绘图、实验记录导出和实验数据导出等开关，可以根据需要打开。实验记录导出和实验数据导出效果见1.1和1.2展示。

被控对象类的微分方程函数ode_plant(self,state0=[0,0,0,0])包含分别考虑空气阻力和重力的四种情况，可以根据需要取消注释对应微分方程的代码，共可以观察四种系统的特性（包括线性系统和非线性系统）。

1.1 实验数据导出效果

1.2 实验记录导出效果

1.3 完整程序

本文使用python=3.10环境，需要提前安装如下第三方库
pip install numpy、pandas、matplotlib、openpyxl
之后如下完整代码可直接运行

frommatplotlibimportpyplotaspltimportnumpyasnpimporttimeimportpandasaspdimportos# 2025.12.4# 针对考虑空气阻力的竖直方向运动质点的位置控制问题，设计了一种条件积分的PID控制器# 使用RK4方法进行全连续系统仿真classPlant:"""被控对象：考虑空气阻力的竖直方向运动质点模型"""def__init__(self,state0=[0,0],params=None,dt=0.001)->None:self.m=params['m']# 质量self.g=params['g']# 重力加速度self.k=params['k']# 空气阻力系数self.dt=dt# 时间步长# 初始状态self.h=state0[0]# 初始竖直方向位置self.v=state0[1]# 初始竖直方向速度# TODO: 调整传感器噪声水平self.sensor_h=Sensor(noise_level=params['noise_level'])# TODO: 调整 PID 参数 和 控制器kp=params['kp']ki=params['ki']kd=params['kd']amplitude_limit=params['amplitude_limit']# 控制输入限幅开关，False 表示不启用限幅rate_limit=params['rate_limit']# 控制输入变化率限幅开关，False 表self.controller_h2u=ControllerPID(kp=kp,ki=ki,kd=kd,amplitude_limit=amplitude_limit,rate_limit=rate_limit,dt=dt)defode_plant(self,state0=[0,0,0,0]):""" 常微分方程组 (ODE)：质点的动力学模型（h, v）+ 控制器状态（i_e, e） F_net = F_control - F_gravity + F_air_resistance dv/dt = (u - m*g + f_air) / m """# 状态变量初值h,v=state0[0],state0[1]i_e,e=state0[2],state0[3]# 控制器微分方程组de=-v di_e,u=self.controller_h2u.control(e,i_e,de)# 空气阻力计算air_resist=-self.k*v*abs(v)# 质点运动微分方程组dh=v# 考虑重力和空气阻力dv=(u-self.m*self.g+air_resist)/self.m# 不考虑重力，只考虑空气阻力# dv = (u + air_resist) / self.m# 不考虑空气阻力，只考虑重力# dv = (u - self.g * self.m) / self.m# 重力和空气阻力都不考虑# dv = u / self.mdstate=np.array([dh,dv,di_e,de],dtype=float)returndstate.reshape((len(dstate),1))defupdate(self):"""使用 RK4 方法更新系统状态"""self.controller_h2u.e=self.controller_h2u.state_desired-self.sensor_h.measure(self.h)state=np.array([self.h,self.v,self.controller_h2u.i_e,self.controller_h2u.e],dtype=float)state=state.reshape((len(state),1))k1=self.ode_plant(list(state.flatten()))state1=state+k1*self.dt/2k2=self.ode_plant(list(state1.flatten()))state2=state+k2*self.dt/2k3=self.ode_plant(list(state2.flatten()))state3=state+k3*self.dt k4=self.ode_plant(list(state3.flatten()))d_state=(k1+2*k2+2*k3+k4)/6state+=d_state*self.dt# state 和 d_state 是 (4,1) 的 numpy 数组，直接取标量 floatself.h=float(state[0,0])self.v=float(state[1,0])self.controller_h2u.i_e=float(state[2,0])self.controller_h2u.e=float(state[3,0])returnstateclassControllerPID:"""PID 控制器，采用条件积分法实现抗积分饱和 (Anti-Windup)"""def__init__(self,kp=6,ki=0,kd=0,amplitude_limit=False,rate_limit=False,dt=0.1)->None:self.kp=kp self.ki=ki self.kd=kd self.dt=dt self.state_desired=0self.i_e=0self.e=0self.amplitude_limit=amplitude_limit self.rate_limit=rate_limit self.u_raw=0self.u=0defamplitude_restriction(self,u):"""控制信号幅值限制"""upper_bound=self.amplitude_limit lower_bound=-self.amplitude_limit flag_restrict=