CN105652665A - 一种微型燃气轮机冷热电三联供系统的协调控制方法 - Google Patents

Abstract

本发明公开了微型燃气轮机冷热电三联供系统的协调控制方法，包括离线阶段和在线阶段；离线阶段包括：步骤1-1，简化被控系统；步骤1-2，在MGT-CCHP系统处于稳定运行状态时，依次在各输入变量中加入高斯白噪声扰动，采集输入变量和输出变量的数据作为模型辨识的数据；步骤1-3，根据扰动实验的数据，利用子空间辨识的方法辨识出该三输入三输出系统的离散状态空间模型；在线阶段包括：步骤2-1，观测当前时刻k的状态向量x(k)的值；步骤2-2，据当前时刻k的状态向量x(k)、k-1时刻的输入向量u(k-1)以及辨识出的状态空间模型，推导出未来P个时刻的输出向量的表达式；步骤2-3，优化性能指标；步骤2-4，利用二次规划的方法，求解优化性能指标最小的优化问题，得到当前时刻的控制增量。

Description

一种微型燃气轮机冷热电三联供系统的协调控制方法

技术领域

本发明涉及一种工业控制方法，特别是一种基于预测控制算法的微型燃气轮机冷热电三联供系统的协调控制方法。

背景技术

基于微型燃气轮机的冷热电三联供系统(MicroGasTurbineBasedCombinedCoolingHeatingandPowersystem,以下简称为MGT-CCHP系统)是一种具有代表性的分布式能源的形式，由微型燃气轮机、供热系统、制冷系统组成。动力系统利用燃料燃烧产生的热量来发电，发电后产生的烟气余热被送至供热系统和制冷系统实现余热利用。传统的火力发电的能量转换效率在30-35％左右，而三联供系统对低能量品位的烟气余热再利用，其平均能量转换效率可高达80％，具有广阔的应用前景。

MGT-CCHP系统的控制包括微型燃气轮机、双效吸收式制冷机和集成系统三个方面，即分别研究微型燃气轮机侧、供热制冷系统侧和整体的动态特性并建模控制。燃气轮机的控制比较成熟，燃气轮机的各重要设备厂商均为各自的设备提供了相对完善的控制系统，学者们也提出使用PID控制器以及一些先进控制器，如神经网络控制、模糊控制、预测控制等来解决燃气轮机控制问题。然而对于微型燃气轮机控制系统的研究基本上是对燃气轮机的横向延伸，忽略了微型然轮机中余热利用需求、变转速运行状况以及回热器惯性等需求。吸收式制冷机控制系统主要采用PID控制器。对于双效溴化锂吸收式制冷机控制系统的研究存在以下问题：一是吸收式制冷机大多数孤立运行，缺乏与微型燃气轮机组成联供系统时的研究；二是大部分吸收式制冷机的负荷输出单一(仅供冷或供热)，调节量少(仅热源调节)，控制系统普遍为单输入单输出结构；三是吸收式制冷机中有很多热交换设备，热惯性较大。上述研究主要针对各装置独立运行，而针对各装置组成的MGT-CCHP集成系统的控制研究较少。组合后的MGT-CCHP系统是一个有多个输出负荷需求(发电、供热、制冷)的多变量系统，有较强的耦合度，因此应当从协调控制的角度来考虑。针对强耦合的特点，有学者提出了解耦PID控制的方法，即先将该多变量系统解耦成多个单变量系统，并分别对解耦后的单变量系统采用PID控制。解耦PID可以很好地解决多变量的问题。根据分析，MGT-CCHP系统是一个多输入多输出的大惯性系统，有输入量约束。而传统的PID控制方法在大惯性以及输入约束的处理方面较弱。

发明内容

发明目的：本发明所要解决的技术问题是针对传统PID控制在处理微型燃气轮机冷热电三联供系统中被控量约束以及大惯性问题的不足，提供一种微型燃气轮机冷热电三联供系统的协调控制方法。

为了解决上述技术问题，本发明公开了一种一种微型燃气轮机冷热电三联供系统的协调控制方法，包括离线阶段和在线阶段两个部分；

离线阶段包括：

步骤1-1，简化被控系统；

步骤1-2，在MGT-CCHP系统处于稳定运行状态时，依次在各输入变量中加入高斯白噪声扰动，采集输入变量和输出变量的数据作为模型辨识的数据；

步骤1-3，根据扰动实验的数据，利用子空间辨识的方法辨识出该三输入三输出系统的离散状态空间模型；

在线阶段包括：

步骤2-1，利用卡尔曼滤波器来观测当前时刻k的状态向量x(k)的值；

步骤2-2，据当前时刻k的状态向量x(k)、k-1时刻的输入向量u(k-1)以及辨识出的状态空间模型，推导出未来P个时刻的输出向量的表达式；

步骤2-3，优化性能指标；

步骤2-4，利用二次规划的方法，求解优化性能指标最小的优化问题，得到当前时刻的控制增量。

步骤1-1中，简化被控系统包括：

保持冷媒水流量、活热水流量以及相应入口温度三者恒定；保证冷剂泵始终全开；制冷机中溶液循环泵处于自动运行模式；微型燃气轮机切换到带变速负载运行；

