王乘熙， 李泽*， 崔国增， 许王瑶

(1.苏州科技大学电子与信息工程学院， 苏州 215009； 2.苏州科技大学苏州智慧城市研究院， 苏州 215009)

中国每年建筑能耗占能源总消耗量的30%， 多区域暖通空调(heating,ventilation and air conditioning, HVAC)又占建筑能耗的55%。近年来虽然HVAC空调的工艺设计方法和能耗部件结构已得到了较大的改进与优化，但是大多仍采用较为传统的控制算法，容易造成室内人员体感舒适度与空调日常运行能耗无法兼顾的问题，至今依旧是建筑节能领域一项较大的挑战[1-3]。

空调房间是整个HVAC系统最重要的终端工作区域，一般可针对房间内部自身存在的热特性进行建模。现有的模型大多是单区域模型，即一个空气节点代表整个房间空气状态，与实际运行工况差别较大。针对这种情况，文献[4]采用计算流体力学(computational fluid dynamics,CFD)方法，研究室内空气流动和温湿度分布情况，将一个房间按环境差异特性分为3个区域，建立了多区域空调房间多输入多输出温湿度动态模型。

模型预测控制(model predictive control, MPC)为先进控制中的一种典型算法，因其能处理各类约束，控制器设计灵活多变等优势，在中外空调节能控制领域已有较为丰硕的成果[5-8]。文献[6]建立了一种直膨式空调机组的降阶模型，通过开环控制对能耗进行优化得到系统稳态，并采用MPC控制空气流量、风机转速以跟踪参考值。文献[7]依据建筑热特性建立了HVAC系统一阶电容-电阻网络模型，采用MPC算法并考虑天气预测扰动因素，在某地办公楼进行了实地试验。文献[8]对多个空调房间用于调节风量的空气处理单元进行建模，采用考虑通风和节能的多目标优化MPC策略以控制供气流量和温度。

然而，上述研究都只是基础MPC算法的应用，无法对外部扰动进行处理。而空调动态模型建立过程中存在的参数变化，运行工况中室内人员的走动，室外天气状况都会给系统带来不确定性影响。文献[9]设计了一种经典的鲁棒控制器，对空调的送风风机进行控制，仿真验证了外部加性扰动能较好得到处理且闭环系统渐进稳定。文献[10]建立了空调机组的传递函数模型，基于线性矩阵不等式理论，设计了H∞鲁棒控制器，对接入电网的空调机组电源进行调频控制，有效减少电源负荷。文献[11]用鲁棒不变集描述了变风量空调系统不确定性，通过设计离线椭球不变集算法，控制冷水机组中变频水泵，最后仿真验证其鲁棒性。将鲁棒控制与MPC相结合，采用改进鲁棒模型预测控制(robust model predictive control, RMPC)策略，能够明确处理系统中存在的扰动的特性，因此已有一些RMPC 算法用于建筑物内气候控制。文献[12]提出了一种基于管式 RMPC方法的HVAC空调室内环境控制策略，并引入状态观测器，以优化能源使用成本，并且能有效处理建筑热负荷变化带来的非线性扰动。文献[13]通过对建筑模型不确定性参数进行在线估计，将加性扰动与建筑模型相结合，设计了具有自适应功能的RMPC控制器，最后在某测试大楼实地进行仿真得出了该算法能使大楼内部保持较佳的热舒适性性能。然而，上述文献所采用的RMPC策略大多只是用于单区域的空调控制，而多区域空调因其各区域之间复杂耦合性，很适合用多胞模型刻画系统的不确定性，并运用RMPC算法解决[14-15]。

因此，现考虑到多区域空调房间建模时存在的不确定性扰动，基于RMPC算法设计一种多区域空调房间温湿度跟踪控制策略。采用计算流体力学方法建立的传统空调房间状态空间模型，经过处理拓展为多胞不确定性模型。为了兼顾温湿度跟踪控制效果以及空调运行工况下的稳定性约束条件，最终设计一个无穷时域带终端约束的鲁棒目标函数。利用线性矩阵不等式(linear matrix inequality, LMI)将性能指标函数最终转换为求解状态反馈控制率的问题，最终仿真验证所提出算法的可行性与有效性。

1 系统描述

1.1 多区域HVAC空调房间模型

