CN104898426A - 基于梯度下降法和广义预测控制的室温回路控制方法 - Google Patents

Abstract

本发明公开了一种基于梯度下降法和广义预测控制的室温回路控制方法，包括如下步骤：（1）、确定预设室温的设定值；（2）、采用梯度下降法在线辨识出变风量空调系统的室温回路的模型参数；（3）、在变风量空调系统中加入广义预测控制器，对模型参数进行在线估计，实现自适应控制。本发明该控制方式采用了梯度下降法能够在线辨识出空调系统室温回路的模型参数，并采用广义预测控制方式进行控制，易于在线估计参数，实现自适应控制。

Description

基于梯度下降法和广义预测控制的室温回路控制方法

技术领域

本发明涉及一种基于梯度下降法和广义预测控制的室温回路控制方法。

背景技术

随着建筑能耗的加剧，空调系统的节能已经越来越受到人们的关注，20世纪60年代变风量空调系统诞生于美国，相比于定风量空调系具有节能、舒适等优点，因此获得了广泛的推广。与此同时，人们也意识到空调系统的能耗与其控制系统的性能息息相关，但由于变风量空调系统具有非线性、大时滞等特性在控制方面有较大难度，因此开发适合于变风量空调系统的先进控制算法已成为当前较为热门的研究方向。

广义预测控制算法于1987年被Clarke等人提出，能够克服广义最小方差、极点配置等各种自适应控制方法的不足，具有较强的鲁棒性，并能够适用于具有纯时延、时变参数的控制过程，因此很快被应用于空调系统的控制领域。

在广义预测控制算法中提到要针对预测模型进行控制，但必须要首先知道预测模型中的A(z^-1)和B(z^-1)才可以进行运算，尤其是当模型受到外界扰动具有时变特性时，这就需要对系统模型进行辨识，先在线估计出，然后用参数估计值代替真实值进行控制律的推导计算。所以，模型反馈校正的优劣直接影响广义预测自适应控制的性能，需要选择适当的在线辨识算法，来及时修正模型误差，以提高预测精度，保证控制效果最佳。

广义预测控制一般采用递推最小二乘法在线估计系统参数，使用最小二乘法的过程中不可避免的遇到矩阵求逆的情况，但矩阵求逆有可能会存在数值不稳定的情况(比如对希尔伯特矩阵求逆就几乎是不可能的)，因而使用最小二乘法具有一定而局限性。

发明内容

本发明提出了一种基于梯度下降法和广义预测控制的室温回路控制方法。

技术方案如下：

一种基于梯度下降法和广义预测控制的室温回路控制方法，包括如下步骤：

(1)、确定预设室温的设定值；

(2)、采用梯度下降法在线辨识出变风量空调系统的室温回路的模型参数；

(3)、在变风量空调系统中加入广义预测控制器，对模型参数进行在线估计，实现自适应控制。

上述步骤(3)中广义预测控制器的控制方法包括：

(21)、预测模型；

(22)、滚动优化；

(23)、反馈校正。

上述步骤(21)中预测模型采用受控自回归积分滑动平均模型，模型如公式(1)所示；

A(z^-1)y(k)＝B(z^-1)u(k-d)+C(z^-1)ξ(k)/Δ    (1)

式中A(z^-1)、B(z^-1)、C(z^-1)分别是n、m和n阶的z^-1的多项式；

$A\left(z^{-1}\right)=1+a_1{z}^{-1}+...+a_{n_a}{z}^{-n_a}$

$B\left(z^{-1}\right)=b_0+b_1{z}^{-1}+...+a_{n_b}{z}^{-n_b}$

$C\left(z^{-1}\right)=1+c_1{z}^{-1}+...+c_{n_c}{z}^{-n_c}$

其中Δ＝1-z^-1；y(k)、u(k)、ξ(k)分别表示输出、输入和均值为零的白噪声序列，d为系统的时滞，C(z^-1)＝1。

上述步骤(22)中滚动优化的步骤如下：

步骤401：设定目标函数，采用如下目标函数，

$J=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{[y(k+j)-w(k+j)]}^2+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{\λ}\left(j\right){\left[\mathrm{\Δu}(k+j-1)\right]}^2---\left(2\right)$