将微型燃气轮机冷热电三联供系统方简化为一个三输入变量三输出变量的系统，三个输出变量是汽轮机转速y_n、冷媒水温度y_Tclo以及生活热水温度y_Thwo，三个输入变量分别是微型燃气轮机的燃料流量控制阀门u_mbf、再热阀门开度u_re以及双效溴化锂吸收式制冷机中的高压冷剂阀开度u_hgr；输入变量服从以下约束：0≤u_mbf,u_re,u_hgr≤1；选定微型燃气轮机的燃料流量控制阀门u_mbf、再热阀门开度u_re以及双效溴化锂吸收式制冷机中的高压冷剂阀开度u_hgr这三个输入量作为控制系统的控制量。

步骤1-3中，三输入三输出系统的离散线性定常状态空间模型为：

$\begin{array}{c}x(k+1)=Ax\left(k\right)+Bu\left(k\right)\\ y\left(k\right)=Cx\left(k\right)\end{array}$ 公式(1)，

式中，u(k)∈R^m，y(k)∈R^l分别是k时刻的输入列向量、输出列向量，输入向量u＝[u_mbf,u_re,u_hgr]^T，输出向量y＝[y_n,y_Tclo,y_Thwo]^T，x(k)∈Rⁿ是k时刻的状态向量，x(k+1)∈Rⁿ是k+1时刻的状态向量；m、l和n分别是输入向量、输出向量和状态向量的维数，R^m，R^l和Rⁿ分别是m、l和n维实空间，分别表示全体m、l和n维实向量构成的集；A是系统矩阵，表示了系统内部状态变量之间的联系；B是控制矩阵，表示各个输入变量如何控制状态变量；C为输出矩阵即观测矩阵，输出变量反映状态变量。采用的状态空间模型中的直接转移矩阵为零矩阵，该模型的参数A，B和C矩阵均可通过步骤1-3辨识得到，本发明中采用的模型辨识方法是子空间辨识方法。

步骤2-2中，状态空间模型为公式(1)。

步骤2-2中，推导出未来P个时刻的输出向量的表达式为：

$\hat{Y}\left(k\right)=S_{x}x\left(k\right)+S_{u1}u(k-1)+S_u\mathrm{\Δ}U\left(k\right)$ 公式(2)，

式中，是在k时刻对未来P个时刻(k+1时刻到k+P时刻)的输出向量值的预估值；ΔU(k)＝[Δu(k),Δu(k+1),…,Δu(k+M-1)]^T是未来M个时刻(k时刻到k+M-1时刻)的输入增量值，S_x∈R^lP×n,S_u1∈R^lP×m和S_u∈R^lP×mM分别是由状态空间模型推导出的预测矩阵。S_x∈R^lP×n表示当前时刻(k时刻)的状态向量x(k)如何影响未来P个时刻的输出向量预估值；S_u1∈R^lP×m表示上一时刻(k-1时刻)的输入向量u(k-1)如何影响未来P个时刻的输出向量预估值；S_u∈R^lP×mM表示未来M个时刻(k时刻到k+M-1时刻)的输入增量值如何影响未来P个时刻的输出向量预估值。R^lP×n，R^lP×m和R^lP×mM分别是lP×n,lP×m和lP×mM维实向量空间，分别表示全体lP×n,lP×m和lP×mM维实矩阵构成的集。m,l和n分别是输入向量、输出向量和状态向量的维数。P和M分别是预测时域和控制时域。

$S_x={\left[\begin{array}{cccc}CA& {\mathrm{CA}}^2& ...& {\mathrm{CA}}^P\end{array}\right]}^T,$

