前几年我帮朋友做了一套电加热箱的温控改造最初用的就是常规PID。结果在冬季外温波动大的时候温度经常超调2-3℃而且加热器频繁启停实际能耗并不低。后来我换成了模型预测控制MPC再叠加卡尔曼滤波处理传感器噪声和未知热扰动实测下来温度偏差基本控制在0.3℃以内设备动作次数也大幅下降。这篇就是把我这次“MPC卡尔曼滤波调节空调加热器/室内温度”的完整建模、代码设计和调试过程做个记录适合正在学MPC、做温控仿真、或者准备把Matlab里的控制器落地到实际暖通项目里的朋友参考。先说明一点你在网上搜“MPC”时可能会看到存储管理平台之类的同名缩写和模型预测控制完全是两码事。本文只聊控制领域的Model Predictive Control场景是空调制热/电加热器向房间供热语言是Matlab。目标不是堆一堆公式而是让你看完之后能照着代码自己跑一遍并且知道每个参数为什么这么取。1. 项目思路为什么是MPC配合卡尔曼滤波1.1 温控场景里“好算法”到底解决了什么问题房间温度控制并不是什么新鲜事但把它做精细并不容易。最基本的温控方式是回差控制温度低于下限就开加热高于上限就关。这种方式实现简单缺点也很明显——温度在上下限之间来回摆动传感器噪声还会导致继电器频繁跳动加热器寿命和舒适度都很差。传统PID通过比例、积分、微分来调节加热功率在没有大干扰的情况下表现尚可。但它本质上是“看到误差才反应”外温骤降导致房间开始降温PID要等温度真的偏离目标后才慢慢增加输出因此天然存在滞后。而且PID参数往往在一个工况下标定到了另一个外部温度和房间负荷下又得重新整定。到这里MPC的价值就凸现出来了。MPC会用一个房间热动态模型去预测未来十几分钟甚至更长时间的温度走势然后根据这些预测算出接下来一组最优的加热功率序列每走一步就重新算一遍这就是“滚动优化”。因为控制器拥有了对未来的预测能力它可以提前加大功率来抵抗即将到来的降温而不是事后追误差。更重要的是MPC能非常自然地处理约束——加热器功率不能为负、不能超过最大功率、温度最好不超过舒适区上界这些约束在PID里很难统一处理在MPC里只是优化问题里的几个不等式。1.2 卡尔曼滤波在这里补了什么短板有了模型预测控制器就依赖模型的“当前状态”来推算未来。而实际系统中室温传感器信号是有噪声的而且房间本身还有开门、人员活动、太阳辐射这类无法建模的干扰。如果你直接把带噪声的温度测量值送给MPC控制器会把噪声当成真实温度变化进而产生不必要的功率波动甚至触发振荡。卡尔曼滤波的价值就在这里。它不像低通滤波那样简单地把信号抹平而是结合一个状态模型并对测量噪声和模型误差分别设置协方差参数从而给出一个“最小均方误差意义下最优”的状态估计。更妙的是我们可以把未知热流干扰也扩展成状态的一部分让卡尔曼滤波在线估计出“当前房间多了一股300W的异常散热”MPC拿到这个估计值后就能提前补偿。我把这套结构总结成一句话MPC负责“往后看多远、怎么控制”卡尔曼滤波负责“当前到底处在一个什么状态”两者一配合控制品质才真正上了一个档次。2. 系统建模把房间变成状态空间方程2.1 房间热动态的等效集总参数模型温控对象本身不是什么复杂的机械系统但如果你把墙体内部分布、空气流动、窗户传热全部建模复杂度会爆炸而且对控制算法来讲往往没必要。工程上常用的是集总参数模型也就是把整个房间等效成一个热容C和一个热阻1/UA。具体物理意义是房间内空气、家具、内墙等总蓄热能力用等效热容C表示单位是J/K墙壁、窗户等围护结构对室内外温差的热传导能力用UA表示单位是W/K。根据能量守恒房间温度T的变化满足C·(dT/dt) -UA·(T - T_out) P_heat Q_dist其中T_out是室外温度P_heat是加热器功率Q_dist是未建模的干扰热流比如开门冷风渗透、人员散热量、设备发热等。这个一阶微分方程已经把问题简化得足够好了实际辨识后精度在温控场景下完全够用。在Matlab里我会直接把这个连续模型写成状态空间形式再用c2d离散化。因为MPC和卡尔曼滤波都在离散域里工作采样周期选多大很重要。对于室温这类大惯性对象采样周期取30秒到60秒都很正常本文代码统一用60秒。2.2 状态空间表达式和离散化把模型写成状态方程dx/dt A_c·x B_c·u在这个问题里我取状态x [T]室内温度控制输入u [P_heat]加热器功率单位W外部温度T_out作为可测干扰单独处理未知热流Q_dist作为扩展状态交给卡尔曼滤波估计。连续状态空间是A_c -UA/CB_c 1/C输入除了控制量P_heat还有个外部温度通道系数是UA/C。因为后面要估计未知热流我会把状态扩成两个x1Tx2Q_dist。这样离散化之后的状态矩阵就是2×2的增广矩阵。假设采样时间Ts60sMatlab代码里直接这样离散化C_room 360000; % J/K房间等效热容 UA 80; % W/K围护结构传热系数 A_c -UA / C_room; B_c 1 / C_room; % 连续系统输入通道之一为控制功率另一个为外温前馈 sys_c ss(A_c, B_c, 1, 0); % 离散化 sys_d c2d(sys_c, Ts, zoh); Ad sys_d.a; Bd sys_d.b;扩展后的状态方程写成A_aug [Ad, Bd; 0, 1]; B_aug [Bd; 0]; % 控制输入P_heat C_aug [1, 0]; % 观测量仍是室内温度这里的思路是把未知热流Q_dist当作一个随机游走变量它本身没有固定动态模型每个采样周期都可能变化变化大小由过程噪声协方差决定。这样卡尔曼滤波就能根据温度测量残差不断修正Q_dist的估计值。2.3 一套能直接用的模型参数很多同学卡在第一步就是不知道该填什么数。我给出我这次使用的典型参数你可以直接搬过去做仿真后续再做辨识替换成自己系统的数值。参数数值含义C_room360000 J/K等效热容约相当于60m³房间加家具UA80 W/K围护结构综合传热系数P_max3000 W加热器最大功率Ts60 s控制/采样周期初始室温10 ℃仿真起始温度设定温度22 ℃ → 24 ℃前30分钟22℃之后24℃按照这套参数房间时间常数τ C/UA 4500s也就是75分钟。升温初期3000W满功率时升温速率约0.5K/min符合现实中一间普通房间电加热的响应速度。仿真总时长我设为5400秒也就是90分钟足够看到一次设定值变化和一次未知扰动的全过程。3. 卡尔曼滤波让控制器看到“修正后的真实状态”3.1 增广状态估计的思路如果只做仿真温度当然可以直接读到但实际系统里传感器就是有噪声的。更麻烦的是房间还老有未知散热。让我用一个具体场景说明你正控制房间加热到22℃一切稳定。突然有人推开大门一股冷风进来相当于房间多了一个约500W的散热负荷。这个散热负荷模型里没有MPC也不知道。如果没有任何在线估计控制器只能等温度真正跌下去之后靠着反馈误差慢慢回调这期间室温可能已经掉到21℃了。但如果把Q_dist扩展成状态变量卡尔曼滤波会在每个采样周期根据温度测量偏差不断调整Q_dist的估计值。用不了几个周期它就能估计出“房间当前多了大约400-500W的额外负荷”。MPC在下一次滚动优化时会把这项负荷直接放进预测模型里于是控制器几乎能在温度明显下降前就提高加热功率来补偿。这就是我为什么强调“卡尔曼滤波不是简单平滑一下信号”的原因。它给MPC提供的不光是一个干净的T还有那些模型里缺失的干扰估计值。3.2 卡尔曼滤波递推公式和代码实现标准的离散卡尔曼滤波分为预测和更新两步。状态x是二维的[ T; Q_dist ]观测只有温度T预测步x_pred A_aug·x_est;P_pred A_aug·P_est·A_aug Q_kf;更新步K P_pred·C_aug / (C_aug·P_pred·C_aug R_kf);x_est x_pred K·(y_meas - C_aug·x_pred);P_est (eye(2) - K·C_aug)·P_pred;R_kf是测量噪声方差。温度传感器精度一般在±0.3℃左右转化成方差大概0.1到0.25我取0.1。Q_kf是过程噪声协方差矩阵这个要斟酌Q_kf diag([0.0001, 100]);第一项对应温度状态的模型误差相当于每个采样周期温度本身可能偏移0.01℃左右第二项对应未知热流Q_dist的变化能力方差100表示每个采样周期未知热流可能变化10W左右。这个数值如果设太小滤波对Q_dist的响应就很慢设太大估计容易跟着测量噪声跳。实际操作中我一般会做几次仿真对比看估计曲线在不失真的前提下能多快跟上突变负荷。你还可以直接调用Matlab自带的kalman函数但手写递推能让你更清楚地看到每一步在做什么而且方便以后扩展成无迹卡尔曼滤波或扩展卡尔曼滤波。新手我建议至少手写一遍。3.3 初始状态设置和滤波效果观察初始状态我故意设成和真实值有偏差比如T_est 15℃ 而真实T 10℃ Q_dist_est 0 而真实Q_dist 0初始协方差P_0 diag([4, 400])表示我对温度初始估计误差有信心到2℃左右对热流估计完全没有把握。卡尔曼滤波会根据前几个测量点的残差迅速拉近估计值这个收敛过程在仿真曲线上非常直观。仿真中我还会在30分钟时刻给Q_dist叠加一个-500W的突变负荷代表突然开门等扰动。你会发现卡尔曼滤波估计的Q_dist在大约5到8个采样周期内逐渐接近-500W而室温几乎看不出明显下跌因为MPC配合着把功率提上去了。这正是整套设计最精彩的地方。4. MPC控制器设计滚动优化里的约束和目标4.1 预测模型与约束设置丢开工具箱的手写MPC核心就是把预测方程写出来。假设当前时刻k我取预测时域Np30步控制时域Nc也暂时取30步工程上Nc小于Np可以减少计算量但演示代码里先保持两者相同避免复杂度。已知当前状态x_est和当前外温T_out那么未来30步的温度预测为Y F·x_est G·U F_out·T_out其中U[u(k), u(k1), ..., u(kNp-1)]是未来每个时刻的加热功率序列每行Y对应未来一步的室内温度预测值。矩阵F由状态矩阵A_aug幂次构成G是控制输入到未来输出的响应矩阵F_out是外温前馈项。这些矩阵在Matlab里用循环就能构造不依赖MPC Toolbox。约束方面加热器功率有物理边界0 ≤ u ≤ 3000 W这就是最典型的不等式约束。在quadprog里直接通过lb和ub设置。如果你想加输出约束比如室温不超过26℃那需要做软约束处理否则遇到不可行问题很头痛。演示代码里暂不加输出硬约束只约束控制量最大程度保证问题一定有解。4.2 目标函数怎么设计MPC每一轮要解的优化问题本质上是寻找最合适的未来控制序列U使得下面这个目标函数最小J Σ Q_p·(T_pred - T_ref)² W_r·u²第一项是温度跟踪误差让预测温度尽量贴近设定值第二项是控制量平方项代价很小主要防止优化问题奇异。Q_p 5表示温度偏差每1℃会产生5的惩罚W_r 0.0001对3000W功率来说相当于0.9的惩罚远小于温度误差的贡献所以控制器会优先满足温度跟踪同时也不会无意义地猛加功率。这个目标函数写成关于U的二次型就是标准的min 1/2·U·H·U f·U其中H 2·(G·Q_p·G W_r·I) f -2·G·Q_p·(T_ref_vec - F·x_est - F_out·T_out)这里T_ref_vec就是未来Np步的设定温度向量。因为采样周期60秒设定温度变化是缓慢的我直接把当前设定值或已知未来设定曲线填进去即可。在Matlab里调用quadprog求解U_opt quadprog(H, f, [], [], [], [], lb, ub, U_init);U_init设置成全1×u_prev或者上一步求出的U都可以quadprog的active-set算法在有初始点时收敛更快。4.3 滚动优化和反馈校正得到U_opt后只有第一个元素u(k)真正送给执行器下一时刻全部重新来过。这一步就是“滚动优化”或者叫“后退时域控制”。为何只取第一步因为预测模型并不是完美准确的卡尔曼滤波给出的状态估计也在每步更新。如果我现在一口气把未来30步的功率都执行完中间一旦出现扰动后面29步都是失效的。至少每步重新基于最新状态估计求解一次才能始终保持控制动作贴着真实系统走。在实际代码里主循环的结构是更新真实房间温度仿真中模拟真实对象读取传感器温度真实温度加噪声卡尔曼滤波得到T_est和Q_dist_est利用T_est和Q_dist_est构造MPC预测方程quadprog求解U_opt取U_opt(1)作为当前控制功率送给执行器温度误差和功率曲线都记录后最后画成图分析。5. Matlab代码实现从模型到仿真一条龙5.1 初始化参数与预测矩阵构造先把整个脚本的主干列出来。这段代码可以直接存成一个脚本文件跑通前提是你有Optimization Toolbox因为会用到quadprog。不需要MPC Toolbox这样门槛低很多。%% 参数初始化 clear; clc; close all; Ts 60; % 采样周期 60s C_room 360000; % 等效热容 J/K UA 80; % 传热系数 W/K P_max 3000; % 最大加热功率 W A_c -UA / C_room; B_c 1 / C_room; sys_c ss(A_c, B_c, 1, 0); sys_d c2d(sys_c, Ts, zoh); Ad sys_d.a; Bd sys_d.b; % 增广状态x [T; Q_dist]Q_dist为未知热流 A_aug [Ad, Bd; 0, 1]; B_aug [Bd; 0]; C_aug [1, 0]; B_feed [Bd * UA; 0]; % 外温前馈通道 Np 30; % 预测时域 nx 2; % 构造预测矩阵 F zeros(Np, nx); G zeros(Np, Np); F_out zeros(Np, 1); for i 1:Np F(i, :) C_aug * (A_aug^i); for j 1:i-1 G(i, j) C_aug * (A_aug^(i-j)) * B_aug; F_out(i) F_out(i) C_aug * (A_aug^(i-j)) * B_feed; end end这里构造G矩阵时j从1到i-1表示第j个控制量对未来第i步温度预测的贡献。因为是SISO系统矩阵规模不算大。Np取30quadprog求解一个30维二次规划非常轻松单步求解时间在毫秒级。5.2 卡尔曼滤波参数设置这一块我写成函数形式方便在主循环里调用。这里为了演示清晰我把递推公式直接写在主循环里没有封装成子函数你可以自行整理%% 卡尔曼滤波参数 Q_kf diag([0.0001, 100]); R_kf 0.1; x_est [15; 0]; % 初始估计温度误差故意给大 P_est diag([4, 400]);温度过程的噪声方差0.0001很小因为温度模型本身是一个确定性较强的过程Q_dist的方差给100表示它可能随时变化这样卡尔曼滤波对负荷突变的响应才会够快。传感器R_kf取0.1对应标准差约0.316℃的测量噪声。5.3 仿真场景搭建与控制主循环为了让效果明显我设计了一个90分钟的仿真场景中途包含两类事件第30分钟设定温度从22℃升到24℃第30分钟房间额外出现-500W的扰散热流模拟寒风吹入或多人进入导致的热量流失。室外温度T_out设定为%% 仿真场景 Nsim 90; % 90分钟 T_out 0 3 * sin((1:Nsim) / 30 * pi / 3); % 外温缓慢波动 T_ref [22 * ones(1, 30), 24 * ones(1, Nsim-30)]; Q_true [zeros(1, 30), -500 * ones(1, Nsim-30)];T_true初值取10℃也就是房间一开始比较冷MPC从满功率附近开始加热。主循环代码%% 主循环 T_true 10; u_prev 0; T_hist zeros(1, Nsim); u_hist zeros(1, Nsim); est_hist zeros(2, Nsim); meas_hist zeros(1, Nsim); for k 1:Nsim % 真实房间更新 w_proc sqrt(0.0001) * randn; T_true Ad * T_true Bd * (u_prev Q_true(k)) Bd * UA * T_out(k) w_proc; T_true max(T_true, -10); % 传感器测量含噪声 v_noise sqrt(R_kf) * randn; y_meas T_true v_noise; % 卡尔曼滤波预测更新 x_pred A_aug * x_est; P_pred A_aug * P_est * A_aug Q_kf; K P_pred * C_aug / (C_aug * P_pred * C_aug R_kf); x_est x_pred K * (y_meas - C_aug * x_pred); P_est (eye(nx) - K * C_aug) * P_pred; % MPC滚动求解用估计状态 current_e x_est(1); current_q x_est(2); x_mpc [current_e; current_q]; f_qp -2 * G * Q_p * (T_ref_vec - F * x_mpc - F_out * T_out(k)); lb zeros(Np, 1); ub P_max * ones(Np, 1); options optimoptions(quadprog, Display, off); U_opt quadprog(H_qp, f_qp, [], [], [], [], lb, ub, [], options); u_now U_opt(1); u_prev u_now; % 记录 T_hist(k) T_true; meas_hist(k) y_meas; est_hist(:, k) x_est; u_hist(k) u_now; end注意T_ref_vec在循环前要根据每个时刻k取未来Np步的参考轨迹。如果k Np超过序列长度就截尾并用最后值填充if k Np Nsim T_ref_vec T_ref(k1 : kNp); else T_ref_vec [T_ref(k1:end), T_ref(end)*ones(1, kNp-Nsim)]; end因为我在代码里设Np30Nsim90所以前60步都正常不会触发截尾。5.4 结果可视化绘制三张图室内温度曲线、卡尔曼滤波估计和真实值的对比、加热器功率曲线。figure; subplot(3,1,1); plot((1:Nsim)/60, T_hist, linewidth, 1.5); hold on; plot((1:Nsim)/60, T_ref, --, linewidth, 1); legend(室温, 设定温度); ylabel(温度 / ℃); grid on; subplot(3,1,2); plot((1:Nsim)/60, T_true, linewidth, 1.5); hold on; plot((1:Nsim)/60, meas_hist, color, [0.7 0.7 0.7]); plot((1:Nsim)/60, est_hist(1,:), --, linewidth, 1.5); legend(真实温度, 传感器测量, 卡尔曼估计); ylabel(温度 / ℃); grid on; subplot(3,1,3); stairs((1:Nsim)/60, u_hist, linewidth, 1.5); ylabel(加热功率 / W); xlabel(时间 / min); grid on;跑完之后你应当能看到温度由初始10℃逐渐上升约30分钟左右到达22℃附近30分钟后设定值升到24℃同时Q_true叠加-500W扰动温度会有一个短暂小幅下垂但MPC会在几分钟内把功率提上来室温迅速回到24℃附近。滤波后的曲线明显比原始测量曲线平滑而且估计值能跟踪真实温度。5.5 使用MPC Toolbox的快速实现方式如果你手头有MPC Toolbox代码还能更短。先建好被控对象的离散状态空间模型plant ss(A_aug, B_aug, C_aug, 0, Ts, Ts); mpcobj mpc(plant, Ts); mpcobj.PredictionHorizon 30; mpcobj.ControlHorizon 5; mpcobj.Weights.OutputVariables 5; mpcobj.Weights.ManipulatedVariablesRate 0.5; mpcobj.ManipulatedVariables(1).Min 0; mpcobj.ManipulatedVariables(1).Max P_max;然后把卡尔曼滤波维护的状态est_hist(:, k)传给sim函数或者在每一环手动调用mpcmove。Toolbox更方便的地方在于增量权重、输出软约束、模型失配鲁棒性设置都很成熟。但我还是建议你先跑通手写版本因为只有自己构造过G矩阵你才算真正理解MPC内部在算什么。6. 调试实录参数怎么调才能不翻车6.1 卡尔曼滤波协方差整定的经验卡尔曼滤波里Q_kf和R_kf的比值本质上决定了你更相信“模型预测”还是“传感器测量”。Q/R越小估计越平滑、抗噪越强但对快速突变的响应越慢Q/R越大估计越敏捷但噪声也容易被带进控制量。我在这个场景里的Q_kf取diag([0.0001, 100])R_kf取0.1。如果你发现温度估计有很明显的锯齿多半是R_kf太小或者Q_kf(1,1)太大。如果你发现Q_dist估计响应特别慢甚至无法跟踪500W的突变负荷就应该把Q_kf(2,2)调大比如调成400到1000让滤波器允许热流状态更快变化。另一个容易被忽略的地方是初始协方差P_0。如果P_0给太小卡尔曼滤波一开始会非常信任初始状态估计收敛会很慢。我特意给P_est初值设成diag([4, 400])对应2℃的温度不确定度和20W左右的负荷不确定度这样前几个周期就能快速收敛。6.2 MPC预测时域和控制时域的取值逻辑预测时域Np的选取直接决定控制器“看得多远”。这个距离应该覆盖系统动态的主要时间尺度。对于时间常数75分钟的房间Np30步、步长60秒对应预测30分钟已经覆盖了近半个时间常数跟踪设定值和扰动抑制表现不错。如果Np太小比如只有5步MPC就变成了一个近视眼约束和未来轨迹的作用都会大打折扣。控制时域Nc没必要和Np一样大。实际工程中把Nc设成3到10就够了因为后续的控制量变化对当前性能影响有限减小Nc还能降低优化问题维度。在手写版本里我为了简化让NcNp换成Toolbox时明显看到Nc5和Nc30的结果差别不大但求解速度更高。如果你发现MPC动作很激进、功率频繁波动优先加大ManipulatedVariablesRate权重或者减少Nc。如果你发现跟踪速度太慢、升温太磨叽则要适当调大Q_p或者检查你的模型参数是不是把房间热容设得太大了。6.3 常见问题速查表现象可能原因处理办法卡尔曼估计温度有锯齿Q/R比例偏大滤波过于相信测量增大R_kf或减小Q_kf(1,1)Q_dist估计跟不上突变负荷Q_kf(2,2)太小将Q_kf(2,2)调到400~1000温度超调明显MPC预测时域太短或R_w太小增加Np适当增大控制量权重加热功率高频波动控制增量没有惩罚使用增量权重或将Nc调小quadprog提示不可行约束设置过于严格检查是否加了过紧输出约束改用软约束系统静差迟迟消除不了模型参数有偏状态估计不准检查Q_dist是否被估计必要时做模型辨识6.4 一些值得说的实战细节第一个细节外温前馈的作用不可小觑。虽然卡尔曼滤波能估计Q_dist但它需要时间收敛。如果外温突然从0℃降到-5℃这算可测变化直接通过前馈通道加进预测模型MPC当拍就能反应。相比之下完全靠反馈去抵消外温变化那要等好几个采样周期。第二个细节控制量更新之间最好限幅。即使你的MPC在模型里约束了0到3000W实际加热器也未必能在1秒内完成满功率输出。真正的执行器有变化率限制。工程上我会额外加一个每步最大变化量比如300W/步实现方式是在优化问题里加入线性不等式约束或者简单点对quadprog求出的U_opt(1)做一次限速处理u_now max(min(U_opt(1), u_prev 300), u_prev - 300); u_now max(0, min(u_now, P_max));第三个细节仿真里虽然可以用真实温度做反馈但实际系统必须用卡尔曼滤波估计值。我代码里MPC的输入是est_hist(:,k)也就是x_est而不是T_true。这一点看着不起眼却是整套控制方案能不能落地到现场的关键。另一个我实际踩过的坑是模型参数不准确情况下直接上MPC稳态会存在静差。这个问题靠卡尔曼滤波扩展状态基本能兜住模型不准导致的偏差会被看成一种“等效扰动”Q_dist估计到一定值后MPC会自动补上对应的功率偏移。如果发现还是差我会把室内温度测量和Q_dist估计一起拿去做个简单的最小二乘回归重新更新一下UA和C_room的数值。说到底MPC的好用程度永远逃不开模型和状态估计准不准这两件事。最后说一个我最近常用的习惯仿真做差不多了我会刻意把传感器噪声方差调大、把Q_true的负荷突变调猛先看看这套MPC加卡尔曼滤波在“恶劣条件”下能不能撑住。如果这种苛刻工况下还能把温度控制在0.5℃以内那真正装到现场时就比较放心了。你有兴趣的话可以在此基础上试着把Q_true改成随机游走型扰动或者加上一个阶跃负荷后马上又恢复看看卡尔曼滤波的跟踪会不会出现滞回。这个拓展实验做下来你对状态估计和预测控制的理解会比只看书本里的公式深得多。