速冻库温度控制:BP神经网络PID自整定与MATLAB/MCU落地
2026/9/17 12:09:31 网站建设 项目流程

简介:这份资源是一篇面向控制工程、机电一体化与嵌入式方向学习者的技术论文,主题为基于BP神经网络PID算法的快速冷冻系统设计,适合希望将神经网络与经典控制理论结合落地的本科高年级学生、研究生及工程技术人员参考。压缩包内仅含1个PDF文件,整体约1.71MB,即该论文的完整排版版本,便于打印精读与标注。论文以STM32F429为控制核心,系统讲解温度采集电路(NTC热敏电阻、桥式电路与AD12864组合)的设计与电桥输出电压推导,剖析比例、积分、微分三环节的作用,并说明BP神经网络如何在线自整定KP、KI、KD参数,取代人工试错;同时给出PWM占空比调节压缩机实现快速降温与恒温保持的完整控制逻辑。内容包含系统结构图、温度检测公式与PID原理图,可作为课程设计、毕业设计或制冷温控项目的算法建模与硬件选型参考。目前已有82人学习下载,适合需要掌握智能温控算法实现思路的读者。

1. 从速冻库温度曲线说起:为什么普通 PID 在快速冷冻上会“翻车”

一条典型的速冻隧道温度曲线是这样垮掉的:进料瞬间热负荷从 3 kW 跳到 18 kW,库温从 -32℃ 抬到 -21℃,普通 PID 看到偏差直接满输出,压缩机拉到 100%,库温冲过头砸到 -38℃,然后回弹、再冲、再回弹,一条曲线拉出三个波峰。操作工只能把 Ki 调小,把 Kp 压到 0.3,结果是降温变慢,食品中心温度在 30 分钟内穿不过 -18℃,速冻变成了慢冻,细胞里的冰晶长得又大又粗。

问题不在于 PID 本身,而在于快速冷冻这个对象在工况内是强非线性、大滞后、参数漂移的:蒸发器结霜让换热系数掉三成,库内风道被货架遮挡让纯滞后从 20 s 变成 40 s,压缩机变频下限和电子膨胀阀的步进特性又给执行器加了死区。一组固定 Kp、Ki、Kd 不可能同时覆盖空库预冷、满载速冻和除霜后复温三种状态。

BP 神经网络 PID 的做法是把 PID 的三个参数从“常数”变成“在线整定的函数”:用偏差 e 和偏差变化率 ec 作为输入,网络输出 Kp、Ki、Kd,控制量仍由增量式 PID 公式算出来。这样既保留了 PID 在工程上可解释、可限幅、可回退的优点,又让参数随工况自己漂。这篇文章讲清楚三件事:对象怎么建模、MATLAB 里最小可跑的代码长什么样、下位机部署和调参时哪些参数一旦设错就必然振荡。

2. 快速冷冻被控对象建模与 BP 神经网络 PID 的结构拆解

2.1 速冻库温度为什么能用一阶惯性加纯滞后近似

速冻库的热过程可以拆成三层:蒸发器侧制冷量、库内空气与货架的对流换热、食品内部的导热。对控制回路来说,把库温作为输出、把制冷量(电子膨胀阀开度或压缩机频率)作为输入,在 ±5℃ 的小偏差范围内,动态特性非常接近一阶惯性加纯滞后:

$$G(s)=\frac{K}{Ts+1}e^{-\tau s}$$

工程上常用的取值区间是:时间常数 T 取 120~300 s,纯滞后 τ 取 20~40 s,对象增益 K 为负——制冷量增大,库温下降。这里的 τ 不是理论推导出来的,是冷风机出口到回风温度传感器之间的空气混合时间加上传感器本身的响应时间。传感器如果装在回风口而离蒸发器只有 30 cm,τ 会小到 10 s,但测到的是局部温度,整库温度反而更滞后,这是现场最容易忽略的一处。

离散化直接用零阶保持器近似:

参数计算式典型值(Ts=1 s)
极点系数 aexp(-Ts/T)T=180 s 时 a≈0.9945
输入系数 bK·(1-a)K=-0.85 时 b≈-0.0047
滞后步数 dround(τ/Ts)τ=25 s 时 d=25