$S_{u1}={\left[\begin{array}{cccc}CB& CB+CAB& ...& \underset{h=0}{\overset{P-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B\end{array}\right]}^T,$

P和M分别是预测时域和控制时域。表示将h₁＝0,1，…，P-1时的CA^hB相加得到的数值，h₁＝0,1，…，P-1，表示将h₂＝0,1，…，P-M时的CA^hB相加得到的数值，h₂＝0,1，…，P-M。

推导未来P个时刻的输出向量的具体过程如下：

根据公式(1)，可推导出：

k时刻的输出向量估计值为

x(k+1)＝Ax(k)+Bu(k)。

k+1时刻的输出向量估计值和表示为：

$\begin{array}{c}\hat{y}(k+1|k)\\ =Cx(k+1)\\ =C\left(Ax\right(k)+Bu(k\left)\right)\\ =CAx\left(k\right)+CBu\left(k\right)\end{array}.$

k+2时刻的状态向量值为：

x(k+2)＝Ax(k+1)+Bu(k+1)

＝A(A_x(k)+B_u(k))+B_u(k+1)。

＝A²x(k)+ABu(k)+Bu(k+1)

k+2时刻的输出向量估计值为：

$\begin{array}{c}\hat{y}(k+2|k)\\ =Cx(k+2)\\ =C\left(A^{2}x\right(k)+ABu(k)+Bu(k+1\left)\right)\\ ={\mathrm{CA}}^{2}x\left(k\right)+CABu\left(k\right)+CBu(k+1)\end{array}.$

k+3时刻的状态向量值为：

x(k+3)＝Ax(k+2)+Bu(k+2)

＝A(A²x(k)+ABu(k)+Bu(k+1))+Bu(k+2)。

＝A³x(k)+A²Bu(k)+ABu(k+1)+Bu(k+2)

k+2时刻的输出向量估计值为：

$\begin{array}{c}\hat{y}(k+3|k)\\ =Cx(k+3)\\ =C\left(A^{3}x\right(k)+A^{2}Bu(k)+ABu(k+1)+Bu(k+2\left)\right)\\ ={\mathrm{CA}}^{3}x\left(k\right)+{\mathrm{CA}}^{2}Bu\left(k\right)+CABu(k+1)+CBu(k+2)\end{array}.$

依此类推，在k时刻对k+i时刻的输出向量估计值为：

$\begin{array}{c}\hat{y}(k+i|k)\\ =Cx(k+i)\\ =C\[A^{i}x\left(k\right)+\underset{h=0}{\overset{i-1}{\mathrm{\Σ}}}{A}^{i-h-1}\left(Bu\right(k+h\left)\right)\]\\ ={\mathrm{CA}}^{i}x\left(k\right)+\underset{h=0}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{i-h-1}\left(Bu\right(k+h\left)\right)\\ ={\mathrm{CA}}^{i}x\left(k\right)+\underset{h=0}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{i-h-1}\left(Bu\right(k-1)+\underset{j=0}{\overset{h}{\mathrm{\Σ}}}\mathrm{\Δ}u(k+j\left)\right)\\ ={\mathrm{CA}}^{i}x\left(k\right)+\underset{h=0}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{i-h-1}\left(B\right(u\left(k-1\right)+\mathrm{\Δ}u\left(k\right)+\mathrm{\Δ}u\left(k+1\right)+...+\mathrm{\Δ}u\left(k+h\right)\left)\right)\\ ={\mathrm{CA}}^{i}x\left(k\right)+\underset{j=0}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{i-j-1}Bu(k-1)+\underset{j=0}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{i-j-1}B\mathrm{\Δ}u\left(k\right)+\underset{j=0}{\overset{i-2}{\mathrm{\Σ}}}{\mathrm{CA}}^{i-j-2}B\mathrm{\Δ}u(k+1)+...+\underset{j=0}{\overset{0}{\mathrm{\Σ}}}{\mathrm{CA}}^{0}B\mathrm{\Δ}u(k+i-1)\\ ={\mathrm{CA}}^{i}x\left(k\right)+\underset{h=0}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}Bu(k-1)+\underset{h=0}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B\mathrm{\Δ}u\left(k\right)+\underset{h=0}{\overset{i-2}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B\mathrm{\Δ}u(k+1)+...+CB\mathrm{\Δ}u(k+i-1)\\ ={\mathrm{CA}}^{i}x\left(k\right)+\underset{h=0}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}Bu(k-1)+\[\begin{array}{cccc}\underset{h=0}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B& \underset{h=0}{\overset{i-2}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B& ...& CB\end{array}\]\left[\begin{array}{c}\mathrm{\Δ}u\left(k\right)\\ \mathrm{\Δ}u(k+1)\\ .\\ .\\ .\\ \mathrm{\Δ}u(k+i-1)\end{array}\right]\end{array}$

则在k时刻对未来P个时刻(k+1时刻到k+P时刻)的输出向量值的预估值可以表示为：

$\begin{array}{c}\hat{Y}\left(k\right)=\left[\begin{array}{c}\hat{y}(k+1|k)\\ \hat{y}(k+2|k)\\ .\\ .\\ .\\ \hat{y}(k+P|k)\end{array}\right]\\ =\left[\begin{array}{c}CA\\ {\mathrm{CA}}^2\\ .\\ .\\ .\\ {\mathrm{CA}}^P\end{array}\right]x\left(k\right)+\left[\begin{array}{c}CB\\ CB+CAB\\ .\\ .\\ .\\ \underset{h=0}{\overset{P-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B\end{array}\right]u(k-1)+\left[\begin{array}{c}CB\\ \begin{array}{cc}CB+CAB& CB\end{array}\\ .\\ .\\ .\\ \begin{array}{cccc}\underset{h=0}{\overset{P-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B& ...& ...& \underset{h=0}{\overset{P-M}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B\end{array}\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}u\left(k\right)\\ \mathrm{\Δ}u(k+1)\\ .\\ .\\ .\\ \mathrm{\Δ}u(k+M-1)\end{array}\right]\end{array}.$

步骤2-3中，优化性能指标使该系统在未来P个时刻的输出预测值最大程度接近输出向量设定值R_P；R_P＝[R(k+1)^T,R(k+2)^T,…,R(k+P)^T]^T。

在性能指标中加入对ΔU(k)的软约束，性能指标取作：0

$\min J={\left(\hat{Y}\right(k)-R_P)}^T{Q}_y\left(\hat{Y}\right(k)-R_P)+\mathrm{\Δ}U{\left(k\right)}^T{W}_u\mathrm{\Δ}U\left(k\right)$ 公式(3)，

基于阀门开度的约束0≤u_mbf,u_re,u_hgr≤1，ΔU(k)满足以下约束：

$\left[\begin{array}{c}I\\ I\\ .\\ .\\ .\\ I\end{array}\right](u_{min}-u(k-1\left)\right)\≤\mathrm{\ζ}\mathrm{\Δ}U\left(k\right)\≤\left[\begin{array}{c}I\\ I\\ .\\ .\\ .\\ I\end{array}\right](u_{max}-u(k-1\left)\right)$ 公式(4)，

式中，umin＝[0,0,0]^T，u_max＝[1,1,1]^T，

步骤2-4中，确定当前时刻的控制量：ΔU(k)的求解即以ΔU(k)为优化变量，在公式(2)下使公式(3)最小，用二次规划法求解，得到未来M个时刻的输入增量值ΔU(k)，当前时刻的控制量可确定为u(k)＝u(k-1)+Δu(k)。

本发明基于子空间辨识和模型预测控制的微型燃气轮机冷热电联供系统的协调控制方法，能够处理被控量约束，大延迟以及多变量系统等控制方面的问题，能够进一步提高对用户负荷的响应速度、降低动态偏差，从而提高供电、供热及制冷的品质。

本发明具有以下优点：由于在协调预测控制的优化性能指标中包含对控制增量的软约束，各阀门开度变化平缓；协调预测控制算法很好地解决了PID难以处理的输入量约束问题，各控制量均在合理范围内；由于解耦PID控制器在处理系统大惯性和延迟方面存在不足，所以被控量会出现超调，造成汽轮机转速y_n(r/min)、冷媒水温度y_Tclo(℃)以及生活热水温度y_Thwo(℃)的波动，而协调预测控制算法在预测模型中考虑了系统惯性和延迟的影响，可以良好地跟踪用户侧负荷。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1MGT-CCHP系统的原理图。

图2基于协调预测控制(coordinateMPC)的微型燃气轮机的冷热电三联供系统(MGT-CCHP)控制框图。

图3微型燃气轮机冷热电联供机组输入输出结构图。

图4基于预测控制算法和基于解耦PID控制算法的MGT-CCHP控制系统性能比较(输出量)示意图。

图5基于预测控制算法和基于解耦PID控制算法的MGT-CCHP控制系统性能比较(输入量)示意图。

具体实施方式

实施例1

本发明的控制方法应用于一台80kW回热微型燃气轮机以及一台425kW双效溴化锂吸收式制冷机组成的MGT-CCHP系统的仿真模型，如图1所示。此三联供系统可以输出电能，同时还能提供60-80℃的生活热水和7-12℃的冷媒水(用来制冷)。控制目标是实现对用户侧供电、供热以及制冷负荷的良好跟踪

本发明提出一种微型燃气轮机冷热电三联供系统的协调预测控制方法，能够处理被控量约束，大延迟以及多变量系统等控制方面的问题，能够进一步提高对用户负荷的响应速度、降低动态偏差，从而提高供电、供热及制冷的品质。包括如下步骤：

(一)离线步骤：

1)简化被控系统：保持冷媒水流量，生活热水流量以及相应入口温度都保持恒定；保证冷剂泵始终全开；制冷机中溶液循环泵处于自动运行模式，以保证高低压发生器水位都位于安全运行范围；微型燃气轮机切换到带变速负载运行，如此燃料阀可以切换至手动状态。在保证上述前提下，微型燃气轮机冷热电三联供系统方可简化为一个三输入变量三输出变量的系统。三个输出变量是汽轮机转速y_n(r/min)、冷媒水温度y_Tclo(℃)以及生活热水温度y_Thwo(℃)，三个输入变量分别是微型燃气轮机的燃料流量控制阀门u_mbf、再热阀门开度u_re以及双效溴化锂吸收式制冷机中的高压冷剂阀开度u_hgr。输入变量服从以下约束：0≤u_mbf,u_re,u_hgr≤1。选定微型燃气轮机的燃料流量控制阀门u_mbf、再热阀门开度u_re以及双效溴化锂吸收式制冷机中的高压冷剂阀开度u_hgr这三个输入量作为控制系统的控制量。

2)在MGT-CCHP系统处于稳定运行状态时，依次在各输入变量中加入高斯白噪声扰动，采集输入输出变量的数据作为模型辨识的数据。

3)根据扰动实验的数据，利用子空间辨识的方法辨识出该三输入三输出系统的离散状态空间模型：

$\begin{array}{c}x(k+1)=Ax\left(k\right)+Bu\left(k\right)\\ y\left(k\right)=Cx\left(k\right)\end{array}$ (1)

式中，u(k)∈R^m，y(k)∈R^l分别是k时刻的输入向量、输出向量。输入向量u＝[u_mbf,u_re,u_hgr]^T，输出向量y＝[y_n,y_Tclo,y_Thwo]^T。A、B、C是系统矩阵。x(k)∈Rⁿ是k时刻的状态向量，由于子空间辨识方法辨识出的状态量是没有实际的物理意义的，其数值不可以通过测量得到，可以设计状态观测器观测得到。辨识出的状态空间模型是线性的。

(二)在线步骤：

1)在每一时刻，利用卡尔曼滤波器来观测当前时刻(假设是k时刻)的状态向量x(k)的值。

2)在每一时刻，根据当前时刻(假设是k时刻)的状态向量x(k)、k-1时刻的输入向量u(k-1)以及辨识出的状态空间模型(1)，推导出未来P个时刻的输出向量的表达式：

$\hat{Y}\left(k\right)=S_{x}x\left(k\right)+S_{u1}u(k-1)+S_u\mathrm{\Δ}U\left(k\right)---\left(2\right)$

式中， $\hat{Y}\left(k\right)={\[\hat{y}(k+1|k),\hat{y}(k+2|k),...,\hat{y}(k+P|k)\]}^T$ 是未来P个时刻的输出向量值；ΔU(k)＝[Δu(k),Δu(k+1),…,Δu(k+M-1)]^T是未来M个时刻的输入增量值。P和M分别是预测时域和控制时域。S_x∈R^lP×n,S_u1∈R^lP×m和S_u∈R^lP×mM分别是由状态空间模型推导出的预测矩阵。

S_x＝[CACA²…CA^P]^T,

$S_{u1}={\left[\begin{array}{cccc}CB& CB+CAB& ...& \underset{h=0}{\overset{P-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B\end{array}\right]}^T,$

表示将h₁＝0,1，…，P-1时的CA^hB相加得到的数值，h₁＝0,1，…，P-1，表示将h₂＝0,1，…，P-M时的CA^hB相加得到的数值，h₂＝0,1，…，P-M。

3)优化性能指标：控制目标是实现对用户侧供电、供热以及制冷负荷的良好跟踪。因此，优化性能指标应当使该系统在未来P个时刻的输出预测值尽可能地接近给定的期望值R_P。R_P＝[R(k+1)^T,R(k+2)^T,…,R(k+P)^T]^T，是未来P个时刻的输出向量设定值。考虑到不希望控制增量变化剧烈，在性能指标中加入对ΔU(k)的软约束。因此性能指标取作：

$\min J={\left(\hat{Y}\right(k)-R_P)}^T{Q}_y\left(\hat{Y}\right(k)-R_P)+\mathrm{\Δ}U{\left(k\right)}^T{W}_u\mathrm{\Δ}U\left(k\right)---\left(3\right)$

考虑到阀门开度的约束0≤u_mbf,u_re,u_hgr≤1，ΔU(k)应满足以下约束：

$\left[\begin{array}{c}I\\ I\\ .\\ .\\ .\\ I\end{array}\right](u_{min}-u(k-1\left)\right)\≤\mathrm{\ζ}\mathrm{\Δ}U\left(k\right)\≤\left[\begin{array}{c}I\\ I\\ .\\ .\\ .\\ I\end{array}\right](u_{max}-u(k-1\left)\right)---\left(4\right)$

式中，u_min＝[0,0,0]^T，u_max＝[1,1,1]^T。

4)确定当前时刻的控制量：预测控制是一种以优化确定控制策略的算法。ΔU(k)的求解即以ΔU(k)为优化变量，在动态模型(2)下使性能指标(3)最小的优化问题。根据分析，该优化问题是一个典型的带约束的二次规划问题，用二次规划法求解，得到未来M个时刻的输入增量值ΔU(k)。当前时刻的控制量可确定为u(k)＝u(k-1)+Δu(k)。

实施例2

下面结合附图和具体实施例，进一步阐述本发明。针对一台80kW回热微型燃气轮机以及一台425kW双效溴化锂吸收式制冷机组成的MGT-CCHP系统的仿真模型，本发明设计了一个三输入三输出的协调预测控制系统，使微型燃气轮机侧的燃料流量控制阀门u_mbf和双效溴化锂吸收式制冷机侧的再热阀门开度u_re、高压冷剂阀开度u_hgr能够协调工作，同时有效地处理制冷机侧的大惯性，使三个输出变量(汽轮机转速y_n(r/min)、冷媒水温度y_Tclo(℃)以及生活热水温度y_Thwo(℃))能够良好地追踪用户侧的负荷指令，提高调节品质。

这个协调预测控制方法的控制框图如图2：

图2中的MGT-CCHP模块代表MGT-CCHP系统，此处是一台80kW回热微型燃气轮机以及一台425kW双效溴化锂吸收式制冷机组成的MGT-CCHP系统的机理模型。CoordinatedMPC是指本发明采用的协调预测控制方法，具体算法详见下文的具体实施步骤。

本发明中的协调预测控制方法选取MGT-CCHP微型燃气轮机侧的燃料流量控制阀门u_mbf和双效溴化锂吸收式制冷机侧的再热阀门开度u_re、高压冷剂阀开度u_hgr作为控制变量同时协调工作，来调整三个输出变量：汽轮机转速y_n(r/min)、冷媒水温度y_Tclo(℃)以及生活热水温度y_Thwo(℃)的数值，使得输出的电、热、冷负荷能够满足用户侧的需求汽轮机转速设定值r_n(r/min)、冷媒水温度设定值r_Tclo(℃)以及生活热水温度设定值r_Thwo(℃)。

在每个控制周期内，测量三个输出变量的实际值y_n(r/min)、y_Tclo(℃)和y_Thwo(℃)以及用户侧的设定值r_n(r/min)、r_Tclo(℃)和r_Thwo(℃)，送至协调预测控制器模块(CoordinatedMPC模块)。经协调预测控制模块处理计算得到当前时刻的控制量(燃料流量控制阀门u_mbf和双效溴化锂吸收式制冷机侧的再热阀门开度u_re、高压冷剂阀开度u_hgr)。将计算得到的控制量值应用到对MGT-CCHP模块，经过一个周期运行，得到下一控制周期的实际输出值，再检验实际输出量是否已经满足用户侧的负荷指令，如不满足，重复之前的步骤。

这个协调预测控制方法的具体实施步骤分为离线步骤和在线实时步骤，如下：

(一)离线步骤

1)简化被控系统：保持冷媒水流量，生活热水流量以及相应入口温度都保持恒定；保证冷剂泵始终全开；制冷机中溶液循环泵处于自动运行模式，以保证高低压发生器水位都位于安全运行范围；微型燃气轮机切换到带变速负载运行，如此燃料阀可以切换至手动状态。在保证上述前提下，微型燃气轮机冷热电三联供系统方可简化为一个三输入变量三输出变量的系统，如图3。

