应用IMOABC空调房间CO2浓度二自由度内模分数阶PI控制性能的研究

杨 睿 李绍勇 魏明松 王少波

（1.兰州理工大学土木工程学院 兰州 730050；2.广州市恒盛建设工程有限公司 广州 510080）

0 引言

伴随世界经济快速增长，人们的生活质量得到了明显改善。同时，对建筑环境的舒适度、空气品质的关注度也与日俱增，使得中央空调系统在各类建筑环境中得到了广泛应用。由于定风量空调系统（Constant Air Volume,CAV）具备结构简单、初始投资低、易操作和送风量稳定等特点，在高大空间和室内人员密度大的建筑场所，例如大型会议室、商场和办公大厅等，广泛应用[1]。对于CAV系统中的关键设备CAVAHU而言，新风量qoa、回风量qra以及排风量qea经空调工艺设计后，往往保持不变[2]。当CAVAHU实际运行，常出现空调房间内的额定人员超员或不足的状况，导致qoa不够或过量，产生室内CO2浓度超标（Cn＞Cn,set）或处理新风所消耗的冷、热负荷增加的问题，使得室内空气品质下降或增大CAVAHU运行能耗。

室内CO2浓度变化呈现非线性、多干扰等特性，采用传统控制方式很难得到一个满意的控制效果。因此，室内CO2浓度的控制问题，是当前暖通空调自动化领域的一个有研究价值的课题。Kusuda[3]针对传统的机械通风系统，首次提出基于室内CO2平均浓度的需求控制通风策略，将房间CO2平均浓度作为调节新风量大小的控制指标。Igor 等[4]在通风系统上设计了室内CO2浓度内模控制系统，对输送到室内的新风量qoa进行调节。结果显示，qoa能够根据室内CO2浓度的变化快速做出响应。文献[5]利用直接反馈线性化设计了基于室内CO2浓度变化的机械通风控制系统。结果表明，该控制系统对于室内CO2浓度控制的速度和稳定性等方面都优于传统PID 控制系统。

分数阶PID（Fractional Order PID Controller,FOPID）控制器比整数阶PID 控制器[6]更具有设计自由度广、自适应好和强鲁棒性等特性[7]，在许多工业领域得到了应用，呈现出了更好的控制性能[8-10]。将内模控制与FOPID 控制相结合能使得控制系统获得更好的设定值跟踪和抗干扰等性能[11,12]。此外，控制器的参数整定对于控制系统的稳定性和快速性具有显著意义。目前，许多研究工作都集中在控制器参数整定的单目标优化问题[9-13,14]。但在实际应用中，控制系统所输出的被控参数的多个性能指标之间往往是相互制约的[15]。因此，解决控制器参数整定的多目标优化问题，显得十分迫切。文献[15]提出了一种基于改进多目标极值优化算法的自动稳压器系统FOPID 控制器设计方法。通过积分绝对误差，绝对稳态误差和调节时间三个性能指标，建立多目标优化函数来整定FOPID 控制器参数的最优值。结果表明，该方法使自动稳压器系统的控制精度和鲁棒性方面有明显提升。针对旋转式倒立摆系统在稳定阶段的控制问题，文献[16]选择摆角绝对积分时间误差和转角绝对积分时间误差两个目标函数，建立多目标优化问题，并设计了一种基于改进多目标粒子群算法来实施FOPID 控制器参数的整定。通过试验验证了该算法的有效性，使得倒立摆系统具有良好的稳定性和动态品质。