式中n为最大预测长度；m称为控制长度(m≤n)，y(k+j)为(k+j)时的实际室温，w(k+j)为k+j时刻的室温设定值，Δu(k+j-1)为k+j-1时刻的模型输入与前一时刻的模型输入的差值，即Δu(k+j-1)＝u(k+j-1)-u(k+j；λ(j)是大于零的控制加权系数，取λ(j)为常数λ；

采用柔化控制，跟踪参考轨迹线，参考轨迹线如公式3所示；

w(k+j)＝a^jy(k)+(1-a^j)y_r     (3)

式中y_r、y(k)和w(k+j)分别为设定值、输出值和参考轨迹线；a为柔化系数，0＜a＜1；

步骤402：输出预测，

根据广义预测理论，引入丢番图方程：

1＝E_j(z^-1)A(z^-1)Δ+z^-JF_j(z^-1)     (4)

E_j(z^-1)B(z^-1)＝G_j(z^-1)+z^-jH_j(z^-1)      (5)

式中E_j(z^-1)＝e₀+e₁z^-1+…+e_j-1z^-j+1

$F_j\left(z^{-1}\right)={f_0}^j+{f_1}^j{z}^{-1}+...+{f_{n_a}}^j{z}^{{-n}_a}$

G_j(z^-1)＝g₀+g₁z^-1+…+g_j-1z^-j+1

$H_j\left(z^{-1}\right)={h_0}^j+{h_1}^j{z}^{-1}+...+{f_{n_b-1}}^j{z}^{{-n}_b-1}$

将公式(1)等号两侧同时乘以E_j(z^-1)Δ，并将公式(4)、(5)带入可得：

y(k+j)＝G_jΔu(k+j-1)+F_jy(k)+H_jΔu(k-1)+E_jξ(k+j)   (6)

公式(6)中的最后一项表示外界的白噪声序列，在预测未来时刻的输出时不考虑，公式(6)即可表示为：

$\hat{y}(k+j)=G_j\mathrm{\Δu}(k+j-1)+F_{j}y\left(k\right)+H_j\mathrm{\Δu}(k-1)---\left(7\right)$

将公式(6)写成向量形式，即为：

Y＝GΔU+Fy(k)+HΔu(k-1)+E     (8)

表示对未来j时刻的预测，将其分为已知量和未知量两部分，用f(k)表示已知量，写成矩阵形式即为：

f＝Fy(k)+HΔu(k-1)       (9)

由公式(8)得知：

$\hat{Y}=\mathrm{G\ΔU}+f---\left(10\right)$

式中 $\hat{Y}{[\hat{y}(k+1),\hat{y}(k+2),...,\hat{y}(k+n)]}^T;$

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

f＝[f(k+1),f(k+2),…,f(k+n)]^T；

F＝[F₁,F₂,…,F_n]^T；

H＝[H₁,H₂,…,H_n]^T；

E＝[E₁w(k+1),E₂w(k+2),…,E_nw(k+n)]；

步骤403：求取最优控制率

若令W＝[w(k+1),w(k+2),…,w(k+n)]^T

则公式(2)表示成：

J＝(Y-W)^T(Y-W)+λΔU^TΔU   (11)

用Y的最优预测值来代替Y，将公式(10)带入公式(11)，并令 $\frac{\∂J}{\∂\mathrm{\ΔU}}=0,$

得到：

ΔU＝(G^TG-λI)^-1G^T(W-f)   (12)

每次将第一个分量加入系统，即

u(k)＝u(k-1)+g^T(W-f)   (13)

其中g^T表示(G^TG-λI)^-1G^T中的第一行向量。

上述步骤(2)中梯度下降法的步骤如下：

选取过程模型如下：

$G_m=\frac{b_m{z}^{-1}}{1-a_m{z}^{-1}}{z}^{-d}---\left(14\right)$

室温回路的模型参数中需要进行辨识的参数θ＝[a_m,b_m,d]^T，即广义预测控制中的A(z^-1)＝1-a_mz^-1，B(z^-1)＝b_mz^-1，和时间延迟d；

过程模型的输出与广义预测控制的误差为：

e(k)＝y(k)-y_m(k)＝y(k)-[a_my(k-1)-z^-db_mu(k-1)]   (15)

定义目标函数为：

$J\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\frac{1}{2}{e}^2\left(i\right)---\left(16\right)$

验证J(θ)的扶梯度方向不断修正θ(k)，直至J(θ)获得最小值，梯度下降法可以成：

$\mathrm{\θ}(k+1)=\mathrm{\θ}\left(k\right)-R\left(k\right)\·\mathrm{grad}\left[J\left(\mathrm{\θ}\right)\right]{|}_{{\mathrm{\θ}}_m\left(k\right)}---\left(17\right)$

式中R(k)为加权矩阵，表示J(θ)关于θ_m(k)的梯度，根据目标函数可知：

$\mathrm{grad}\left[J\left(\mathrm{\θ}\right)\right]{|}_{{\mathrm{\θ}}_m\left(k\right)}={[\frac{\∂J}{\∂a_m},\frac{\∂J}{\∂b_m},\frac{\∂J}{\∂d}]}^T\·e\left(k\right);$

根据公式(15)和公式(16)可知：

$\frac{\∂J}{\∂a_m}=-y(k-1)---\left(18\right)$

$\frac{\∂J}{\∂b_m}=-z^{-d}\·u(k-1)---\left(19\right)$

根据Z变换的实域位移定理，推导得

$\frac{\∂J}{\∂b_m}=-u(k-d_m\left(k\right)-1)$

由拉普拉斯变换与Z变换的关系z＝e^ST可知：

$\frac{\∂J}{\∂d}=z^{-d}\·\ln z\·b_m\·u(k-1)=z^{-d}\·s\·T\·b_m\·u(k-1)$

在上式中s为拉普拉斯变换算子，T为采样周期，采用欧拉线性变换式对上式进行化解，推导获得

$\frac{\∂J}{\∂d}==[u(k-d-1)-u(k-d-2)]\·b_m---\left(20\right)$

则

$\underset{\mathrm{\θ}}{\mathrm{grad}}\left[J\left(\mathrm{\θ}\right)\right]{|}_{{\mathrm{\θ}}_m\left(k\right)}=\left[\begin{array}{c}-y(k-1)\\ -u(k-d-1)\\ (u(k-d-1)-u(k-d-2))\·b_m\end{array}\right]\·e\left(k\right)$

根据公式(17)获得在最速下降法情况下空调系统参考对象模型特性参数向量θ_m(k)的递推公式，此时

${\mathrm{\θ}}_m(k+1)={\mathrm{\θ}}_m\left(k\right)-R\left(k\right)\·\left[\begin{array}{c}-y(k-1)\\ -u(k-d-1)\\ (u(k-d-1)-u(k-d-2))\·b_m\end{array}\right]\·e\left(k\right)---\left(21\right)$

令

$h\left(k\right)=\left[\begin{array}{c}-y(k-1)\\ -u(k-d-1)\\ (u(k-d-1)-u(k-d-2))\·b_m\end{array}\right]$

h(k)是整个辨识过程的输入数据向量，加权矩阵R(k)的作用是用来控制各输入数据分量对参数估计值的影响程度，最速下降法的收敛性能直接取决于该加权阵的选择；

选用加权矩阵：

$R^{*}\left(k\right)=\frac{1}{\mathrm{\Σ}{\mathrm{\Λ}}_i\left(k\right)h_i^2\left(k\right)}\mathrm{diag}[{\mathrm{\Λ}}_1\left(k\right),{\mathrm{\Λ}}_2\left(k\right),...,{\mathrm{\Λ}}_N\left(k\right)]---\left(22\right)$

式中的Λ₁(k),Λ₂(k),...,Λ_N(k)表示在加权矩阵中对角线上的变量，在公式(22)中N＝3；

辨识出的时延参数d用整数表示，对d进行取整操作，即：

$d(k+1)=\left\{\begin{array}{c}d\left(k\right)+1(d(k+1)=d\left(k\right)\&GreaterEqual;0.5\\ d\left(k\right)|d(k+1)-d\left(k\right)|<0.5\\ d\left(k\right)-1d(k+1)-d\left(k\right)\&le;-0.5\end{array}\right.$

最终辨识出的1-a_mz^-1，b_mz^-1，d即为广义预测控制器算法中预测模型的A(z^-1)，B(z^-1)和d。

本发明所达到的有益效果：

本发明广义预测控制一般采用递推最小二乘法在线估计系统参数，而本次发明使用了梯度下降法，因为使用最小二乘法的过程中不可避免的遇到矩阵求逆的情况，但矩阵求逆有可能会存在数值不稳定的情况(比如对希尔伯特矩阵求逆就几乎是不可能的)，因而使用最小二乘法具有一定而局限性。而相比之下，梯度下降法虽然有一些弊端，迭代的次数可能也比较高，但是相对来说计算量并不是特别大。而且在处理较大数据量的时候，梯度下降法能够提高预测精度和系统稳定性。

附图说明

图1是本发明的方法示意图；

图2是本发明的流程图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

如图1所示，一种基于梯度下降法和广义预测控制的室温回路控制方法，包括如下步骤：

(1)、确定预设室温的设定值；

(2)、采用梯度下降法在线辨识出变风量空调系统的室温回路的模型参数；

(3)、在变风量空调系统中加入广义预测控制器，对模型参数进行在线估计，实现自适应控制。

上述步骤(3)中广义预测控制器的控制方法包括：

(21)、预测模型；

在广义预测控制理论中，需要有一个描述系统动态行为的基础模型，成为预测模型。他应具有预测功能，即能够根据系统的历史数据和未来的输入，预测系统未来的输出值。广义预测控制一般使用CARIMA模型作为预测模型，即“受控自回归积分滑动平均模型”，模型如公式(1)所示；

A(z^-1)y(k)＝B(z^-1)u(k-d)+C(z^-1)ξ(k)/Δ   (1)

式中A(z^-1)、B(z^-1)、C(z^-1)分别是n、m和n阶的z^-1的多项式；

$A\left(z^{-1}\right)=1+a_1{z}^{-1}+...+a_{n_a}{z}^{-n_a}$

$B\left(z^{-1}\right)=b_0+b_1{z}^{-1}+...+a_{n_b}{z}^{-n_b}$

$C\left(z^{-1}\right)=1+c_1{z}^{-1}+...+c_{n_c}{z}^{-n_c}$

其中Δ＝1-z^-1；y(k)、u(k)、ξ(k)分别表示输出、输入和均值为零的白噪声序列，d为系统的时滞，C(z^-1)＝1。

(22)、滚动优化

步骤401：设定目标函数，

首先需要设定一个目标函数，为了增强系统的鲁棒性，需要在目标函数中考虑到现在时刻的控制量u(k)对系统未来时刻的影响，采用如下目标函数：

$J=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{[y(k+j)-w(k+j)]}^2+\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\mathrm{\&lambda;}\left(j\right){\left[\mathrm{\&Delta;u}(k+j-1)\right]}^2---\left(2\right)$

式中n为最大预测长度；m称为控制长度(m≤n)，y(k+j)为(k+j)时的实际室温，w(k+j)为k+j时刻的室温设定值，Δu(k+j-1)为k+j-1时刻的模型输入与前一时刻的模型输入的差值，即Δu(k+j-1)＝u(k+j-1)-u(k+j；λ(j)是大于零的控制加权系数，取λ(j)为常数λ；

采用柔化控制，跟踪参考轨迹线，参考轨迹线如公式3所示；

w(k+j)＝a^jy(k)+(1-a^j)y_r    (3)

式中y_r、y(k)和w(k+j)分别为设定值、输出值和参考轨迹线；a为柔化系数，0＜a＜1；

目标函数中后一项的目的是压制过于剧烈的控制增量，防止系统超出限制范围或发生剧烈震荡。因此广义预测控制的问题也可以归结为求Δu(k)，Δu(k+1)，…，Δu(k+m-1)使得目标函数达到最小值的问题。

步骤402：输出预测

根据广义预测理论，引入丢番图方程：

1＝E_j(z^-1)A(z^-1)Δ+z^-JF_j(z^-1)     (4)

E_j(z^-1)B(z^-1)＝G_j(z^-1)+z^-jH_j(z^-1)    (5)

式中E_j(z^-1)＝e₀+e₁z^-1+…+e_j-1z^-j+1

$F_j\left(z^{-1}\right)={f_0}^j+{f_1}^j{z}^{-1}+...+{f_{n_a}}^j{z}^{{-n}_a}$

G_j(z^-1)＝g₀+g₁z^-1+…+g_j-1z^-j+1

$H_j\left(z^{-1}\right)={h_0}^j+{h_1}^j{z}^{-1}+...+{f_{n_b-1}}^j{z}^{{-n}_b-1}$

将公式(1)等号两侧同时乘以E_j(z^-1)Δ，并将公式(4)、(5)带入可得：

y(k+j)＝G_jΔu(k+j-1)+F_jy(k)+H_jΔu(k-1)+E_jξ(k+j)   (6)

公式(6)中的最后一项表示外界的白噪声序列，在预测未来时刻的输出时不考虑，公式(6)即可表示为：

$\hat{y}(k+j)=G_j\mathrm{\&Delta;u}(k+j-1)+F_{j}y\left(k\right)+H_j\mathrm{\&Delta;u}(k-1)---\left(7\right)$

将公式(6)写成向量形式，即为：

Y＝GΔU+Fy(k)+HΔu(k-1)+E   (8)

表示对未来j时刻的预测，将其分为已知量和未知量两部分，用f(k)表示已知量，写成矩阵形式即为：

f＝Fy(k)+HΔu(k-1)    (9)

由公式(8)得知：

$\hat{Y}=\mathrm{G\&Delta;U}+f---\left(10\right)$

式中 $\hat{Y}{[\hat{y}(k+1),\hat{y}(k+2),...,\hat{y}(k+n)]}^T;$

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

f＝[f(k+1),f(k+2),…,f(k+n)]^T；

F＝[F₁,F₂,…,F_n]^T；

H＝[H₁,H₂,…,H_n]^T；

E＝[E₁w(k+1),E₂w(k+2),…,E_nw(k+n)]；

步骤403：求取最优控制率

若令W＝[w(k+1),w(k+2),…,w(k+n)]^T

则公式(2)表示成：

J＝(Y-W)^T(Y-W)+λΔU^TΔU   (11)

用Y的最优预测值来代替Y，将公式(10)带入公式(11)，并令 $\frac{\&PartialD;J}{\&PartialD;\mathrm{\&Delta;U}}=0,$

得到：

ΔU＝(G^TG-λI)^-1G^T(W-f)    (12)

由于上述内容已经给出G、I、W和f均为矩阵，而他们的计算结果ΔU也为矩阵，为了获得最优的控制结果，防止多次迭代后发散，因此取矩阵ΔU的第一个向量，由公式10中可知即为Δu(k)，等号右边也得到的相应结果即为g^T(W-f)。

得到Δu(k)＝g^T(W-f)，即可得到公式(13)。

u(k)＝u(k-1)+g^T(W-f)    (13)

其中g^T表示(G^TG-λI)^-1G^T中的第一行向量。

与通常的最优控制不同，广义预测控制采用滚动优化策略，随时间的推移优化的目标也会发生变化。即在每一时刻都会提出一个当前时刻的局部优化目标，而不是采用不变的全局优化目标，这也说明优化过程是在线不断进行的，当预测模型存在失配、非线性扰动影响时，可以及时进行弥补和修整，保持最优的控制效果。

(23)、反馈校正

由于广义预测控制算法在滚动优化的过程中已经强调了在控制的每一时刻都会将实际值与预测值进行比较，取得局部优化目标，所以在优化过程中已经采取了反馈校正的机制，提高了系统的鲁棒性，在实际的工业过程中具有较强现实意义。

上述步骤(2)中梯度下降法的步骤如下：

首先，变风量空调的室温回路一般可以取一阶惯性加延迟的模型来表示，在离散情况下通常使用z^-1来表示时延，选取过程模型如下：

$G_m=\frac{b_m{z}^{-1}}{1-a_m{z}^{-1}}{z}^{-d}---\left(14\right)$

室温回路的模型参数中需要进行辨识的参数θ＝[a_m,b_m,d]^T，即广义预测控制中的A(z^-1)＝1-a_mz^-1，B(z^-1)＝b_mz^-1，和时间延迟d；

过程模型的输出与广义预测控制的误差为：

e(k)＝y(k)-y_m(k)＝y(k)-[a_my(k-1)-z^-db_mu(k-1)]   (15)

定义目标函数为：

$J\left(\mathrm{\&theta;}\right)=\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}\frac{1}{2}{e}^2\left(i\right)---\left(16\right)$

验证J(θ)的扶梯度方向不断修正θ(k)，直至J(θ)获得最小值，梯度下降法可以成：

$\mathrm{\&theta;}(k+1)=\mathrm{\&theta;}\left(k\right)-R\left(k\right)\&CenterDot;\mathrm{grad}\left[J\left(\mathrm{\&theta;}\right)\right]{|}_{{\mathrm{\&theta;}}_m\left(k\right)}---\left(17\right)$

式中R(k)为加权矩阵，表示J(θ)关于θ_m(k)的梯度，根据目标函数可知：

$\mathrm{grad}\left[J\left(\mathrm{\&theta;}\right)\right]{|}_{{\mathrm{\&theta;}}_m\left(k\right)}={[\frac{\&PartialD;J}{\&PartialD;a_m},\frac{\&PartialD;J}{\&PartialD;b_m},\frac{\&PartialD;J}{\&PartialD;d}]}^T\&CenterDot;e\left(k\right);$

根据公式(15)和公式(16)可知：

$\frac{\&PartialD;J}{\&PartialD;a_m}=-y(k-1)---\left(18\right)$

$\frac{\&PartialD;J}{\&PartialD;b_m}=-z^{-d}\&CenterDot;u(k-1)---\left(19\right)$

根据Z变换的实域位移定理，推导得

$\frac{\&PartialD;J}{\&PartialD;b_m}=-u(k-d_m\left(k\right)-1)$

由拉普拉斯变换与Z变换的关系z＝e^ST可知：

$\frac{\&PartialD;J}{\&PartialD;d}=z^{-d}\&CenterDot;\ln z\&CenterDot;b_m\&CenterDot;u(k-1)=z^{-d}\&CenterDot;s\&CenterDot;T\&CenterDot;b_m\&CenterDot;u(k-1)$

在上式中s为拉普拉斯变换算子，T为采样周期，采用欧拉线性变换式对上式进行化解，推导获得

$\frac{\&PartialD;J}{\&PartialD;d}==[u(k-d-1)-u(k-d-2)]\&CenterDot;b_m---\left(20\right)$

则

$\underset{\mathrm{\&theta;}}{\mathrm{grad}}\left[J\left(\mathrm{\&theta;}\right)\right]{|}_{{\mathrm{\&theta;}}_m\left(k\right)}=\left[\begin{array}{c}-y(k-1)\\ -u(k-d-1)\\ (u(k-d-1)-u(k-d-2))\&CenterDot;b_m\end{array}\right]\&CenterDot;e\left(k\right)$

根据公式(17)获得在最速下降法情况下空调系统参考对象模型特性参数向量θ_m(k)的递推公式，此时

${\mathrm{\&theta;}}_m(k+1)={\mathrm{\&theta;}}_m\left(k\right)-R\left(k\right)\&CenterDot;\left[\begin{array}{c}-y(k-1)\\ -u(k-d-1)\\ (u(k-d-1)-u(k-d-2))\&CenterDot;b_m\end{array}\right]\&CenterDot;e\left(k\right)---\left(21\right)$

令

$h\left(k\right)=\left[\begin{array}{c}-y(k-1)\\ -u(k-d-1)\\ (u(k-d-1)-u(k-d-2))\&CenterDot;b_m\end{array}\right]$

h(k)是整个辨识过程的输入数据向量，加权矩阵R(k)的作用是用来控制各输入数据分量对参数估计值的影响程度，最速下降法的收敛性能直接取决于该加权阵的选择；

选用加权矩阵：

$R^{*}\left(k\right)=\frac{1}{\mathrm{\&Sigma;}{\mathrm{\&Lambda;}}_i\left(k\right)h_i^2\left(k\right)}\mathrm{diag}[{\mathrm{\&Lambda;}}_1\left(k\right),{\mathrm{\&Lambda;}}_2\left(k\right),...,{\mathrm{\&Lambda;}}_N\left(k\right)]---\left(22\right)$

式中的Λ₁(k),Λ₂(k),...,Λ_N(k)表示在加权矩阵中对角线上的变量，在公式(22)中N＝3；

辨识出的时延参数d用整数表示，对d进行取整操作，即：

$d(k+1)=\left\{\begin{array}{c}d\left(k\right)+1(d(k+1)=d\left(k\right)\&GreaterEqual;0.5\\ d\left(k\right)|d(k+1)-d\left(k\right)|<0.5\\ d\left(k\right)-1d(k+1)-d\left(k\right)\&le;-0.5\end{array}\right.$

最终辨识出的1-a_mz^-1，b_mz^-1，d即为广义预测控制器算法中预测模型的A(z^-1)，B(z^-1)和d。

如图2所示，本发明的流程如下：

1.置初值，选取m，n，控制加权系数λ，柔化因子a等；

2.采集到设定值y_r、输出值y(k)和参考轨迹线w(k+j)；

3.使用梯度下降法来辨识模型参数，得到a_m，b_m，d作为广义预测控制的输入参数；

4.通过丢番图方程来计算E_j(z^-1)、F_j(z^-1)、G_j(z^-1)、H_j(z^-1)；

5.计算出矩阵(G^TG-λI)^-1G^T和g^T；

6.求得输入量的增量Δu(k)，即可得知下一时刻的输入量；

7.将k+1赋给k。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (5)

1.一种基于梯度下降法和广义预测控制的室温回路控制方法，其特征在于包括如下步骤：

(1)、确定预设室温的设定值；

(2)、采用梯度下降法在线辨识出变风量空调系统的室温回路的模型参数；

(3)、在变风量空调系统中加入广义预测控制器，对模型参数进行在线估计，实现自适应控制。

2.根据权利要求1所述的基于梯度下降法和广义预测控制的室温回路控制方法，其特征在于：所述步骤(3)中广义预测控制器的控制方法包括：

(21)、预测模型；

(22)、滚动优化；

(23)、反馈校正。

3.根据权利要求1所述的基于梯度下降法和广义预测控制的室温回路控制方法，其特征在于：所述步骤(21)中预测模型采用受控自回归积分滑动平均模型，模型如公式(1)所示；

A(z^-1)y(k)＝B(z^-1)u(k-d)+C(z^-1)ξ(k)/Δ   (1)

式中A(z^-1)、B(z^-1)、C(z^-1)分别是n、m和n阶的z^-1的多项式；

$A\left(z^{-1}\right)=1+a_1{z}^{-1}+...+a_{n_a}{z}^{-n_a}$

$B\left(z^{-1}\right)=b_0+b_1{z}^{-1}+...+b_{n_b}{z}^{{-n}_b}$

$C\left(z^{-1}\right)=1+c_1{z}^{-1}+...+c_{n_c}{z}^{-n_c}$

其中Δ＝1-z^-1；y(k)、u(k)、ξ(k)分别表示输出、输入和均值为零的白噪声序列，d为系统的时滞，C(z^-1)＝1。

4.根据权利要求1所述的基于梯度下降法和广义预测控制的室温回路控制方法，其特征在于：所述步骤(22)中滚动优化的步骤如下：

步骤401：设定目标函数，采用如下目标函数，

$J=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{[y(k+j)-w(k+j)]}^2+\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\mathrm{\&lambda;}\left(j\right){\left[\mathrm{\&Delta;u}(k+j-1)\right]}^2---\left(2\right)$

式中n为最大预测长度；m称为控制长度(m≤n)，y(k+j)为(k+j)时的实际室温，w(k+j)为k+j时刻的室温设定值，Δu(k+j-1)为k+j-1时刻的模型输入与前一时刻的模型输入的差值，即Δu(k+j-1)＝u(k+j-1)-u(k+j；-2λ)(j)是大于零的控制加权系数，取λ(j)为常数λ；

采用柔化控制，跟踪参考轨迹线，参考轨迹线如公式3所示；

w(k+j)＝a^jy(k)+(1-a^j)y_r   (3)

式中y_r、y(k)和w(k+j)分别为设定值、输出值和参考轨迹线；a为柔化系数，0＜a＜1；

步骤402：输出预测，

根据广义预测理论，引入丢番图方程：

1＝E_j(z^-1)A(z^-1)Δ+z^-JF_j(z^-1)   (4)

E_j(z^-1)B(z^-1)＝G_j(z^-1)+z^-jH_j(z^-1)   (5)

式中E_j(z^-1)＝e₀+e₁z^-1+…+e_j-1z^-j+1

$F_j\left(z^{-1}\right)={f_0}^j+{f_1}^j{z}^{-1}+...+{f_{n_a}}^j{z}^{{-n}_a}$

$G_j\left(z^{-1}\right)=g_0+g_1{z}^{-1}+...+g_{j-1}{z}^{-j+1}$

$H_j\left(z^{-1}\right)={h_0}^j+{h_1}^j{z}^{-1}+...+{h_{n_b-1}}^j{z}^{{-n}_b-1}$

将公式(1)等号两侧同时乘以E_j(z^-1)Δ，并将公式(4)、(5)带入可得：

y(k+j)＝G_jΔu(k+j-1)+F_jy(k)+H_jΔu(k-1)+E_jξ(k+j)   (6)

公式(6)中的最后一项表示外界的白噪声序列，在预测未来时刻的输出时不考虑，公式(6)即可表示为：

$\hat{y}(k+j)=G_j\mathrm{\&Delta;u}(k+j-1)+F_{j}y\left(k\right)+H_j\mathrm{\&Delta;u}(k-u)---\left(7\right)$

将公式(6)写成向量形式，即为：

Y＝GΔU+Fy(k)+HΔu(k-1)+E   (8)

表示对未来j时刻的预测，将其分为已知量和未知量两部分，用f(k)表示已知量，写成矩阵形式即为：

f＝Fy(k)+HΔu(k-1)   (9)

由公式(8)得知：

$\hat{Y}=\mathrm{G\&Delta;U}+f---\left(10\right)$

式中 $\hat{Y}={[\hat{y}(k+1),\hat{y}(k+2),...,\hat{y}(k+n)]}^T;$

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

f＝[f(k+1),f(k+2),…,f(k+n)]^T；

F＝[F₁,F₂,…,F_n]^T；

H＝[H₁,H₂,…,H_n]^T；

E＝[E₁w(k+1),E₂w(k+2),…,E_nw(k+n)]；

步骤403：求取最优控制率

若令W＝[w(k+1),w(k+2),…,w(k+n)]^T

则公式(2)表示成：

J＝(Y-W)^T(Y-W)+λΔU^TΔU   (11)

用Y的最优预测值来代替Y，将公式(10)带入公式(11)，并令 $\frac{\&PartialD;J}{\&PartialD;\mathrm{\&Delta;U}}=0,$

得到：

ΔU＝(G^TG-λI)^-1G^T(W-f)   (12)

每次将第一个分量加入系统，即

u(k)＝u(k-1)+g^T(W-f)   (13)

其中g^T表示(G^TG-λI)^-1G^T中的第一行向量。

5.根据权利要求1所述的基于梯度下降法和广义预测控制的室温回路控制方法，其特征在于：所述步骤(2)中梯度下降法的步骤如下：

选取过程模型如下：

$G_m=\frac{b_m{z}^{-1}}{1-a_m{z}^{-1}}{z}^{-d}---\left(14\right)$

室温回路的模型参数中需要进行辨识的参数θ＝[a_m,b_m,d]^T，即广义预测控制中的A(z^-1)＝1-a_mz^-1，B(z^-1)＝b_mz^-1，和时间延迟d；

过程模型的输出与广义预测控制的误差为：

e(k)＝y(k)-y_m(k)＝y(k)-[a_my(k-1)-z^-db_mu(k-1)]   (15)

定义目标函数为：

$J\left(\mathrm{\&theta;}\right)=\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}\frac{1}{2}{e}^2\left(i\right)---\left(16\right)$

验证J(θ)的扶梯度方向不断修正θ(k)，直至J(θ)获得最小值，梯度下降法可以成：

$\mathrm{\&theta;}(k+1)=\mathrm{\&theta;}\left(k\right)-R\left(k\right)\&CenterDot;\mathrm{grad}\left[J\left(\mathrm{\&theta;}\right)\right]{|}_{{\mathrm{\&theta;}}_m\left(k\right)}---\left(17\right)$

式中R(k)为加权矩阵，表示J(θ)关于θ_m(k)的梯度，根据目标函数可知：

$\mathrm{grad}\left[J\left(\mathrm{\&theta;}\right)\right]{|}_{{\mathrm{\&theta;}}_m\left(k\right)}={[\frac{\&PartialD;J}{\&PartialD;a_m},\frac{\&PartialD;J}{\&PartialD;b_m},\frac{\&PartialD;J}{\&PartialD;d}]}^T\&CenterDot;e\left(k\right);$

根据公式(15)和公式(16)可知：

$\frac{\&PartialD;J}{\&PartialD;a_m}=-y(k-1)---\left(18\right)$

$\frac{\&PartialD;J}{\&PartialD;b_m}=-z^{-d}\&CenterDot;u(k-1)---\left(19\right)$

根据Z变换的实域位移定理，推导得

$\frac{\&PartialD;J}{\&PartialD;b_m}=-u(k-d_m\left(k\right)-1)$

由拉普拉斯变换与Z变换的关系z＝e^ST可知：

$\frac{\&PartialD;J}{\&PartialD;d}=z^{-d}\&CenterDot;\ln z\&CenterDot;b_m\&CenterDot;u(k-1)=z^{-d}\&CenterDot;s\&CenterDot;T\&CenterDot;b_m\&CenterDot;u(k-1)$

在上式中s为拉普拉斯变换算子，T为采样周期，采用欧拉线性变换式对上式进行化解，推导获得

$\frac{\&PartialD;J}{\&PartialD;d}==[u(k-d-1)-u(k-d-2)]\&CenterDot;b_m---\left(20\right)$

则

$\underset{\mathrm{\&theta;}}{\mathrm{grad}}\left[J\left(\mathrm{\&theta;}\right)\right]{|}_{{\mathrm{\&theta;}}_m\left(k\right)}=\left[\begin{array}{c}-y(k-1)\\ -u(k-d-1)\\ (u(k-d-1)-u(k-d-2))\&CenterDot;b_m\end{array}\right]\&CenterDot;e\left(k\right)$

根据公式(17)获得在最速下降法情况下空调系统参考对象模型特性参数向量θ_m(k)的递推公式，此时

${\mathrm{\&theta;}}_m(k+1)={\mathrm{\&theta;}}_m\left(k\right)-R\left(k\right)\&CenterDot;\left[\begin{array}{c}-y(k-1)\\ -u(k-d-1)\\ (u(k-d-1)-u(k-d-2))\&CenterDot;b_m\end{array}\right]\&CenterDot;e\left(k\right)---\left(21\right)$

令

$h\left(k\right)=\left[\begin{array}{c}-y(k-1)\\ -u(k-d-1)\\ (u(k-d-1)-u(k-d-2))\&CenterDot;b_m\end{array}\right]$

h(k)是整个辨识过程的输入数据向量，加权矩阵R(k)的作用是用来控制各输入数据分量对参数估计值的影响程度，最速下降法的收敛性能直接取决于该加权阵的选择；

选用加权矩阵：

$R^{*}\left(k\right)=\frac{1}{\mathrm{\&Sigma;}{\mathrm{\&Lambda;}}_i\left(k\right)h_i^2\left(k\right)}\mathrm{diag}[{\mathrm{\&Lambda;}}_1\left(k\right),{\mathrm{\&Lambda;}}_2\left(k\right),...,{\mathrm{\&Lambda;}}_N\left(k\right)]---\left(22\right)$

式中的Λ₁(k),Λ₂(k),...,Λ_N(k)表示在加权矩阵中对角线上的变量，在公式(22)中N＝3；

辨识出的时延参数d用整数表示，对d进行取整操作，即：

$d(k+1)=\left\{\begin{array}{cc}d\left(k\right)+1& (d(k+1)=d\left(k\right)\&GreaterEqual;0.5\\ d\left(k\right)& |d(k+1)-d\left(k\right)|<0.5\\ d\left(k\right)-1& d(k+1)-d\left(k\right)\&le;-0.5\end{array}\right.$

最终辨识出的1-a_mz^-1，b_mz^-1，d即为广义预测控制器算法中预测模型的A(z^-1)，B(z^-1)和d。