2)离线辨识预测模型。

①在MGT-CCHP系统处于稳定运行状态时，依次在各输入变量中加入高斯白噪声扰动，采集输入输出变量的数据作为模型辨识的数据，采样时间为2s，扰动信号的变化时间为100s，整个扰动过程持续30000s。得到输入向量和输出向量的时间序列X和Y，均包含15000组数据，为了方便辨识，将各数据标称化。

②采用子空间辨识的方法对系统进行辨识。利用子空间辨识的方法辨识出该三输入三输出系统的离散线性状态空间模型：

x(k+1)＝Ax(k)+Bu(k)

y(k)＝Cx(k)

式中，u(k)∈R^m，y(k)∈R^l分别是k时刻的输入向量、输出向量。输入向量u＝[u_mbf,u_re,u_hgr]^T，输出向量y＝[y_n,y_Tclo,y_Thwo]^T。A、B、C是系统矩阵。x(k)∈Rⁿ是k时刻的状态向量，由于子空间辨识方法辨识出的状态量是没有实际的物理意义的，其数值不可以通过测量得到，可以设计状态观测器观测得到。辨识出的状态空间模模型验证。比较辨识出的预测模型和仿真模型的阶跃响应进行比较，检验预测模型的动态和稳态特性。

在本案例，辨识出的状态空间模型的阶数n为10。由于是三输入三输出系统，l＝m＝3。因为阶数比较高，所以系统矩阵A、B、C在此不提供具体数值。