综上分析，空调房间CO2浓度与输配到每个空调房间的新风量密切相关。基于保持CAVAHU 输出的总送风量qsa恒定，同时动态调节新风量qoa和新风比m，确保Cn=Cn,set的目标，本文提出一种室内CO2浓度二自由度内模分数阶PI 控制策略，发挥控制器参数少，系统调节灵活，快速跟踪和鲁棒性强等优势。此外，对于室内CO2浓度二自由度内模分数阶PI 控制器参数的整定问题，依据人工蜂群算法，对雇佣蜂阶段引入自适应惯性权重，同时对观察蜂阶段施加精英组策略，分别进行非线性递减和柯西变异的演变，并结合观察蜂搜索特性，将最小粒子角度[17]引入外部档案集，获取相应的Pareto 解集，设计得到IMOABCA。考虑调节时间tc、稳态误差绝对值ESS和绝对积分时间误差（Integrated Time Absolute Error,ITAE），选择min（ITAE,ESS,tc）作为IMOABCA 目标函数，设计室内CO2浓度二自由度内模分数阶PI 控制器参数整定算法。仿真结果表明：该参数整定算法可获取3 个控制器参数和3 个目标函数的Pareto 最优解。从而提升室内CO2浓度调节品质，确保Cn=Cn,set和满足室内空气品质需求。同时，二自由度内模分数阶PI 控制器作用下的CAVAHU，其输出的新风量qoa是动态变化的，可适应空调房间内的额定人员超员或不足的实际情况。

1 CAVAHU 作用下的空调房间CO2浓度控制

对于新风+一次回风的空调方式，CAVAHU 输送的总送风量qsa包括新风量qoa和一次回风量qra，即qsa=qoa+qra，新风比m=qoa/qsa。同时为了维持室内的微正压，空调工艺的设计使得新风量等于排风量qea，即qoa=qea，且qsa＞qra。由上述可知，CAVAHU输配的qoa对空调区域的CO2浓度Cn影响显著。所以综合考虑CAVAHU 输出的qsa恒定和空调房间内的额定人员超员或不足的实际状况，本文提出一种室内CO2浓度二自由度内模分数阶PI 控制策略，能够动态调节新风量qoa，实现Cn=Cn,set的目的。相应的空调工艺测控流程图，如图1所示。

图1 CAVAHU 运行工艺测控流程图Fig.1 Flowchart of measurement and control for CAVAHU operation process

当CAVAHU 运行空调冬/夏季工况时，设置在一次回风管道上CT 将反映n个空调房间平均CO2浓度测量值信号Cn上传到CC，与经过CF 的室内CO2浓度设定值Cn,set进行比较，求偏差e=Cn,set-Cn，对其进行内模分数阶PI 运算后，输出3 路控制信号p1、p2和p3，分别到新风电动调节阀、一次回风电动调节阀和排风电动调节阀。这里，p1=p3，p2=pmax+pmin-p1。p1、p2、p3∈[pmin，pmax]，pmin，pmax-CC 输出的控制信号的下限、上限。动态改变输送到空调房间的新风量qoa，经过通风置换，实现Cn=Cn,set的目的，适应空调房间内的额定人员超员或不足的实际情况。同时CAVAHU 以变新风比（m＜100%）运行，输出的总送风量qsa是恒定的；且维持室内的微正压。这是由于当e＜0（＞0）时，接受控制指令p1和p3，新风阀开度v1和排风阀开度v3是同向增大（减小），即v1=v3=a%（某个开度值），qoa=qea；而接受控制指令p2，回风阀开度v2（100%-a%）是反向减小（增大），即v1+v2=100%，qoa+qra=qsa。

这样，该室内CO2浓度二自由度内模分数阶PI 控制策略，可以实现CAVAHU 对室内CO2浓度的实时调整，既动态调节新风量qoa，实现Cn=Cn,set的目的；又使得CAVAHU 以输出的qsa是恒定的和变新风比（m＜100%）运行。

2 室内CO2浓度二自由度内模分数阶PI 控制系统

2.1 空调房间CO2 浓度的数学模型

分析文献[5]，可得空调房间CO2浓度质量平衡方程如下：

式中，qoa为新风量，m3/s；Cw和Cn为室外和室内CO2浓度，%；mn为室内CO2发生量，m3/s；V为房间体积，m3。