由于单区域的HVAC空调房间系统存在一定局限性，没有考虑各区域之间的耦合，控制效果以及人体舒适性较差。采用计算流体力学方法，建立多区域的HVAC空调房间系统热力学模型，建模的具体过程如下。

假设1：房间可考虑分成3个典型区域，即空气供应，工作和回风区。每个工作区域完全融合且用一个状态表示，如图1和图2所示。若考虑更多区域，则会造成模型维度过高的问题。

假设2：每个内外墙的表面温度都由一个集总值来表示。外墙的传热只影响工作区的空气且在空气供应区和空气返回区之间没有热传递和质量传递。

假设3：在空气供应区，灯泡是唯一热源，在瞬态响应模拟过程中表面温度保持恒定。且空气供应区不存在水分来源。

假设4：在工作区，热源是电气设备和室内的人，水分来源是来自皮肤和室内人体呼吸的蒸发水。在动态响应模拟过程中，假设所有热源的表面温度恒定。

假设5：回风区没有热源和湿源。

假设6：房间内墙壁和物体之间的传热可以忽略不计。在动态响应模拟期间，墙或热源与相邻空气之间的对流传热系数认为是恒定的。

图1 多区域HVAC空调房间Fig.1 Multi-zone HVAC air conditioned room

基于上述假设，多区域HVAC空调房间模型依据能量和质量守恒定律，得式(1)～式(13)。

图2 多区域HVAC空调工作原理Fig.2 Working principle of multi-zone HVAC air conditioning

对于空气供应区[4]：

(1)

(2)

(3)

式中：ρo和ρnl分别为空气和内墙的密度；Vo,g、ho,g、To,g、Wo,g分别为供气区的体积、焓、温度、湿度；Vnl,g、cnl,g、Onl,g、Znl,g、Tnl,g分别为供气区与内墙之间空气夹层的体积、质量比热容、质量流动速率、区域面积、温度；Oo,i、ho,i、To,i、Wo,i分别为空气供应区入口的质量流动速率、焓、温度、湿度。

对于工作区[4]：

(4)

(5)

(6)

(7)

式中：ρns和λns、ζns分别为外墙的密度和导热系数；Vo,w、ho,w、To,w、Wo,w分别为工作区的体积、焓、温度、湿度；Vnl,w、cnl,w、Onl,w、Znl,w、Tnl,w分别为工作区与内墙之间空气夹层的体积、质量比热容、质量流动速率、区域面积、温度；Vns,w、cns,w、Ons,w、Zns,w、Tns,w分别为工作区与外墙之间空气夹层的体积、质量比热容、质量流动速率、平面面积、温度。

对于回风区[4]：

(8)

(9)

(10)

式中：Vo,b、ho,b、Oo,b、To,b、Wo,b分别为回风区的体积、焓、质量流动速率、温度、湿度；Vnl,b、cnl,b、Onl,b、Znl,b、Tnl,b分别为回风区与内墙之间空气夹层的体积、质量比热容、质量流动速率、区域面积、温度；Δqo,g、Δqo,w和Δqo,b为对应3个区域的空气热增长率[4]。

Δqo,g=αnl,gZnl,g(Tnl,g-To,g)+αog-owZog-ow(To,w-To,g)+αsa,gZsa,g(Tsa,g-To,g)

(11)

Δqa,r=αnl,bZnl,b(Tnl,b-To,b)+αow-obZow-ob(To,w-To,b)

(12)

Δqo,w=αnl,wZnl,w(Tnl,w-To,w)+αns,wZns,w(Tns,w-

To,w)+αog-owZog-ow(To,g-To,w)+

αow-obZow-ob(To,b-To,w)+

(13)

式中：αnl,g、αnl,w、αnl,b分别为内墙面与供气区、内墙面与工作区、内墙面与回风区的对流换热系数；αsa,g、αsa,w分别为室内热源与供气区、室内热源与工作区的对流换热系数；αog-ow、Zog-ow为供气区与工作区之间空气夹层的对流换热系数、区域面积；αow-ob、Zow-ob为回风区与工作区之间空气夹层的对流换热系数、区域面积；Tsa,g、Tsa,w为供气区热源；工作区热源的温度；i为工作区域内第i个热源。

