风机溶解氧串级控制仿真 —— 主环DO/副环风压双闭环调节“那年污水曝气池DO老是忽高忽低单PID怎么调都稳不住。后来改成串级控制外环盯DO内环跟风压风机频率跟着风压走DO立马稳在±0.2mg/L以内。那一刻我才明白单闭环是‘头痛医头’串级才是‘标本兼治’。”—— 哈尔滨工程大学《工业过程控制》课程核心思想延伸一、实际应用场景描述在城市污水处理厂、水产养殖、发酵工程等需氧生物处理场景中溶解氧DO浓度是决定微生物活性的关键参数┌──────────────────────────────────────────────┐│ 风机溶解氧串级控制系统 ││ ││ [生物反应池] ││ • 活性污泥微生物 ││ • 需氧量随水质/水量波动 ││ • DO目标: 2.0±0.3 mg/L ││ │ ││ ▼ 空气扩散 ││ ┌────────────────────────────┐ ││ │ 曝气系统 (Diffusers) │ ││ │ • 微孔曝气头 │ ││ │ • 氧转移效率受水温/堵塞影响 │ ││ │ • 大滞后 (30~60s) │ ││ └────────────┬───────────────┘ ││ │ 空气流量 ││ ▼ ││ ┌────────────────────────────┐ ││ │ 罗茨/离心风机 │ ││ │ • 变频调速 (20~50Hz) │ ││ │ • 风压-流量非线性特性 │ ││ │ • 响应较快 (5~10s) │ ││ └────────────┬───────────────┘ ││ │ 风压反馈 ││ ▼ ││ ┌────────────────────────────┐ ││ │ 副PID控制器 (Inner Loop) │ ││ │ • 设定值: 主环输出 │ ││ │ • 测量值: 风压变送器 │ ││ │ • 输出: 风机频率 │ ││ │ • 快速响应抑制内扰 │ ││ └────────────┬───────────────┘ ││ │ 频率指令 ││ ▼ ││ ┌────────────────────────────┐ ││ │ 主PID控制器 (Outer Loop) │ ││ │ • 设定值: DO目标值 │ ││ │ • 测量值: DO分析仪 │ ││ │ • 输出: 风压设定值 │ ││ │ • 慢响应克服外扰 │ ││ └────────────┬───────────────┘ ││ │ DO反馈 ││ ▼ ││ ┌────────────────────────────┐ ││ │ DO在线分析仪 │ ││ │ • 荧光法/电化学法 │ ││ │ • 测量滞后 (10~30s) │ ││ │ • 易受污泥干扰 │ ││ └───────────────────────────┘ ││ ││ 核心: 主环稳DO 副环抗风压波动 串级解耦 节能优化 │└──────────────────────────────────────────────┘单PID vs 串级控制维度 单PID控制 串级控制抗扰能力 ❌ 风压波动直接影响DO ✅ 副环快速抑制风压波动响应速度 ❌ 滞后大超调明显 ✅ 副环快速响应主环慢稳稳态精度 ❌ 受非线性影响大 ✅ 主环消除余差参数整定 ❌ 互相牵制难整定 ✅ 先副后主分步整定节能效果 ❌ 风机常开大能耗高 ✅ 按需供氧节能15~30%二、引入痛点2.1 现场的真实困境场景 现场发生了什么 根因“DO 过山车” “DO 在 1.0~4.0mg/L 之间大幅震荡” 单PID无法克服风压波动“污泥发黑” “好氧池局部缺氧污泥发黑发臭” DO控制不稳局部死区“能耗爆表” “风机常年50Hz满转电费惊人” 无按需供氧过量曝气“出水超标” “COD/BOD 出水超标环保罚款” DO控制精度不足“风机损坏” “风机频繁启停轴承烧坏” 控制振荡导致设备疲劳2.2 核心矛盾生物处理的核心是“稳氧”而曝气的难点是“风压波动氧转移滞后”。单PID试图用一个控制器同时解决“快速抗扰”和“慢速稳氧”本质上是“一脑两用”的矛盾。串级控制通过主环管“稳”、副环管“快”**的分工实现了控制性能的质的飞跃。2.3 我们要解决什么用一段精简的 Python 程序构建一个 风机溶解氧串级控制仿真系统实现1. 双闭环结构 —— 主环DO、副环风压2. 串级耦合 —— 主环输出作为副环设定值3. 非线性风机 —— 频率-风压-流量特性4. 氧转移滞后 —— 曝气-溶解-传质全过程5. 抗积分饱和 —— 主副环双重保护6. 性能量化 —— ISE/IAE/TV 评价指标三、核心逻辑讲解3.1 理论基础串级控制系统本工具基于哈工程《工业过程控制》第六章“复杂控制系统”和第四章“PID控制器”① 串级控制结构┌──────────┐ ┌──────────┐ ┌──────────┐│ 主给定 │─────▶│ 主PID │─────▶│ 副给定 ││ DO_sp │ │ (Outer) │ │ P_sp │└──────────┘ └────┬─────┘ └────┬─────┘│ │┌──────────┐ ┌────▼─────┐ ┌────▼─────┐│ 副反馈 │◀─────│ 副PID │◀─────│ 执行器 ││ P_pv │ │ (Inner) │ │ 风机频率 │└──────────┘ └────┬─────┘ └────┬─────┘│ │┌──────────┐ ┌────▼─────┐ ┌────▼─────┐│ 主反馈 │◀─────│ 被控对象 │◀─────│ 风压/流量││ DO_pv │ │ (曝气池) │ │ (非线性) │└──────────┘ └──────────┘ └──────────┘② 主副环分工环路 任务 响应速度 参数特点主环DO 克服水质/水量外扰保证稳态精度 慢τ30~60s Kp小、Ti大、Td适中副环风压 克服风机内扰快速跟踪设定 快τ5~10s Kp大、Ti小、Td小③ 抗积分饱和策略主环抗饱和:if (主PID输出 副环设定上限) AND (主误差 0):停止主积分累积if (主PID输出 副环设定下限) AND (主误差 0):停止主积分累积副环抗饱和:if (风机频率 50Hz) AND (副误差 0):停止副积分累积if (风机频率 20Hz) AND (副误差 0):停止副积分累积④ 串级控制优势- 抑制内扰副环快速抑制风机风压波动- 改善动态副环加速响应减小主环滞后影响- 适应非线性副环可补偿风机部分非线性- 提高鲁棒性主环参数对副环变化不敏感3.2 被控对象建模① 风机模型非线性P K_p \cdot f^2 K_f \cdot Q^2其中- P 风压kPa- f 频率Hz- Q 流量隐含- K_p, K_f 风机特性系数② 氧转移模型大滞后\frac{dC}{dt} K_La \cdot (C^* - C) - OUR其中- C 溶解氧浓度mg/L- C^* 饱和溶解氧随温度/盐度变化- K_La 氧转移系数与风压正相关- OUR 微生物耗氧速率随水质波动四、代码讲解面向对象设计4.1 类结构总览类名 职责 设计模式DOControlConfig DO控制配置dataclass 值对象BlowerConfig 风机配置dataclass 值对象PIDController 抗饱和PID控制器 模板方法BlowerModel 风机非线性模型 领域模型AerationTankModel 曝气池氧转移模型 领域模型CascadeController 串级控制器聚合根 聚合根DisturbanceGenerator 扰动发生器 工厂模式PerformanceAnalyzer 性能分析器 封装Visualizer 可视化工具 封装4.2 核心代码完整可运行完整源码约 350 行包含 8 个类、串级PID、非线性风机、氧转移滞后、性能分析。以下为精简核心版可直接复制运行。detailssummary 完整源码点击展开/折叠/summary风机溶解氧串级控制仿真 —— 主环DO/副环风压双闭环参考哈尔滨工程大学《工业过程控制》第六章复杂控制系统from dataclasses import dataclass, fieldfrom typing import List, Tuple, Optional, Dict, Dequefrom enum import Enum, autoimport numpy as npimport matplotlib.pyplot as pltfrom collections import dequeimport mathfrom datetime import datetime# # 1. 基础数据结构值对象# dataclassclass DOControlConfig:DO串级控制配置 —— 值对象# 主环DO环参数do_sp: float 2.0 # DO设定值 (mg/L)do_kp: float 0.8 # 主环比例增益do_ki: float 0.02 # 主环积分增益do_kd: float 0.1 # 主环微分增益do_output_min: float 10.0 # 主环输出下限 (kPa)do_output_max: float 50.0 # 主环输出上限 (kPa)# 副环风压环参数press_kp: float 1.5 # 副环比例增益press_ki: float 0.08 # 副环积分增益press_kd: float 0.05 # 副环微分增益press_output_min: float 20.0 # 副环输出下限 (Hz)press_output_max: float 50.0 # 副环输出上限 (Hz)# 抗积分饱和anti_windup: bool True# 过程参数sampling_time: float 1.0 # 采样周期 (s)max_sim_time: float 3600.0 # 最大仿真时间 (s)dataclassclass BlowerConfig:风机配置 —— 值对象rated_power: float 55.0 # 额定功率 (kW)max_freq: float 50.0 # 最大频率 (Hz)min_freq: float 20.0 # 最小频率 (Hz)pressure_coeff: float 0.02 # 风压-频率系数 (P ∝ f²)flow_coeff: float 0.05 # 流量-频率系数response_time: float 3.0 # 风机响应时间常数 (s)efficiency: float 0.75 # 风机效率dataclassclass AerationConfig:曝气池配置 —— 值对象volume: float 5000.0 # 池容 (m³)depth: float 5.0 # 水深 (m)sat_do: float 9.0 # 饱和溶解氧 (mg/L, 20℃)kl_a_base: float 0.15 # 氧转移系数基准 (1/min)temp_coeff: float 1.024 # 温度系数 (θ)alpha: float 0.85 # 污水修正系数beta: float 0.95 # 盐度修正系数our_base: float 0.05 # 基准耗氧速率 (mg/L/min)# # 2. PID控制器模板方法# class PIDController:抗积分饱和PID控制器 —— 模板方法def __init__(self, kp: float, ki: float, kd: float,output_min: float, output_max: float,anti_windup: bool True, sample_time: float 1.0):self.kp kpself.ki kiself.kd kdself.output_min output_minself.output_max output_maxself.anti_windup anti_windupself.ts sample_time# 状态变量self.integral 0.0self.prev_error 0.0self.prev_output 0.0self.prev_pv 0.0self._last_reset_time 0.0def compute(self, setpoint: float, process_variable: float,external_limit: Optional[Tuple[float, float]] None) - float:计算PID输出error setpoint - process_variable# 比例项p_term self.kp * error# 积分项带抗饱和if self.anti_windup:# 检查是否被外部限制if external_limit is not None:ext_min, ext_max external_limitif (self.prev_output ext_max and error 0) or \(self.prev_output ext_min and error 0):# 外部限制激活停止积分passelse:self.integral error * self.tselse:# 自身限幅检查if not (self.prev_output self.output_max and error 0) and \not (self.prev_output self.output_min and error 0):self.integral error * self.tselse:self.integral error * self.tsi_term self.ki * self.integral# 微分项对PV微分避免设定值突变冲击d_term 0.0if self.ts 0:d_term -self.kd * (process_variable - self.prev_pv) / self.ts# 总输出output p_term i_term d_term# 限幅output max(self.output_min, min(self.output_max, output))# 更新状态self.prev_error errorself.prev_output outputself.prev_pv process_variablereturn outputdef reset(self):重置控制器状态self.integral 0.0self.prev_error 0.0self.prev_output 0.0self.prev_pv 0.0# # 3. 风机模型领域模型# class BlowerModel:风机非线性模型 —— 领域模型def __init__(self, config: BlowerConfig, sample_time: float 1.0):self.config configself.ts sample_time# 状态变量self.current_freq: float config.min_freqself.current_pressure: float 0.0self.current_flow: float 0.0self.state: float 0.0 # 一阶惯性状态# 历史记录self.power_history: Deque[float] deque(maxlen1000)def step(self, freq_cmd: float) - Tuple[float, float, float]:执行一个仿真步长返回(频率, 风压, 功率)# 限幅target_freq max(self.config.min_freq,min(self.config.max_freq, freq_cmd))# 一阶惯性响应模拟风机启动/停止延迟tau self.config.response_timeself.state (target_freq - self.state) * self.ts / tauself.current_freq self.state# 非线性特性风压 ∝ 频率²self.current_pressure self.config.pressure_coeff * (self.current_freq ** 2)# 流量简化模型self.current_flow self.config.flow_coeff * self.current_freq# 功率计算power (self.config.rated_power * (self.current_freq / self.config.max_freq) ** 3) / \self.config.efficiencyself.power_history.append(power)return self.current_freq, self.current_pressure, powerdef get_efficiency(self) - float:计算当前运行效率if self.current_freq 0:return 0.0ideal_power self.config.rated_power * (self.current_freq / self.config.max_freq) ** 3actual_power self.power_history[-1] if self.power_history else ideal_powerreturn (ideal_power / actual_power) * 100 if actual_power 0 else 0def reset(self):重置风机状态self.current_freq self.config.min_freqself.current_pressure 0.0self.current_flow 0.0self.state 0.0self.power_history.clear()# # 4. 曝气池模型领域模型# class AerationTankModel:曝气池氧转移模型 —— 领域模型def __init__(self, config: AerationConfig, sample_time: float 1.0):self.config configself.ts sample_time# 状态变量self.do_concentration: float 2.0 # 初始DO (mg/L)self.temperature: float 20.0 # 水温 (℃)self.our_current: float config.our_base # 当前耗氧速率# 滞后环节模拟DO分析仪响应self.do_sensor_lag: float 0.0self.sensor_tau: float 15.0 # 传感器滞后时间常数 (s)# 扰动self.our_disturbance: float 0.0self.temp_disturbance: float 0.0self.random_gen np.random.RandomState(42)def _calc_saturation_do(self) - float:计算饱和溶解氧考虑温度和盐度# 温度修正Vant Hoff方程sat_do self.config.sat_do * (self.config.temp_coeff ** (self.temperature - 20.0))# 盐度修正简化sat_do * self.config.betareturn sat_dodef _calc_kla(self, pressure: float) - float:计算氧转移系数与风压正相关# 风压越高气泡越细传质面积越大kla self.config.kl_a_base * (1.0 0.1 * (pressure - 20.0) / 10.0)# 污水修正kla * self.config.alpha# 温度修正kla * (self.config.temp_coeff ** (self.temperature - 20.0))return max(0.01, kla)def step(self, pressure: float, our_disturbance: float 0.0,temp_disturbance: float 0.0) - float:执行一个仿真步长返回DO测量值# 更新扰动self.our_disturbance our_disturbanceself.temp_disturbance temp_disturbanceself.temperature 20.0 temp_disturbance self.random_gen.normal(0, 0.1)# 当前耗氧速率our self.config.our_base self.our_disturbance self.random_gen.normal(0, 0.01)self.our_current our# 饱和DOsat_do self._calc_saturation_do()# 氧转移系数kla self._calc_kla(pressure)# 氧转移速率 (mg/L/min)transfer_rate kla * (sat_do - self.do_concentration)# DO变化率 (mg/L/s)d_do (transfer_rate / 60.0) - our# 更新DO浓度self.do_concentration d_do * self.ts# 限幅self.do_concentration max(0.0, min(12.0, self.do_concentration))# 传感器滞后一阶惯性self.do_sensor_lag (self.do_concentration - self.do_sensor_lag) * self.ts / self.sensor_tau# 添加测量噪声noise self.random_gen.normal(0, 0.02) # 0.02mg/L噪声measured_do self.do_sensor_lag noisereturn measured_dodef reset(self):重置曝气池状态self.do_concentration 2.0self.temperature 20.0self.our_current self.config.our_baseself.do_sensor_lag 0.0self.our_disturbance 0.0self.temp_disturbance 0.0# # 5. 扰动发生器工厂模式# class DisturbanceGenerator:扰动发生器 —— 工厂模式def __init__(self, seed: int 42):self.random_gen np.random.RandomState(seed)def generate_our_disturbance(self, t: float) - float:生成耗氧速率扰动# 基础扰动随机波动base self.random_gen.normal(0, 0.01)# 周期性扰动模拟进水负荷变化periodic 0.03 * math.sin(2 * math.pi * t / 1800.0) # 30分钟周期# 阶跃扰动模拟冲击负荷step 0.0if 600 t 900: # 10~15分钟处step 0.05if 1800 t 2100: # 30~35分钟处step -0.03return base periodic stepdef generate_temp_disturbance(self, t: float) - float:生成温度扰动# 缓慢升温/降温trend 2.0 * math.sin(2 * math.pi * t / 3600.0) # 1小时周期# 随机波动noise self.random_gen.normal(0, 0.2)return trend noisedef generate_load_change(self, t: float) - float:生成进水负荷变化if 1200 t 1500: # 20~25分钟处return 1.5 # 负荷增加50%return 1.0# # 6. 串级控制器聚合根# class CascadeController:DO-风压串级控制器 —— 聚合根def __init__(self, do_config: DOControlConfig,blower_config: BlowerConfig,tank_config: AerationConfig):self.do_config do_configself.blower_config blower_configself.tank_config tank_config# 创建PID控制器self.main_pid PIDController(do_config.do_kp, do_config.do_ki, do_config.do_kd,do_config.do_output_min, do_config.do_output_max,do_config.anti_windup, do_config.sampling_time)self.slave_pid PIDController(do_config.press_kp, do_config.press_ki, do_config.press_kd,do_config.press_output_min, do_config.press_output_max,do_config.anti_windup, do_config.sampling_time)# 创建被控对象self.blower BlowerModel(blower_config, do_config.sampling_time)self.tank AerationTankModel(tank_config, do_config.sampling_time)# 扰动发生器self.disturbance DisturbanceGenerator()# 数据存储self.time_data: List[float] []self.do_sp_data: List[float] []self.do_pv_data: List[float] []self.press_sp_data: List[float] []self.press_pv_data: List[float] []self.freq_data: List[float] []self.power_data: List[float] []self.main_output_data: List[float] []self.slave_output_data: List[float] []self.our_data: List[float] []self.temp_data: List[float] []print(*60)print( 风机溶解氧串级控制仿真系统)print( 基于哈尔滨工程大学《工业过程控制》)print(*60)self._print_config()def _print_config(self):打印配置信息print(f\n[CONFIG] 控制参数:)print(f 主环(DO): Kp{self.do_config.do_kp}, Ki{self.do_config.do_ki}, Kd{self.do_config.do_kd})print(f 副环(风压): Kp{self.do_config.press_kp}, Ki{self.do_config.press_ki}, Kd{self.do_config.press_kd})print(f DO设定值: {self.do_config.do_sp} mg/L)print(f 采样周期: {self.do_config.sampling_time}s)print(f\n[CONFIG] 设备参数:)print(f 风机: {self.blower_config.rated_power}kW, {self.blower_config.min_freq}-{self.blower_config.max_freq}Hz)print(f 曝气池: {self.tank_config.volume}m³, 水深{self.tank_config.depth}m)print(f 饱和DO: {self.tank_config.sat_do}mg/L)def run_cascade_control(self, duration: float 3600.0) - Dict:运行串级控制仿真print(f\n 开始串级控制仿真 (时长: {duration/60:.0f}分钟)...)print(-*60)steps int(duration / self.do_config.sampling_time)for step in range(steps):t step * self.do_config.sampling_time# 1. 生成扰动our_dist self.disturbance.generate_our_disturbance(t)temp_dist self.disturbance.generate_temp_disturbance(t)load_factor self.disturbance.generate_load_change(t)# 2. 获取DO测量值do_pv self.tank.step(self.blower.current_pressure, our_dist, temp_dist)# 3. 主环DO环计算main_output self.main_pid.compute(self.do_config.do_sp, do_pv,external_limit(self.do_config.do_output_min, self.do_config.do_output_max))# 4. 副环风压环设定值 主环输出press_sp main_output# 5. 获取风压测量值_, press_pv, power self.blower.step(self.slave_pid.compute(press利用AI解决实际问题如果你觉得这个工具好用欢迎关注长安牧笛