由于空调CO2浓度变化是非线性的，利用式（1）对空调房间CO2浓度模型在平衡点（稳定点）建立线性过程，得到相应的传递函数为[4]：

式中，a=-(qoa*/V)，b=(Cw-Cn*)/V。

这里，qoa*、Cn*-qoa、Cn的平衡点。两者的关系如下：

式中，m*-qoa*对应的新风比。

2.2 空调房间CO2 浓度测量变送器传递函数

室内CO2浓度测量变送器的作用是实时测量空调房间内的CO2浓度Cn大小，输送标准信号到CC控制器与Cn,set进行比较，其传递函数如下式[18]：

2.3 新风量、排风量与一次回风量电动调节阀传递函数

由于新风量qoa可直接影响房间CO2浓度Cn，所以设置新风量电动调节阀，接受输出的控制指令p1（0～10mA·DC 或4～20mA·DC），通过改变其开度v1来控制qoa大小。对于排风量qea和一次回风量qra，由第1 节可知qea=qoa，qoa+qra=qsa，且p1=p3，p2=pmax+pmin-p1。所以，三个风阀本文均选择直线流量特性的电动调节阀，其可看作放大环节处理[19]：

2.4 室内CO2浓度二自由度内模分数阶PI 控制器传递函数

从上文可知，新风量电动调节阀与空调房间所构成的广义被控对象为一阶惯性数学模型：

式中，K=10b/-a；T=b/-a。

分析文献[11]可以得出，内模分数阶PI 控制器的分数特性是通过施加的分数阶滤波器体现的，其传递函数如下：

式中，T1、α为分数阶滤波器的时间常数和阶次的非整数项（0＜α＜1）。

从上式可以看出，内模分数阶PI 控制器是由一个整数阶PI 控制器与一个分数阶积分器串联而成，其中KP=T/KT1，KI=1/KT1。

此外，为了提升参数控制系统对设定值的跟踪特性，加快响应速度和降低超调量，引入二自由度内模控制理论[20]，构建设定值前馈滤波器，其传递函数如下：

式中，T2为一阶低通滤波器的时间常数；T1、α同上。

所以，考虑内模分数阶PI 控制器和设定值前馈滤波器的优势，本文提出的室内CO2浓度二自由度内模分数阶PI 控制器，其结构如图2所示。

图2 室内CO2 浓度二自由度内模分数阶PI 控制器结构Fig.2 Block diagram of structure for two-degree of freedom internal model fractional order PI controller for the indoor CO2 concentration

如上所述，该控制器可输出3 路控制信号p1、p2和p3，分别作用于新风量、一次回风量和排风量电动调节阀，在保证送风量qsa不变的情况下，动态调节新风量qoa，实现Cn=Cn,set的目的。

因此，本文提出的基于IMOABCA 空调房间CO2浓度二自由度内模分数阶PI 控制系统，如图3所示。

图3 基于IMOABCA 空调房间CO2 浓度二自由度内模分数阶PI 控制系统方框图Fig.3 Block diagram of two degree of freedom internal model fractional order PI control system based on IMOABCA for the concentration of CO2 in air-conditioning room

图3为二自由度内模分数阶PI 控制器作用下的空调房间CO2浓度控制系统方框图，并采用IMOABCA 对控制器参数进行整定。考量ITAE，ESS和tc，确定min（ITAE,ESStc）作为IMOABCA的目标函数，进而对控制器的三个参数[T1,T2,α]进行整定，获取相应的Pareto 最优解：[T1*,T2*,α*]和[ITAE*,ESS*tc*]，提升调节质量，使得Cn=Cn,set，满足室内空气品质要求；且动态调节新风量qoa，减少CAVAHU 处理新风所消耗的冷、热负荷。

3 控制器参数整定

3.1 多目标人工蜂群算法改进方法

多目标优化问题是指当在解决一个问题时，需要对多个目标进行最大或最小优化，但这几个目标都是相互制约，相互矛盾的。因此，在多目标优化过程中会在多个目标问题之间寻求一个平衡值，从而得出这些目标函数的Pareto 解集[17]。