(二)在线步骤

1)在每一时刻，利用卡尔曼滤波器来观测当前时刻(假设是k时刻)的状态向量x(k)的值。利用卡尔曼滤波器观测状态值的具体步骤如下：

离线部分：

确定过程白噪声和测量白噪声的协方差矩阵Q_w和R_v。

初始时刻的状态估计值和误差协方差矩阵P(0|0)。

本案例中，令P(0|0)＝O_n×n。令Q_w＝I_n×n，R_v＝I_l×l

在线实时部分：

①获取上一时刻(k-1时刻)状态估计值输入向量u(k-1)、误差协方差矩阵P(k-1|k-1)和当前时刻的输出y(k)；

②根据状态空间方程，先初步估计当前时刻(k时刻)的状态估计值和误差协方差矩阵P(k|k-1)。

$\hat{x}\left(k\right|k-1)=A\hat{x}(k-1|k-1)+Bu(k-1),$

P(k|k-1)＝AP(k-1|k-1)A^T+Q_w。

③计算当前时刻的实际输出y(k)和估计值的偏差e(k)：

$e\left(k\right)=y\left(k\right)-C\hat{x}\left(k\right|k-1).$

④重新估计当前时刻(k时刻)的状态，在步骤②中的状态估计值的基础上加上修正项。