呼出气体的焓hout与呼出气体的温度Tout、湿度Wout比率有关。在典型室内温度25 ℃下，温湿度比率计算公式[4]为

Tout=32.6+0.066To,w

(14)

Wout=0.029 33+0.2Wo,w

(15)

(16)

(17)

式中：E、M和F分别为人体的代谢率、质量和高度。

工作区的空气湿度增长率ΔeWo,w、外墙表面温度Tns,e以及空气的焓ho[4]计算公式为

ΔeWo,w=Gres(Wout-Wo,w)

(18)

Tns,e=To,out+εS/αns,e

(19)

ho=1 005To+2.5×106Wo

(20)

式中：αns,e为外墙的对流换热系数；S为太阳辐射强度；ε为太阳辐射吸收系数；To,out为外界温度。

1.2 多区域HVAC空调房间多胞体模型

为了运用MPC算法进行室内温湿度控制，必须将式(1)～式(10)转换为状态空间模型，可以忽略高阶项，随后新的状态变量变为稳态点(x0)加上状态变量变化值(δx)，于是此时定义状态变量δx=δTo,g、δWo,g、δTnl,g、δTo,w、δWo,w、δTnl,w、δTns,w、δTo,b、δWo,b、δTnl,b,定义输入变量δu=δTo,i、δWo,i、δOo,i、δTo,out、δS。式(1)～式(10)即可写为状态空间形式，即

(21)

系统矩阵Ac、Bc中的对流换热系数是影响模型精度的关键参数，受供气区、工作区、回风区室内人员行为变化，冷热源、外部环境影响[16]，式(11)～式(13)给出的空气热增长率模型已对其含义和单位进行了详细描述。因此通过刻画对流换热系数的上下界范围，原系统可拓展为顶点数为23=8的凸多胞体模型，在各个顶点处的不确定性可表示为下列变化参数所有可能的组合：

αnl,g,w,b∈[(αnl,g,αnl,w,,αnl,b)min(αnl,g,αnl,w,,αnl,b)max]

(22)

αow,ob,og∈[(αow-ob,αog-ow)min(αow-ob,αog-ow)max]

(23)

αsa,g,w∈[(αsa,g,αsa,w)min(αsa,g,αsa,w)max]

(24)

式(22)表示各区域与房间内墙面之间的对流换热系数αnl,g,w,b的范围；式(23)则表示供气区与工作区，回风区与工作区之间空气夹层的对流换热系数αow,ob,og范围；式(24)则表示房间内的室内热源与供气区、工作区的对流换热系数αsa,g,w范围。此时控制对象就由原来线性时不变的系统转化为工作在多个多胞顶点的，且每个顶点处可视为一个线性参变模型的系统。多区域HVAC空调房间多胞不确定性模型为

(25)

式(25)中：Co为凸包。状态矩阵以及输入矩阵[AzsBzs](其中s=1,2,…,8)此时包含于一个多胞不确定集ηz，工作在多胞体的不同顶点处，且βl为非负常数。

以采样时间为步长，离散后的控制模型最终表示为

(26)

2 鲁棒模型预测反馈控制策略

2.1 基于LMI的min-max RMPC策略

对于上文建立的多区域HVAC空调房间多胞体模型，提出了一种基础的BLMI-RMPC策略，该无限时域的鲁棒模型预测控制问题可表述为在多胞不确定集φz上求解最小化目标函数最坏情况的问题，系统整体结构如图3所示。最终要实现的目标是通过控制空调出风口冷空气流量变化也即压缩机工作频率变化使各区域的温湿度尽可能快地跟踪上参考值，且要保证误差较小。同时还需满足压缩机转速范围的输入以及室内温湿度范围的输出约束条件。该最小化最坏情况的min-max问题可由式(27)表示：

图3 整体系统框图Fig.3 Overall system block diagram

(27)

式(27)中：Qx和Ru分别为状态权重和控制权重。

式(27)为一个无穷时域的优化问题，并不能直接求解，因此需推导得出目标函数的上界，将之变为一个单纯控制变量数目有限的极小化问题。定义一个二次函数[17]:

V(δx)=δxTPδx，P>0,P=PT

(28)

假设V(δx)在采样时间k，对任意[Ak+t,Bk+t]∈φz,t≥0满足以下鲁棒性约束：