人工蜂群算法是一种简单、高效的群体仿生优化算法，具有收敛快、鲁棒性强等特点[21]。但由于采用随机选择邻居的策略，导致Akbari 等[23]所设计的MOABCA在雇佣蜂阶段的前期全局搜索能力弱，后期局部探索能力不足。因此，加入余弦自适应权重因子w，利用余弦函数在[π/2,π]呈现单调递减变化，使得雇佣蜂：（1）前期具有较强的全局搜索能力。（2）随着算法迭代次数的增加，后期也具有充足的局部探索能力。这样，雇佣蜂在全局搜索任务中的能力得以提升。相应的w及改进后的雇佣蜂搜索公式如下：

式中，it、itmax为算法的当前迭代次数和最大迭代次数；φi,j为[-1.1]之间的随机数；xi,j为当前蜜源位置；xk,j为随机蜜源位置；vi,j为新蜜源位置。

当观察蜂对较好蜜源进行挖掘时，MOABCA往往存在陷入局部最优或早熟的问题。因此，引入柯西算子Cauchy对其步长进行变异，使得观察蜂：（1）当出现陷入局部最优的趋向时，利用较长的步长跳出局部极值。（2）当在最优解附近挖掘时，利用较小的步长加速收敛。标准的柯西分布的概率密度表达式如下：

与此同时，本文对MOABCA 的外部档案（External Files,EF）的网格方法应用最小粒子角度[17]，并基于精英组策略，利用观察蜂“优中寻优”的特点对网格中的Pareto 解集进行优化，进而提升Pareto 解的质量。为了保证算法在搜索过程中的多样性，加入控制参数MR[22]进行修正。改进后的网格结构（基于两目标函数优化）和观察蜂搜索公式分别如图4和公式（12）所示：

图4 基于最小粒子角度的网格结构图Fig.4 Grid structure based on minimal particle angles

式中，xi,j为雇佣蜂发送的优质蜜源；leader为xi,j对应EF 中的最优蜜源；repk,j为EF 中随机生成的蜜源；Cauchy为柯西算子；MR为控制参数，一般取0.5[22]。

分析图4，可知观察蜂从雇佣蜂搜索的蜜源中选择较优的蜜源，进行了开采。与此同时，通过计算所选的开采蜜源与网格中的蜜源角度，选取EF中与其蜜源角度最小的蜜源，作为最优蜜源进行柯西变异，达到优化EF 中的Pareto 解集目的。

因此，基于改进的网格结构以及搜索公式，本文在MOABCA 的基础上构建出了IMOABCA。

3.2 室内CO2浓度二自由度内模分数阶PI 控制器参数整定算法

本节选择调节时间tc、稳态误差绝对值ESS和绝对积分时间误差ITAE[9]作为IMOABCA 的三目标函数，即F(x)=minJ=min(ITAE,ESS,tc)，来设计室内CO2浓度二自由度内模分数阶PI 控制器参数整定算法。基于IMOABCA 改进内容，相应运算流程如图5所示。

图5 控制器参数整定算法流程Fig.5 Parameter tuning algorithm flow of the controller

Step1：设置MALAB/Simlink 工具组态中的室内CO2浓度二自由度内模分数阶PI 控制器参数变量[T1,T2,α]，并设定T1、T2和α的上限、下限。

Step2：初始化：设置决策空间维数D，种群规模N，雇佣蜂数目Ne=观察蜂数目Ns=N，阈值limit，外部档案规模nEF，最大迭代次数itmax以及多目标函数F(x)=minJ。

Step3：输入单位闭环负反馈室内CO2浓度二自由度内模分数阶PI 控制系统的反馈信号e(t)，基于min(ITAE,ESS,tc)，计算每个初始蜜源的适应度值F，完成网格非支配排序。