$\hat{x}\left(k\right|k)=\hat{x}\left(k\right|k-1)+Ke\left(k\right)$

式中，K＝P(k|k-1)C^T(CP(k|k-1)C^T+R_v)^-1。

⑤重新估计当前时刻(k时刻)的误差协方差矩阵，在步骤②中的误差协方差矩阵P(k|k-1)的基础上加上修正项。

P(k|k)＝(I-KC)P(k|k-1)

2)在每一时刻，根据当前时刻(假设是k时刻)的状态向量x(k)(取作)、k-1时刻的输入向量u(k-1)以及辨识出的状态空间模型，推导出未来P个时刻的输出向量的表达式：

$\hat{Y}\left(k\right)=S_{x}x\left(k\right)+S_{u1}u(k-1)+S_u\mathrm{\Δ}U\left(k\right)$

式中， $\hat{Y}\left(k\right)={\[\hat{y}(k+1|k),\hat{y}(k+2|k),...,\hat{y}(k+P|k)\]}^T$ 是未来P个时刻的输出向量值；

ΔU(k)＝[Δu(k),Δu(k+1),…,Δu(k+M-1)]^T是未来M个时刻的输入增量值。P和M分别是预测时域和控制时域。S_x∈R^lP×n,S_u1∈R^lP×m和S_u∈R^lP×mM分别是由状态空间模型推导出的预测矩阵。

S_x＝[CACA²…CA^P]^T,

$S_{u1}={\left[\begin{array}{cccc}CB& CB+CAB& ...& \underset{h=0}{\overset{P-1}{\mathrm{\Σ}}}{\mathrm{CA}}^{h}B\end{array}\right]}^T,$

选取P＝100,M＝20，由于预测矩阵的阶数较大，此处不提供各状态空间矩阵的具体数值。表示将h₁＝0,1，…，P-1时的CA^hB相加得到的数值，h₁＝0,1，…，P-1，表示将h₂＝0,1，…，P-M时的CA^hB相加得到的数值，h₂＝0,1，…，P-M。

3)在每一时刻，获取用户侧的负荷指令R_P＝[R(k+1)^T,R(k+2)^T,…,R(k+P)^T]^T，列出优化性能指标： $\min J={\left(\hat{Y}\right(k)-R_P)}^T{Q}_y\left(\hat{Y}\right(k)-R_P)+\mathrm{\Δ}U{\left(k\right)}^T{W}_u\mathrm{\Δ}U\left(k\right)$

其中输出权矩阵和输入增量权矩阵分别取作 $Q_y=I_P\⊗diag(\[1,1,1\]),W_u=I_M\⊗diag(\[0.5,$ $0.5,0.5\]).$

4)在每一时刻，确定当前时刻的控制量。预测控制是一种以优化确定控制策略的算法。ΔU(k)的求解即以ΔU(k)为优化变量，在动态模型(2)下使性能指标(3)最小的优化问题。根据分析，该优化问题是一个典型的带约束的二次规划问题，用二次规划法求解，得到未来M个时刻的输入增量值ΔU(k)。当前时刻的控制量可确定为u(k)＝u(k-1)+Δu(k)。本发明中，利用Matlab中的二次规划工具箱对步骤5中的优化性能指标进行求解。

本实施例为了比较本发明中的协调预测控制方法(CoordinatedMPC)和解耦PID控制方法的控制效果，做了三组阶跃响应实验。仿真过程共计2500s。在0s时，MGT-CCHP系统稳态运行在额定工况，额定工况时的稳态值如下所示：