V[δx(k+t+1|k)]-V[δx(k+t|k)]≤

-[δxT(k+t|k)Qxδx(k+t|k)+

δuT(k+t|k)Ruδu(k+t|k)]

(29)

当无穷时刻控制系统状态无限趋于零即满足式(30)，系统才符合闭环渐进稳定的要求。

(30)

对式(29)从t=0开始做累加到∞得

V[δx(∞|k)]+…+V[δx(k+2|k)]-

V[δx(k+1|k)]+V[δx(k+1|k)]-

V[δx(k|k)]≤-[δxT(k|k)Qxδx(k|k)+

δuT(k|k)Ruδu(k|k)]-[δxT(k+1|k)×

Qxδx(k+1|k)+δuT(k+1|k)×

Ruδu(k+1|k)]-…-δxT(∞|k)×

Qxδx(∞|k)+δuT(∞|k)Ruδu(∞|k)

(31)

由式(30)得到的结果可知在无穷时刻式(31)不等式左边第一项趋近于0，此时式(31)可化简为

-V[δx(k|k)]≤-J∞(k)

(32)

此时最坏情况上界即可设为

maxJ∞(k)≤V[δx(k|k)]=τ

(33)

记P=τθ-1和ψ=Hθ-1，对应于使用二次规划的常规MPC算法仅把当前控制量δu(k)作用于系统。BLMI-RMPC算法以最小化上界τ为目标，在当前采样时刻求解LMI不等式得出唯一的状态反馈控制律，作用空调系统实现反馈控制，然后在下一时刻重新求解，计算公式为

(34)

可将式(27)表示的min-max问题转换为下列LMI不等式：

(35)

式(35)将每一时刻系统状态保持在椭圆不变集θ内，其中*代表同一矩阵中其对角元素的转置。

(36)

式(36)满足了强制V递减的条件且找出了目标函数的上界为τ。

输入和输出约束分别为

(37)

(38)

因此，BLMI-RMPC算法最终转化为

(39)

2.2 基于改进LMI的min-max RMPC策略

然而，2.1节设计的鲁棒模型预测反馈控制器依然有些缺陷，系统状态从初始时刻到终端时刻始终保持在控制不变集θ内，该集合也是终端不变集，因此无法通过缩小终端集来改善控制性能。

在多区域HVAC空调房间温湿度调节控制中，可能因终端集过大，导致在HVAC送风作用结束后3个区域环境温湿度区分不明显，人体体感不舒适。此外，若是考虑系统矩阵中存在更复杂的参数不确定性，多胞集顶点数目将成倍数增长，此时，基础RMPC算法因控制不变集的局限性难以胜任复杂的LMI求解，甚至造成空调系统运行故障。

文献[18]提出了一种相比2.1节基础BLMI-RMPC 算法具有显著优势的切换ILMI-RMPC策略。该算法在一开始将N个自由控制变量作用于系统，而后采用单一状态反馈控制率[18]即

(40)

将式(27)的目标函数分为两部分：

(41)

(42)

对于任意t≥0满足式(27)的约束条件。

此时定义新的时变李雅普诺夫函数：

V(t,k)=δx(k+t)TP(t,k)δx(k+t),t≥N

(43)

假设函数对任意[Ak+t,Bk+t]∈φz,t≥0满足下列鲁棒不等式约束：

(44)

因此，同样对式(44)从t=N做累加到∞得

V(∞,k)+V(N+2,k)-V(N+1,k)+V(N+1,k)-

(45)

不等式(45)右边项即为-J2(k)，再对式(45)进行化简得

-V(N,k)≤-J2(k)

(46)

此时式(42)表示的无穷时域最优问题即等同于最小化V(N,k)的上界：

maxJ2(k)≤V(N,k)=δx(k+N)TP(t,k)δx(k+N)

(47)

结合式(41)目标函数[式(27)]也完成了无穷时域到有限时域的转变：

(48)

可把不等式(48)右边最后一项看作常规MPC算法中的终端约束，因此最终改进的ILMI-RMPC策略采用下列性能指标函数：

(49)

通过式(26)给出的状态空间多胞模型，对未来N步的状态预测矩阵可表示为

(50)

整理成分块的形式:

(51)

式(51)中：