Step4：进入雇佣蜂阶段，根据公式（10）生成新蜜源vi,j，计算目标函数值F，并与旧蜜源进行支配关系判断：若新蜜源支配旧蜜源，则替代旧蜜源，成为非劣解，雇佣蜂阶段结束。非劣解被加入EF 中进行网格非支配排列，并删除网格中被支配蜜源。

Step5：利用轮盘赌方式，在雇佣蜂阶段获取的蜜源F中，为观察蜂选择开采对象。通过计算开采对象与EF 中的蜜源之间的角度，选取角度最小的蜜源作为最优蜜源。将开采对象和最优蜜源同时发送给观察蜂。

Step6：进入观察蜂阶段，根据公式（12）生成新蜜源vi,j，计算目标函数值F，并与收到的开采对象、最优蜜源进行支配关系判断：若新蜜源支配开采对象或最优蜜源，则替代它，成为非劣解，观察蜂阶段结束。非劣解被加入EF 中进行网格非支配排列，并删除被支配蜜源。

Step7：进入侦察蜂阶段，若某蜜源的开采次数超过了阈值limit，则随机生成新蜜源替代旧蜜源，并计算目标函数值F。反之，则继续保留旧蜜源。

Step8：若满足it=itmax，输出控制器参数最优值[T1*,T2*,α*]，ITAE、tc和ESS的Pareto 最优解minJ*以及相应的闭环负反馈控制系统动态响应曲线。否则，返回Step4。

需要说明的是，借助MATLAB 软件，将上述的控制器参数整定流程进行编程，保存为一个独立程序，命名为IMOABCA.m。该IMOABCA 是与后面章节中的控制系统组态Simulink 模型同步运行的。

4 数值模拟

4.1 室内CO2浓度二自由度内模分数阶PI 控制系统的组态与仿真

对于兰州市一个大型办公区域，其容积V=1200m3，额定人员90 人，人员的活动强度为轻度，人均CO2发生率为0.0173m3/h。根据GB/T 18883-2022 标准，室内CO2浓度不应超过0.1%[25]，故本文选择室内CO2浓度设定值Cn,set=0.1%，相应的CAVAHU 总送风量qsa=4.1m3/s。对该CAVAHU配置本文设计的室内CO2浓度二自由度内模分数阶PI控制系统。基于图3和MATLAB中的Simulink工具，对该控制系统进行组态，相应的Simulink模型如图6所示。

图6 基于IMOABCA 的室内CO2 浓度二自由度内模分数阶PI 控制系统组态图Fig.6 Configuration diagram of two-degree of freedom internal model fractional order PI control system for the indoor CO2 concentration based on IMOABCA

根据空调工艺要求[25]，Cn,set=0.1%，假定室内CO2浓度初始值Cn,0=0.11%，IMOABCA.m 中的相关参数设置如下：N=100，nEF=3，控制器参数整定的数目D=3，itmax=200，limit=300，w∈[0,1]，T1∈[90,110]，T2∈[10,12]，α∈[0.4,0.6]，仿真时间设为500s。

在MATLAB 中的Command Window 界面和Simulink 环境，同步运行IMOABCA.m 和图6中的Simulink 模型，可获得室内CO2浓度的动态响应Cn(t)，min(ITAE,ESS,tc)Pareto 前沿变化和3 组控制器参数Pareto 最优值[T1*,T2*,α*]，分别如图7，图8和表1所示。

图7 三组控制器参数Pareto 最优解下的室内CO2 浓度动态响应Fig.7 Dynamic responses of indoor CO2 concentration under three Pareto optimal solutions of parameters of the controller

图8 min（ITAE，ESS，tc）Pareto 前沿变化Fig.8 Pareto front variation for min(ITAE,ESS,tc)

分析图7可知，基于IMOABCA 得到的控制器三组最优参数均能使得Cn=Cn,set，说明IMOABCA 对控制器的参数整定是可行的。从局部放大图可见，输出三组Cn(t)在超调量，调节时间表现出了不同的差异，也与表1中的三个目标函数值所表征的相呼应。