输入量额定工况稳态值：y＝[y_n,y_Tclo,y_Thwo]^T＝[68000r/min,7℃,80℃]^T。

输出额定工况稳态值：u＝[u_mbf,u_re,u_hgr]^T＝[0.7465,0.3,0.2]^T。

在100s、900s和1700s，分别令系统的负荷设定值r_n(r/min)、r_Tclo(℃)和r_Thwo(℃)产生阶跃(+/-10％的幅度)，每次都在系统稳定后才进行下一次阶跃实验。比较协调预测和传统解耦PID控制的MGT-CCHP系统的控制效果，被控量变化见图4，控制量变化见图5。解耦PID控制方法的控制结构和控制器参数均参考文献张俊礼.微型燃气轮机冷热电联供系统集成优化与动态运行规律研究[D].东南大学,2015.。

由图4、图5可以看出：

当转速设定值在100s阶跃-10％时，协调预测控制回路的转速经历小幅震荡(±300r/min)50s后达到稳定，冷媒水温度和生活热水温度在小幅波动(±0.1℃)达到稳定。解耦PID控制回路的转速经历衰减震荡(±2500r/min)，80s达到稳定，，冷媒水温度和生活热水温度在小幅波动(±0.15℃)达到稳定。在转速设定值阶跃下降的瞬间，燃料量迅速减少；同时回热阀迅速关小，减缓燃料量增加对下游制冷机的影响。

当冷媒水温度设定值在900s阶跃10％时，协调预测控制回路的冷媒水温度平稳上升，268s后达到稳定值，转速值在下降500r/min后回升达至稳定值，生活热水温度保持不变。解耦PID控制回路的冷媒水温度上升后回落，超调量达到0.08℃后，共计600s达到稳定，转速值在上升1300r/min后回升达至稳定值，生活热水温度保持不变。有两种手段可以使冷媒水温度上升，一是关小高压冷剂阀(导致生活热水温度上升)，二是减少制冷机进口的烟气流量(导致生活热水温度下降)。为了使得在冷媒水温度上升的同时保持生活热水温度不变，两种控制器均在关小高压冷剂阀的同时开大回热阀和减少燃料流量以减少制冷机进口的烟气流量，也使得转速值在小幅震荡后恢复稳定。由于制冷机的热惯性较大，而且PID控制器对于大惯性系统难以同时兼顾稳定性和快速性，所以PID控制器的效果不佳。

当生活热水温度设定值在1700s阶跃-10％时，协调预测控制回路的生活热水温度平稳下降，150s后达到稳定值，转速值先下降340r/min再回升至稳定状态达到稳定值，冷媒水温度先下降0.3℃后回升到稳定。解耦PID控制回路的转速先下降后回升，超调量达到0.7℃后，共计200s达到稳定，转速值经历大幅震荡(-1200r/min至7000r/min)至稳定值，冷媒水温度经历小幅震荡(-0.3℃至0.13℃)后稳定。有两种手段可以使生活热水温度下降，一是开大高压冷剂阀(导致冷媒水温度下降)，二是减少制冷机进口的烟气流量(导致冷媒水温度上升)。两种控制器同时采用两种手段使生活热水温度降温，同时保持冷媒水温度保持不变。由于制冷机的热惯性比较大，生活热水温度对各控制量变化的响应很慢，PID在处理大惯性系统方面有不足，而协调预测控制器在每个控制周期都对未来P个时刻的输出值进行预估，考虑了系统的惯性，所以控制效果优于解耦PID控制器，过程中没有震荡。

与解耦PID控制器相比，本发明中的协调预测控制算法具有以下优点：由于在协调预测控制的优化性能指标中包含对控制增量的软约束，各阀门开度变化平缓；协调预测控制算法很好地解决了PID难以处理的输入量约束问题，各控制量均在合理范围内；由于解耦PID控制器在处理系统大惯性和延迟方面存在不足，所以被控量会出现超调，造成汽轮机转速y_n(r/min)、冷媒水温度y_Tclo(℃)以及生活热水温度y_Thwo(℃)的波动，而协调预测控制算法在预测模型中考虑了系统惯性和延迟的影响，可以良好地跟踪用户侧负荷。

本发明提供了一种微型燃气轮机冷热电三联供系统的协调控制方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (7)

1.一种微型燃气轮机冷热电三联供系统的协调控制方法，其特征在于，包括离线阶段和在线阶段两个部分；

离线阶段包括：

步骤1-1，简化被控系统；

步骤1-2，在MGT-CCHP系统处于稳定运行状态时，依次在各输入变量中加入高斯白噪声扰动，采集输入变量和输出变量的数据作为模型辨识的数据；

步骤1-3，根据扰动实验的数据，利用子空间辨识的方法辨识出该三输入三输出系统的离散状态空间模型；

在线阶段包括：

步骤2-1，利用卡尔曼滤波器来观测当前时刻k的状态向量x(k)的值；

步骤2-2，据当前时刻k的状态向量x(k)、k-1时刻的输入向量u(k-1)以及辨识出的状态空间模型，推导出未来P个时刻的输出向量的表达式；

步骤2-3，优化性能指标；

步骤2-4，利用二次规划的方法，求解优化性能指标最小的优化问题，得到当前时刻的控制增量。

2.根据权利要求1所述的方法，其特征在于，步骤1-1中，简化被控系统包括：

保持冷媒水流量、活热水流量以及相应入口温度三者恒定；保证冷剂泵始终全开；制冷机中溶液循环泵处于自动运行模式；微型燃气轮机切换到带变速负载运行；