由式(50)获得。

这时，目标函数式(49)改写为

(52)

当且仅当存在L个对称正定矩阵pl=(1,2,…，L)满足下列条件[式(53)]，式(40)和式(44)才能成立。

[Azl+Bzlψ(k)]TPt[Azl+Bzlψ(k)]-Pl+

Qx+ψT(k)Ruψ(k)≤0，l,t=1,2,…,8

(53)

定义新的上界为

(54)

(55)

最终多区域空调房间温湿度优化问题的目标函数可改为

(56)

l,j=1,2,…,8

(57)

(58)

(59)

关于上界τ1的式(54)可表示成下列LMI：

l=1,2,…,8

(60)

同理，也可用新的多胞来表述终端增广矩阵

即

(61)

(62)

关于上界τ2的式(55)表示为下列LMI：

l=1,2,…,8

(63)

同样，输入输出约束可表示为

(64)

(65)

最终，完整的优化问题变为

(66)

s.t.式(37),式(40),式(43),式(44),式(45)

此时ILMI-RMPC算法的状态反馈控制率变为

(67)

3 仿真分析

3.1 不同RMPC策略性能比较

为了测试所提出的RMPC算法的鲁棒性能，对2.2节提出的改进ILMI-RMPC算法与2.1节基础的BLMI-RMPC算法的温湿度跟踪效果与处理模型扰动的能力进行研究。算例1的仿真总时长为6 h，采样间隔To,g为6 min。空调房间系统初始状态如表1所示,权重矩阵设计如表2所示。

因人体感觉较舒适的温湿度随室外环境变化是不相同的，且考虑到空调房间动态热行为特性，供气区、工作区和回风区温湿度呈递增差异[19]，所以在仿真的前3 h，输出参考值设定为[26.4,18.1,29.5,18.3,30.0,18.3]，而后3 h设置为[23.9,17.0,27.7,17.2,28.3,17.2]。并且选取将状态反馈控制率δu应用于系统前的自由控制变量数目N为4，且多胞不变集顶点数目L此时为8，采用yalmlp 工具包对LMI不等式进行求解，并在MATLAB平台完成仿真。

表1 算例1初始条件

表2 算例1权重矩阵设计

两种RMPC算法控制性能比较如图4所示，可明显看出通过多胞描述的模型不确定性扰动并未在空调运行过程中造成较大波动。在BLMI-RMPC算法作用下，供气区、工作区、回风区的温湿度大约36 min才可到达设定的前一组参考值。而后3 h设定值变化，3个区域温湿度依旧能快速跟上参考，达到使人体感到较为舒适的温湿度。相比BLMI-RMPC 算法，ILMI-RMPC算法能在24 min左右更快跟踪上3个区域的湿度，虽然两种算法的温度跟踪效果相差无几，但是若是多胞模型顶点数进一步扩张，BLMI-RMPC算法可能就会失去其可行性。在运行周期内，两种算法都能使系统输出保持在约束范围内，始终满足状态约束条件。

3.2 变化扰动大小的RMPC策略性能比较

为了进一步验证所提出ILMI-RMPC的优越性，考虑分析在模型不确定性扰动进一步复杂的情况下，多区域HVAC空调房间温湿度的鲁棒跟踪性能。随着室外温度，光照强度的变化，式(7)描述的系统模型中的房间外部墙面的导热系数λns也会受到影响，其变化范围[16]为

(68)

图4 各区域温湿度控制仿真效果Fig.4 Simulation results of temperature and humidity control in each zone

此时，在1.2节多胞模型(26)的基础上顶点数拓展为16，鲁棒控制器性能要求更高。该算例初始条件如表3所示，权重矩阵设计如表4所示，为更好地分析控制系统对不同参考量的跟踪效果，新的输出参考前3 h设为[21.8,15.7,26.5,15.9,26.7,15.9]，后3 h设为[24.0,17.2,27.8,17.4,28.2,17.4]。

表3 算例2初始条件

表4 算例2权重矩阵设计

由图5可明显看到，此时BLMI-RMPC算法已经无法胜任控制空调房间湿度的任务。在运行期间，供气区、工作区、回风区湿度在6 min时陡增至大于100 g/(kg dry air)，远远超出运行实际情况。空调出风口容易形成冷凝水并且排出影响建筑立面，造成环境污染和资源浪费，还有可能导致HVAC空调故障[21]。