表1 控制器参数及其三目标函数的Pareto 解Table 1 Pareto solution of parameters of the controller and three objective functions

分析图8，容易观察到3 个目标函数值所对应的不同控制器参数的Pareto 解之间的联系，曲面上的控制器参数Pareto 解的分布位置对应了不同的Cn(t)控制特征。

分析表1，可知控制器参数的三组Pareto 最优解所对应的三个目标函数值之间是相互制约或矛盾的。例如，第1 组中的ITAE 最大，但tc和ESS最小；第3 组中的ITAE 最小，但ESS和tc最大。由此说明IMOABCA 可以使得室内CO2浓度二自由度内模分数阶PI 控制器获得侧重点不同的Cn(t)控制效果，以满足不同控制需求。

相应的新风量qoa、一次回风量qra变化和新风比m变化，如图9所示。

图9 三组Pareto 控制器参数最优解下的响应变化Fig.9 Response changes under three Pareto optimal solutions of the controller

分析图9的（a）和（b），控制器的3 组Pareto参数最优解均能够使控制器根据室内CO2浓度的变化快速做出反应，发送控制信号作用到新风量、一次回风量和排风量电动调节阀完成新风量的调节，达到Cn=Cn,set的目的。图（c）反映了控制器的3 组Pareto 参数最优解的新风比m变化过程，计算得到3 组平均新风比均约为14%，满足m≧10%[25]空调工艺要求。

4.2 与现有算法的比对

此外，为了说明本文IMOABCA 的优越性，分别将IMOABCA 和MOABCA 作用于室内CO2浓度二自由度内模分数阶PI 控制系统对控制器参数进行整定。可获得其控制器参数Pareto 最优解[T1*,T2*,α*]及其三个目标函数值如表2所示。

表2 控制器参数及其相应的三目标函数的Pareto 解Table 2 Pareto solutions of parameters of controllers and the corresponding three objective functions

相应的室内CO2浓度的动态响应Cn(t)，如图10所示。

分析表2，从min(ITAE,ESS,tc)Pareto 最优值不难看出，相比MOABCA，IMOABCA 对本文所设计的控制系统表现出了更好的收敛性和控制精度。同时，从图10可以看出，IMOABCA 作用下的控制系统具有更快的响应速度和更小的超调量，且能够很快的达到稳定状态，说明IMOABCA 在解决多目标优化问题上更具有优越性，可以使系统呈现较强的控制性能，也表明IMOABCA 的设计是可行、有效的。

图10 两算法下的室内CO2 浓度动态响应Fig.10 Dynamic responses of indoor CO2 concentration under two algorithms

5 结论

针对CAVAHU 作用下的空调房间CO2浓度的控制问题，本文提出了一种基于IMOABCA 的二自由度内模分数阶PI 控制策略。依据人工蜂群算法，对雇佣蜂阶段引入自适应惯性权重，同时对观察蜂阶段施加精英组策略，分别进行非线性递减和柯西变异的演变，并结合观察蜂搜索特性，将最小粒子角度引入外部档案集，获取相应的Pareto 解集，设计了IMOABCA，进而对室内CO2浓度二自由度内模分数阶PI 控制器的3 个参数[T1,T2,α]进行整定。结果表明：（1）该IMOABCA 是可行的，可获取控制器的3 个参数相应的最优解[T1*,T2*,α*]，从而提升室内CO2浓度的调节品质，实现Cn=Cn,set的目的。（2）该IMOABCA 可以使得控制器获得侧重点不同的室内CO2浓度控制效果，以满足不同控制需求。（3）与现有算法相比，IMOABCA对本文所设计的系统呈现出了明显的优越性，表明了IMOABCA 的设计是可行、有效的。此外，CAVAHU 是以变新风比m、保持总送风量qsa恒定的方式运行的，为CAVAHU 运行、管理提供了新方式。