采样周期 Ts 的选择有个硬约束:Ts 不能大于对象时间常数的 1/20,否则纯滞后会被采样成 1~2 个步长,仿真里的稳定裕度完全是假的。Ts=1 s 对 T=180 s 来说裕度足够,下位机跑 500 ms 也行,但再快就没必要了,因为阀和变频器的机械响应本身就在 1 s 量级。

2.2 BP 神经网络 PID 回路的信号流与控制量公式

控制回路的信号流是:温度传感器 → 偏差 e(k) 和偏差变化率 ec(k) → 归一化 → BP 网络前向 → 输出 Kp、Ki、Kd → 增量式 PID → 限幅 → 执行器。控制律用增量式而不是位置式:

$$\Delta u(k)=K_p[e(k)-e(k-1)]+K_i e(k)+K_d[e(k)-2e(k-1)+e(k-2)]$$

选增量式有三个具体理由。执行器是电子膨胀阀这类积分型器件,给增量比给绝对开度更自然;输出限幅可以直接作用在 Δu 的累加结果上;Ki 通道即使短时偏大,也不会像位置式那样把积分项累到几百,除霜结束回温时不会出现长时间的积分退饱和。

网络输入为什么只用 e 和 ec 两维,而不是把设定值、货温、阀位都塞进去?因为输入维数每增加一维,在线反传的计算量和样本相关性都上升,而 BP-PID 的整定目标是三个参数随工况的慢漂移,不是精确拟合某个非线性函数。2-5-3 的规模(约 30 个权重和偏置)在 Cortex-M4 上单次前向加反传不到 200 个浮点乘加,1 s 周期跑绰绰有余。

2.3 隐藏层激活函数与输出层非负约束的选型边界

隐藏层用 tanh,不用 ReLU。原因是在线反传时输入样本会落在网络输入空间的边缘(比如化霜后偏差突然到 +15℃),ReLU 在这一侧导数为零,权重会整片“死掉”,而 tanh 在饱和区仍有小导数,网络能慢慢爬回来。

输出层必须做非负约束。Kp、Ki、Kd 在制冷对象里不能取负值,负的 Kp 意味着偏差越大输出越小,直接反向。做法是在输出层加 sigmoid,再乘以各自的上限:

% 输出层:sigmoid 保证非负,再按上限缩放 o = 1./(1+exp(-(w2*h + b2))); % o 为 3x1,范围 (0,1) Kp = KP_MAX * o(1); Ki = KI_MAX * o(2); Kd = KD_MAX * o(3);

三个上限不是拍脑袋定的,来源于经验整定:先用手动 PID 在空库预冷工况下整出一组能稳定的参数,把 Kp 取这组值的 2 倍作为 KP_MAX,Ki 取 1.5 倍,Kd 取 3 倍。这样网络有足够的整定空间,又不会因为一次误更新把参数推到发散区。相对的,如果把输出层写成线性再加 abs(),参数会在零点附近来回跳,实测曲线表现为周期性的小幅抖动。

3. MATLAB 里跑通 BP 神经网络自整定 PID 的最小仿真

3.1 一份可直接运行的 MATLAB 脚本

下面的脚本把对象、增量式 PID、BP 在线整定和限幅串在一个循环里,改几个常数就能对着自己的工况试。