相比基础BLMI-RMPC算法,本文所提出的改进ILMI-RMPC算法就显示了其优越性，如图5(b)～图5(f)所示，3区域湿度虽在初始时间有略微的超调，但是很快就得到修正，且能保持在约束范围内。

3.3 改进RMPC策略与传统MPC策略性能比较

为了更好地体现ILMI-RMPC算法处理模型扰动的能力，与传统的MPC算法[20]进行仿真对比，在后文传统MPC算法统称为B-MPC算法。算例3初始条件和输出参考值与上一节设定相同，权重矩阵如表5所示。B-MPC算法面对模型中存在的扰动，具有天然劣势，无法将系数矩阵参数的变化囊括到控制模型中，只能在每一步使用二次规划求解最优控制输入的时候，由外界施加一个扰动作用于空调的工作状态如温度或湿度，该小节算例采用的是由均值为0，标准差为0.6的高斯分布生成的随机干扰。

仿真结果如图6所示，B-MPC算法虽然有着求解速度较快的优点，但是一旦受到扰动影响，将会缺乏控制的稳定性。当供气区、工作区、回风区达到第一次设定的参考温度后，并不会趋于平滑，而是围绕±0.4 ℃波动，而这导致的后果就是HVAC空调运行机组部分输送冷量的阀门频繁动作，以及制冷机能耗严重。

ILMI-RMPC方法将扰动考虑在控制模型建立中无疑能更好地反映HVAC空调运行实际工况。仿真结果显示，该算法跟踪到参考目标的控制时间虽然要比B-MPC算法长，由图6(b)～图6(f)可看出，3个区域的湿度在控制初期有个约为1.5 g/(kg dry air)的超调量，但是符合约束条件，并且不会持续振荡，且一旦跟踪上参考，3个区域的温、湿度便趋于稳定，模型中扰动的影响已得到较好处理。

由图5和图6可以粗略看出存在扰动情况下，3种算法中各区域温湿度跟踪效果最好的为本文提出的ILMI-RMPC算法，为更直观地表现其优越性，表6给出了各算法的平均跟踪误差(mean tracking error, MTE)[22]。

表6 算例2，算例3中3种算法平均跟踪误差

图6 各区域温湿度控制仿真效果Fig.6 Simulation results of temperature and humidity control in each zone

由表6可见，B-MPC算法3区域湿度平均跟踪误差较小，但是温度误差却达到了±0.4 ℃，BLMI-RMPC算法各区域湿度的平均跟踪误差显然太大，而ILMI-RMPC算法的温湿度误差都很小，综上，ILMI-RMPC算法跟踪效果最好。

4 结论

针对多区域HVAC空调房间温湿度控制系统建模中存在不确定性扰动的问题，提出了一种基于LMI的改进min-max鲁棒模型预测控制策略用于空调房间温湿度跟踪控制。在MATLAB仿真平台上，利用yalmlp工具包对所提出算法中线性矩阵不等式部分进行求解，得到如下结论。

(1)采用计算流体力学方法建立的多区域空调房间动态模型，能够通过传统MPC策略有效运用于温湿度跟踪控制中。

(2)考虑多区域空调房间建模中存在的不确定性扰动影响，BLMI-RMPC算法尽管已经通过多胞化处理，实现了对空调各区域温湿度的控制，几乎消除了扰动影响。然而一旦模型不确定性范围进一步扩大，该算法则无法兼顾状态约束的合理性以及BLMI-RMPC算法的控制性能。

(3)最终提出的ILMI-RMPC算法通过在单一状态反馈控制率前加入N个自由控制变量来改善鲁棒性能。B-MPC算法无法处理扰动因此在仿真实验中围绕参考剧烈波动，而ILMI-RMPC算法相较BLMI-RMPC算法能够应对扰动更加复杂的情况，且温湿度跟踪更加稳定平滑，保证环境舒适度的同时降低能耗。

在下一步工作中，因采用机理方法建立的模型其外部扰动确定存在局限性，拟考虑在微电网中建立HVAC系统能耗负荷模型，在增强模型可靠性的同时还能囊括实时电价的因素，进一步节约电力资源。