将微型燃气轮机冷热电三联供系统方简化为一个三输入变量三输出变量的系统，三个输出变量是汽轮机转速y_n、冷媒水温度y_Tclo以及生活热水温度y_Thwo，三个输入变量分别是微型燃气轮机的燃料流量控制阀门u_mbf、再热阀门开度u_re以及双效溴化锂吸收式制冷机中的高压冷剂阀开度u_hgr；输入变量服从以下约束：0≤u_mbf,u_re,u_hgr≤1；选定微型燃气轮机的燃料流量控制阀门u_mbf、再热阀门开度u_re以及双效溴化锂吸收式制冷机中的高压冷剂阀开度u_hgr这三个输入量作为控制系统的控制量。

3.根据权利要求1所述的方法，其特征在于，步骤1-3中，三输入三输出系统的离散线性定常状态空间模型为：

$\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bu}\left(k\right)\\ y\left(k\right)=\mathrm{Cx}\left(k\right)\end{array}$ 公式(1)，

式中，u(k)∈R^m，y(k)∈R^l分别是k时刻的输入列向量、输出列向量，输入向量u＝[u_mbf,u_re,u_hgr]^T，输出向量y＝[y_n,y_Tclo,y_Thwo]^T，x(k)∈Rⁿ是k时刻的状态向量，x(k+1)∈Rⁿ是k+1时刻的状态向量；m、l和n分别是输入向量、输出向量和状态向量的维数，R^m，R^l和Rⁿ分别是m、l和n维实空间，分别表示全体m、l和n维实向量构成的集；A是系统矩阵，表示了系统内部状态变量之间的联系；B是控制矩阵，表示各个输入变量如何控制状态变量；C为输出矩阵即观测矩阵，输出变量反映状态变量。

4.根据权利要求3所述的方法，其特征在于，步骤2-2中，状态空间模型为公式(1)。

5.根据权利要求3所述的方法，其特征在于，步骤2-2中，推导出未来P个时刻的输出向量的表达式为：

$\hat{Y}\left(k\right)=S_{x}x\left(k\right)+S_{u1}u(k-1)+S_u\mathrm{\Δ}U\left(k\right)$ 公式(2)，

且，式中， $\hat{Y}\left(k\right)={\[\hat{y}(k+1|k),\hat{y}(k+2|k),...,\hat{y}(k+P|k)\]}^T;$

ΔU(k)＝[Δu(k),Δu(k+1),…,Δu(k+M-1)]^T，

S_x＝[CACA²…CA^P]^T,

$S_{u1}={\left[\begin{array}{cccc}CB& CB+CAB& ...& \underset{h=0}{\overset{P-1}{\Σ}}{\mathrm{CA}}^{h}B\end{array}\right]}^T,$

P和M分别是预测时域和控制时域，ΔU(k)是未来M个时刻的输入增量值，S_x∈R^lP×n,S_u1∈R^lP×m和S_u∈R^lP×mM分别是由状态空间模型推导出的预测矩阵，S_x∈R^lP×n表示当前时刻即k时刻的状态向量x(k)影响未来P个时刻的输出向量预估值，S_u1∈R^lP×m表示上一时刻即k-1时刻的输入向量u(k-1)影响未来P个时刻的输出向量预估值；S_u∈R^lP×mM表示未来M个时刻即k时刻到k+M-1时刻的输入增量值影响未来P个时刻的输出向量预估值；R^lP×n，R^lP×m和R^lP×mM分别是lP×n,lP×m和lP×mM维实向量空间，分别表示全体lP×n,lP×m和lP×mM维实矩阵构成的集；m、l和n分别是输入向量、输出向量和状态向量的维数，P和M分别是预测时域和控制时域，表示将h₁＝0,1，…，P-1时的CA^hB相加得到的数值，h₁＝0,1，…，P-1，表示将h₂＝0,1，…，P-M时的CA^hB相加得到的数值，h₂＝0,1，…，P-M。

6.根据权利要求3所述的方法，其特征在于，步骤2-3中，优化性能指标使该系统在未来P个时刻的输出预测值最大程度接近输出向量设定值R_P；

R_P＝[R(k+1)^T,R(k+2)^T,…,R(k+P)^T]^T，

在性能指标中加入对ΔU(k)的软约束，性能指标按照如下公式：

$\min J={\left(\hat{Y}\right(k)-R_P)}^T{Q}_y\left(\hat{Y}\right(k)-R_P)+\mathrm{\Δ}U{\left(k\right)}^T{W}_u\mathrm{\Δ}U\left(k\right)$ 公式(3)，

基于阀门开度的约束0≤u_mbf,u_re,u_hgr≤1，ΔU(k)满足以下约束：

$\left[\begin{array}{c}I\\ I\\ \·\\ \·\\ \·\\ I\end{array}\right](u_{min}-u(k-1\left)\right)\≤\mathrm{\ζ}\mathrm{\Δ}U\left(k\right)\≤\left[\begin{array}{c}I\\ I\\ \·\\ \·\\ \·\\ I\end{array}\right](u_{max}-u(k-1\left)\right)$ 公式(4)，

式中，u_min＝[0,0,0]^T，u_max＝[1,1,1]^T。

7.根据权利要求6所述的方法，其特征在于，步骤2-4中，确定当前时刻的控制量：ΔU(k)的求解即以ΔU(k)为优化变量，在公式(2)下使公式(3)最小，用二次规划法求解，得到未来M个时刻的输入增量值ΔU(k)，当前时刻的控制量确定为：u(k)＝u(k-1)+Δu(k)。