% bp_pid_freezer.m —— 速冻库温度的 BP 神经网络自整定 PID 仿真 clear; clc; close all; % ---------- 1. 被控对象:一阶惯性 + 纯滞后 ---------- Ts = 1; % 采样周期 (s) Tobj = 180; % 对象时间常数 (s) Kobj = -0.85; % 对象增益,负号表示制冷量增大则库温下降 tau = 25; % 纯滞后 (s) N = 1200; % 仿真步数 r = -18; % 设定库温 (℃) a = exp(-Ts/Tobj); b = Kobj*(1-a); d = round(tau/Ts); y = zeros(1, N+1); y(1) = -5; % 初温 -5 ℃ u = zeros(1, N+1); % 控制量,0~100 表示阀开度百分比 % ---------- 2. BP 网络 2-5-3 ---------- sigmoid = @(x) 1./(1+exp(-x)); eta = 0.06; % 学习率 alpha = 0.05; % 动量项,抑制权重振荡 w1 = 0.6*randn(5,2); b1 = 0.2*randn(5,1); % 输入层->隐藏层 w2 = 0.6*randn(3,5); b2 = 0.2*randn(3,1); % 隐藏层->输出层 dw1 = 0; db1 = 0; dw2 = 0; db2 = 0; KP_MAX = 2.0; KI_MAX = 0.03; KD_MAX = 1.0; % 三个参数上限 s = -0.02; % 灵敏度近似系数(含对象增益符号) e1 = 0; e2 = 0; % e(k-1), e(k-2) logKp = zeros(1,N); logKi = logKp; logKd = logKp; for k = 1:N e = y(k) - r; % 正向偏差:库温高于设定为正 ec = e - e1; % ---- 前向:归一化输入 -> 隐藏层 -> 输出层 ---- x = [e/10; ec/10]; h = tanh(w1*x + b1); o = sigmoid(w2*h + b2); Kp = KP_MAX*o(1); Ki = KI_MAX*o(2); Kd = KD_MAX*o(3); % ---- 增量式 PID ---- du = Kp*(e-e1) + Ki*e + Kd*(e-2*e1+e2); u(k) = min(max(u(k-1) + du, 0), 100); % 执行器限幅 0~100% % ---- 被控对象递推 ---- ud = 0; if k-d >= 1, ud = u(k-d); end y(k+1) = a*y(k) + b*ud; % ---- 反向传播:以 0.5*e^2 为代价函数 ---- dJ = [e*(e-e1); e*e; e*(e-2*e1+e2)] * s; % 对 Kp/Ki/Kd 的灵敏度 do = -dJ .* (o .* (1-o)); % 输出层 delta dh = (w2' * do) .* (1 - h.^2); % 隐藏层 delta(tanh 导数) dw2 = eta*(do*h') + alpha*dw2; db2 = eta*do + alpha*db2; dw1 = eta*(dh*x') + alpha*dw1; db1 = eta*dh + alpha*db1; w2 = w2 + dw2; b2 = b2 + db2; w1 = w1 + dw1; b1 = b1 + db1; logKp(k) = Kp; logKi(k) = Ki; logKd(k) = Kd; e2 = e1; e1 = e; end % ---------- 3. 出图 ---------- t = (0:N)*Ts/60; figure; subplot(2,1,1); plot(t, y, 'LineWidth', 1.2); hold on; plot([0 N*Ts/60], [r r], '--'); grid on; xlabel('时间 (min)'); ylabel('库温 (℃)'); title('BP-PID 速冻库温响应'); subplot(2,1,2); plot(t(1:N), logKp, t(1:N), logKi*100, t(1:N), logKd); grid on; legend('Kp','Ki×100','Kd'); xlabel('时间 (min)'); title('在线整定的 PID 参数');

3.2 反向传播里那个 s 系数是怎么来的

dJ这一行是整个脚本最容易写错的地方。代价函数取 J = 0.5·e²,对 Kp 求偏导要用链式法则:

$$\frac{\partial J}{\partial K_p}=e\cdot\frac{\partial y}{\partial u}\cdot\frac{\partial u}{\partial K_p}=e\cdot b\cdot[e(k)-e(k-1)]$$

真实的对象增益 b 在控制器里是未知的,在线辨识又太贵。工程做法是用一个固定的小系数 s 近似 b 的量级,符号按对象方向取。制冷对象的 b 是负的,所以 s 取负值,代码里写成 -0.02。s 的绝对值只影响参数收敛速度,不影响收敛方向,取对象增益量级的 1/50 到 1/20 就够了。

s设得过大,Kp 会在两三个采样周期内顶到上限,然后被限幅卡住来回撞,曲线表现为等幅振荡;设得过小,网络更新量被淹没在浮点误差里,跑十分钟 Kp 只动了 0.01,等于网络没起作用。判断方法很简单:把 logKp 画出来,正常的曲线应该是在阶跃后 30~60 s 内明显抬升,然后随偏差减小缓慢回落,全程平滑。

3.3 仿真里必须调的四个参数

参数作用调大调小
eta(学习率)权重更新步长收敛快,易振荡平稳,跟踪慢
alpha(动量)抑制权重方向抖动抗噪好,易过冲响应干净,易抖
隐藏层节点数拟合能力表达强,易过拟合单工况泛化好,欠拟合
Ts(采样周期)离散精度计算省,纯滞后失真精度高,MCU 负担重

隐藏层节点数从 5 加到 8 通常不会有明显改善,反而在化霜回温这种未见过的工况下参数漂得更厉害,因为多余节点把训练时常遇到的空库预冷工况记死了。我的经验是:只要被控对象的工况不超过 3 种,5 个节点足够。

提示:仿真跑之前先把 KP_MAX、KI_MAX、KD_MAX 用一组手调 PID 参数标定好,否则网络输出的初值就是随机的,前 60 s 的曲线没有参考价值。

4. 从 MATLAB 搬到 MCU:下位机实现的资源账与限幅细节

4.1 用 C 语言重写前向与反传

下位机代码不需要 MATLAB 的矩阵语法,展开成循环更省。核心循环如下,权重用 float 存,Cortex-M4 带 FPU 时单次迭代约 30 μs。

/* bp_pid.c —— MCU 上的 2-5-3 BP-PID,采样周期 1 s */ #include <math.h> #define IN_N 2 #define HID_N 5 #define OUT_N 3 static const float K_MAX[OUT_N] = { 2.0f, 0.03f, 1.0f }; /* Kp/Ki/Kd 上限 */ static float w1[HID_N][IN_N], b1[HID_N]; static float w2[OUT_N][HID_N], b2[OUT_N]; static float sigmoidf(float x) { return 1.0f / (1.0f + expf(-x)); } /* temp: 当前库温, set: 设定温度, out: 返回 0~100 的阀开度 */ void bp_pid_step(float temp, float set, float *out) { static float e1 = 0.0f, e2 = 0.0f, uprev = 0.0f; static uint16_t tick = 0; float e = temp - set; /* 正向偏差,单位 ℃ */ float ec = e - e1; /* --- 前向 --- */ float x[IN_N] = { e * 0.1f, ec * 0.1f }; /* 归一化必须和仿真一致 */ float h[HID_N], o[OUT_N]; for (int i = 0; i < HID_N; i++) { float z = b1[i]; for (int j = 0; j < IN_N; j++) z += w1[i][j] * x[j]; h[i] = tanhf(z); } for (int i = 0; i < OUT_N; i++) { float z = b2[i]; for (int j = 0; j < HID_N; j++) z += w2[i][j] * h[j]; o[i] = sigmoidf(z); } /* --- 增量式 PID --- */ float kp = K_MAX[0]*o[0], ki = K_MAX[1]*o[1], kd = K_MAX[2]*o[2]; float du = kp*(e-e1) + ki*e + kd*(e - 2.0f*e1 + e2); float u = uprev + du; if (u > 100.0f) u = 100.0f; if (u < 0.0f) u = 0.0f; *out = u; uprev = u; /* --- 反向传播:每 5 拍做一次,省 CPU --- */ if (++tick % 5 == 0) { const float s = -0.02f; float dJ[OUT_N] = { e*(e-e1)*s, e*e*s, e*(e-2.0f*e1+e2)*s }; float d_o[OUT_N], d_h[HID_N]; for (int i = 0; i < OUT_N; i++) d_o[i] = -dJ[i] * o[i] * (1.0f - o[i]); for (int i = 0; i < HID_N; i++) { float sum = 0.0f; for (int j = 0; j < OUT_N; j++) sum += w2[j][i] * d_o[j]; d_h[i] = sum * (1.0f - h[i]*h[i]); /* tanh 导数 */ } const float eta = 0.06f, alpha = 0.05f; for (int i = 0; i < OUT_N; i++) { for (int j = 0; j < HID_N; j++) w2[i][j] += eta*d_o[i]*h[j] + alpha*(eta*d_o[i]*h[j]); b2[i] += eta*d_o[i]; } for (int i = 0; i < HID_N; i++) { for (int j = 0; j < IN_N; j++) w1[i][j] += eta*d_h[i]*x[j]; b1[i] += eta*d_h[i]; } } e2 = e1; e1 = e; }

4.2 归一化系数必须和仿真保持一致

x[0] = e * 0.1f这一行看着不起眼,却是仿真能跑、上机就炸的头号原因。MATLAB 里写的是e/10,C 里如果写成e/20,网络的输入尺度变了,权重是按旧尺度训练的,输出层 sigmoid 直接进饱和区,Kp 输出卡在 0 或者上限。改归一化系数等于改网络输入的定义,改完必须重新跑仿真整定,或者干脆在设备上用初始权重重新在线学习。

4.3 MCU 上的资源预算与限幅位置

项目建议值说明
采样周期500 ms ~ 1 s不超过对象时间常数的 1/20
网络规模2-5-3权重+偏置 28 个 float,约 112 B
反传频率每 5 ~ 10 拍一次参数本身慢漂,省下的算力给采样
前向耗时Cortex-M4 @72 MHz 约 8 μs带 FPU,不加反传
PID 输出限幅0 ~ 100%对应阀开度或压缩机频率归一化值
偏差死区±0.2 ℃防止阀在稳态反复步进

限幅必须放在 PID 输出之后、累加之前,不能放在网络输出上。把限幅加在 Kp、Ki、Kd 上,会让网络误以为参数被“接受”了,反传继续往限幅边界推,退出限幅后参数会保持在上限附近,回温工况下表现为输出长时间贴顶。放在 u 上则是物理约束,PID 结构本身能感知到饱和。

另外,阀位死区是执行器的机械特性,不是控制算法能消除的。库里温度稳定在设定值 ±0.2℃ 以内时,直接把 u 保持住,不再累加 Δu,否则电子膨胀阀会以每 1 s 一步的速度来回走,一晚上步进上千次,阀杆磨损比什么都快。

5. 从仿真到实测:验证 BP 神经网络 PID 的三个动作和一组排错对照

仿真曲线好看不代表设备能跑,因为仿真里的对象参数是自己写的。上机之后我一般按三步验证,顺序不能换。

第一步是冻结网络、只用固定 PID 参数跑阶跃。把 eta 置 0,把 Kp、Ki、Kd 手动设成 KP_MAX 的一半,从 -5℃ 预冷到 -18℃,记录超调量和调节时间。这一步的目的是拿到真实对象的 T 和 τ:把响应曲线的 63% 上升点时间和纯滞后段的水平距离量出来,反推 T 和 τ,回填到仿真里。很多现场仿真和实测对不上,就是因为 T 按 180 s 填的,实际是 90 s。

第二步是放开网络但把 eta 压到 0.01。观察 logKp 的上位机曲线,正常行为是在阶跃后 1~2 分钟内 Kp 有可见抬升,然后随偏差收敛缓慢下降。如果 Kp 曲线是一条直线,说明 s 太小或者输入归一化不对;如果 Kp 一路顶到 KP_MAX,先把 eta 减半再试。

第三步才是满载速冻加除霜的完整工况循环。这一步专门看除霜结束复温那 5 分钟,因为库温会从 -35℃ 冲到 -12℃,偏差 20℃ 以上,是网络输入空间的最边缘。

排错时对照下面这张表比盲调快得多。

现象最可能的原因处理
参数一路顶到上限eta 偏大或 s 绝对值偏大eta 减半,s 缩到 1/2
库温等幅振荡、Kp 平滑KP_MAX 定得过高KP_MAX 降到当前 0.6 倍
反向降温慢Ki 上限过小或死区过大KI_MAX 提到 1.5 倍,死区收到 0.1℃
除霜后长时间贴在 -12℃积分项退出饱和慢复温开始时把 u 直接置 0,重置 e1/e2
参数长时间几乎不动学习率太小或浮点精度不够eta 提到 0.1,确认用的是 float 而非定点的残差
仿真好、上机抖归一化系数不一致逐位核对 x 的缩放系数

最后一个技巧关于回退:BP 网络在线整定本质上是把控制器的稳定性押在了一个自适应环节上,工程上必须留后门。我的做法是在上位机上保留一组“验证过的固定 PID 参数”,设备连续 3 分钟偏差超过 3℃ 或者输出被限幅超过 60 s,就自动切回固定参数并把网络权重冻结,同时在 HMI 上报警。这样即使网络因为长期运行漂进了坏区域,也不会把一库货冻坏。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询