基于递推最小二乘法的空调所属建筑物一阶模型实时参数辨识方法

摘要

本发明公开了一种基于递推最小二乘法的空调所属建筑物一阶模型实时参数辨识方法，实时建立负荷聚合商管辖范围内的空调负荷模型，用以进行空调负荷的实时控制与调度；将空调所属建筑物ETP模型等效为标准差分方程，同时确定输入输出序列及待辨识参数序列；通过智能电网的高级计量终端对输入输出序列进行实时数据采集，并向负荷聚合商进行数据传递，负荷聚合商利用递推最小二乘法对空调负荷模型的实时参数进行辨识。本发明弥补了国内空调负荷实时建模方面的研究不足，为其参与需求响应提供了技术支撑。

说明书

技术领域

本发明属于电力系统调度部分的参数辨识技术，尤其涉及负荷聚合商对空调负荷的 实时参数辨识。

背景技术

在硬件支持方面，负荷聚合商可以通过智能电网的高级测量终端和双向通信网络获 取用户空调的开关状态、室内外温度的实时变化情况，并对空调负荷进行实时控制与调 度。在在线辨识方法方面，目前常用的方法有最小二乘法、人工智能法、卡尔曼滤波法 等，但由于人工智能法其本身的复杂性，在实际应用中较少使用，而卡尔曼滤波法在迭 代过程中也要做大量计算工作，因此在多种参数辨识方法中，最小二乘法的原理简单， 计算量少，被广泛应用于各个领域的在线参数辨识。为了减少实时参数辨识过程中的计 算量，一种非常实用的改进的最小二乘法—递推最小二乘法被提了出来。依托于智能电 网的硬件支持和相关的在线辨识技术，负荷聚合商的空调负荷在线辨识对空调负荷参与 需求响应起着至关重要的作用。

空调负荷作为一种重要的需求响应资源，其所属建筑物ETP模型已经被广泛应用于 空调负荷控制的各个领域，但其模型中的一些参数与建筑物的墙体的厚度、窗户面积、 体积大小等因素密切相关，无法通过测量获得，因此需要通过一定的参数辨识手段对空 调负荷模型的参数进行识别。但由于我国需求响应工作起步较晚，目前有关空调负荷模 型的参数辨识还未提及，忽视了空调负荷的实时在线建模。

发明内容

发明目的：为了弥补现有需求响应过程中空调所属建筑物一阶模型参数在线识别技 术的空白，本发明提供一种基于递推最小二乘法的空调模型实时参数辨识方法，将空调 负荷的ETP模型等效为标准差分方程，确定输入输出序列及待辨识序列，通过智能电网 的高级计量终端对输入输出序列进行实时数据采集，并向负荷聚合商进行数据传递，负 荷聚合商利用递推最小二乘法对空调负荷模型的实时参数辨识，实时建立负荷聚合商管 辖范围内的空调负荷模型，用以进行空调负荷的实时控制与调度。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种基于递推最小二乘法的空调所属建筑物一阶模型实时参数辨识方法，实时建立 负荷聚合商管辖范围内的空调负荷模型，用以进行空调负荷的实时控制与调度；将空调 所属建筑物ETP模型等效为标准差分方程，同时确定输入输出序列及待辨识参数序列； 通过智能电网的高级计量终端对输入输出序列进行实时数据采集，并向负荷聚合商进行 数据传递，负荷聚合商利用递推最小二乘法对空调负荷模型的实时参数进行辨识。

上述方法具体包括如下步骤：

(1)将空调所属建筑物ETP模型等效为单输入/输出线性系统，并使用差分方程的 标准形式进行描述，确定k时刻的输入输出序列为[u(k),y(k)]，k时刻的随机变量序列 为e(k)，待辨识参数数量为2n+1，待辨识参数为a₁...a_n和b₀...b_n，负荷聚合商管辖范围 内的空调负荷模型为：

$y\left(k\right)=-\underset{i=1}{\overset{n}{\Sigma }}{a}_{i}y(k-i)+\underset{i=0}{\overset{n}{\Sigma }}{b}_{i}u(k-i)+e\left(k\right)---\left(1\right)$

式中：n为差分方程的阶数；

(2)负荷聚合商通过智能电网的高级计量终端采集输入输出数据，并对采集的输 入输出数据进行预处理形成输入输出序列；

(3)获取k-2n-3～k-3时刻的N＝2n+1个输入输出序列，根据式(1)列出N个 差分方程为：

y(n+1)＝-a₁y(n)-…a_ny(1)+b₀u(n+1)+…+b_nu(1)+e(n+1)

y(n+2)＝-a₁y(n+1)-…a_ny(2)+b₀u(n+2)+…+b_nu(2)+e(n+2)(2)

......

y(n+N)＝-a₁y(n+N-1)-…a_ny(N)+b₀u(n+N)+…+b_nu(N)+e(n+N)

将式(2)写成矩阵的形式为：

Y(N)＝Φ(N)η(N)+E(N)(3)

式中：

Y(N)＝[y(n+1),y(n+2),…,y(n+N)]^T(4)

η(N)＝[a₁,…,a_n,b₀,…,b_n]^T(5)

E(N)＝[e(n+1),e(n+2),…,e(n+N)]^T(6)

式中：Y(N)为输出矩阵，Φ(N)为由输入输出序列组成的矩阵，η(N)为待辨识参 数矩阵，E(N)为误差矩阵；

(4)对采集的数据进行标准最小二乘法的参数估计，获取递推最小二乘法递推过 程的初始值：

P(N)＝[Φ^T(N)Φ(N)]^-1(9)

$\hat{\eta }\left(N\right)={\left[{\Phi }^T\left(N\right)\Phi \left(N\right)\right]}^{-1}{\Phi }^T\left(N\right)Y\left(N\right)---\left(10\right)$

式中：P(N)为中间变量，为经过标准最小二乘法求得的待辨识参数估计值；

(5)负荷聚合商提取k-2时刻的输入输出数据，并对输入输出数据进行预处理得 到k-1时刻的输入输出序列为[u(n+N+1),y(n+N+1)]；

(6)负荷聚合商通过k-2时刻的输入输出序列[u(n+N+1),y(n+N+1)]和P(N)， 对修正系数进行计算：

(7)根据计算更新后的待辨识参数估计值

根据更新式(1)，实时更新空调负荷模型；

(8)计算k-2时刻的P(N+1)：

然后令P(N)＝P(N+1)，k＝k+1，返回步骤(5)。

优选的，所述步骤(2)和步骤(5)中，对采集的输入输出数据进行预处理，具体 包括以下内容：

(a)对于室内/室外温度，按照下述步骤进行预处理：

(a1)由于室内/室外温度变化比较缓慢，根据统计规律，连续采集的五个点数值变 化应该不大，因此按照如下方法计算室内/室外温度的平均值和标准差S：

$\bar{x}=\frac{1}{5}(x_{i-2}+x_{i-1}+x_i+x_{i+1}+x_{i+2})---\left(15\right)$

$S=\sqrt{\frac{1}{5}\underset{i=-2}{\overset{2}{\Sigma }}{(x_i-\bar{x})}^2}---\left(16\right)$

式中：x_i表示连续采集的第i个室内/室外温度；

(a2)判断是否成立：若不成立，则说明x_i是正常数据；否则说明x_i是 坏数据，需要进行修正；

(a3)坏数据修正公式为：

$\widetilde{x_i}=\frac{\alpha }{4}(x_{i-2}+x_{i+2})+\frac{\beta }{4}(x_{i-1}+x_{i+1})---\left(17\right)$

式中：表示修正后的第i个室内/室外温度，α，β为自定义的修正系数，且 α+β＝2；

(b)对于空调负荷的开关状态量，按照下述步骤进行预处理：

在室内温度上升或者下降过程中，空调的开关状态保持不变，只有在到达所设温度 的上下限的时候空调开关才会发生变化；设连续采集的三个开关状态量为{x_i-1,x_i,x_i+1}， 对应采集的三个室内温度为{y_i-1,y_i,y_i+1}，进行如下判断：

①若x_i-1·x_i·x_i+1＝1，则表示x_i是正常数据；

②若x_i-1·x_i+1＝1且x_i-1·x_i·x_i+1＝0，则表示x_i是坏数据，需要进行修正，令x_i＝x_i+1；

③若x_i-1·x_i+1＝0，则：

在(y_i-y_i-1)(y_i+1-y_i)>0时：若x_i-1·x_i＝1，则表示x_i是正常数据；若x_i-1·x_i＝0，则 表示x_i是坏数据，需要进行修正，令x_i＝x_i-1；

在(y_i-y_i-1)(y_i+1-y_i)<0时：若若x_i·x_i+1＝1，则表示x_i是正常数据；若x_i·x_i+1＝0， 则表示x_i是坏数据，需要进行修正，令x_i＝x_i+1；

式中：·表示逻辑与计算。

有益效果：本发明提供的基于递推最小二乘法的空调所属建筑物一阶模型实时参数 辨识方法，可实时构建空调负荷的等效热参数模型，能够为相关部门提供空调负荷的相 关参数，为空调负荷参与电力系统的削峰、调频、备用、抑制可再生能源波动等提供了 依据；同时本发明弥补了国内空调负荷实时建模方面的研究不足，为其参与需求响应提 供了技术支撑。

附图说明

图1为本发明方法的总流程图；

图2为负荷聚合商系统的硬件组成图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

一种基于递推最小二乘法的空调所属建筑物一阶模型实时参数辨识方法，具体实现 流程如图1所示，下面就各个步骤加以具体说明。

(一)将空调所属建筑物ETP模型等效为单输入/输出线性系统，并使用差分方程 的标准形式进行描述，确定k时刻的输入输出序列为[u(k),y(k)]，k时刻的随机变量序 列为e(k)，待辨识参数数量为2n+1，待辨识参数为a₁...a_n和b₀...b_n，负荷聚合商管辖范 围内的空调负荷模型为：

$y\left(k\right)=-\underset{i=1}{\overset{n}{\Sigma }}{a}_{i}y(k-i)+\underset{i=0}{\overset{n}{\Sigma }}{b}_{i}u(k-i)+e\left(k\right)---\left(1\right)$

式中：n为差分方程的阶数。

式(1)的推导过程为如下：

(1.1)空调所属建筑物ETP模型的离散方程为：

s＝0时：$T_i^{t+\Delta t}=T_o^{t+\Delta t}-(T_o^{t+\Delta t}-T_i^t)ϵ---(1-1)$

s＝1时：$T_i^{t+\Delta t}=T_o^{t+\Delta t}+QR-(T_o^{t+\Delta t}+QR-T_i^t)ϵ---(1-2)$

式中：为t时刻室内温度(单位：℃)，为t时刻室外温度(单位：℃)，△t表 示采样时间间隔，ε为散热系数，T_c为时间常数；Q为制冷/制热量(单位： W)，R为等效热阻抗(单位：℃/W)；s为开关状态量，s＝0表示关，s＝1表示开。

(1.2)根据式(1-1)和(1-2)，将空调负荷的差分方程进一步改写为：

$T_i^{t+1}={\mathrm{ϵT}}_i^t+(1-ϵ)T_o^{t+1}+QR(1-ϵ)S\left(t\right)---(1-3)$

式中：S(t)表示t时刻的开关状态量，S(t)＝0表示关，S(t)＝1表示开。

以n＝1的情况为例进行推导，对式(1-3)进行如下的变量替换：

$y\left(k\right)=T_i^{t+1}---(1-4)$

$y(k-1)=T_i^t---(1-5)$

$u\left(k\right)=T_o^{t+1}---(1-6)$

u(k-1)＝S(t)(1-7)

a₁＝-ε(1-8)

b₀＝1-ε(1-9)

b₁＝QR(1-ε)(1-10)

最终得到：

y(k)＝-a₁y(k-1)+b₀u(k)+b₁u(k-1)(1-11)

(二)负荷聚合商通过智能电网的高级计量终端采集输入输出数据，并对采集的输 入输出数据进行预处理形成输入输出序列。

(三)获取k-2n-3～k-3时刻的N＝2n+1个输入输出序列，根据式(1)列出N 个差分方程为：

y(n+1)＝-a₁y(n)-…a_ny(1)+b₀u(n+1)+…+b_nu(1)+e(n+1)

y(n+2)＝-a₁y(n+1)-…a_ny(2)+b₀u(n+2)+…+b_nu(2)+e(n+2)(2)

......

y(n+N)＝-a₁y(n+N-1)-…a_ny(N)+b₀u(n+N)+…+b_nu(N)+e(n+N)

将式(2)写成矩阵的形式为：

Y(N)＝Φ(N)η(N)+E(N)(3)

式中：

Y(N)＝[y(n+1),y(n+2),…,y(n+N)]^T(4)

η(N)＝[a₁,…,a_n,b₀,…,b_n]^T(5)

E(N)＝[e(n+1),e(n+2),…,e(n+N)]^T(6)

式中：Y(N)为输出矩阵，Φ(N)为由输入输出序列组成的矩阵，η(N)为待辨识参 数矩阵，E(N)为误差矩阵。

(四)对采集的数据进行标准最小二乘法的参数估计，获取递推最小二乘法递推过 程的初始值：

P(N)＝[Φ^T(N)Φ(N)]^-1(9)

$\hat{\eta }\left(N\right)={\left[{\Phi }^T\left(N\right)\Phi \left(N\right)\right]}^{-1}{\Phi }^T\left(N\right)Y\left(N\right)---\left(10\right)$

式中：P(N)为中间变量，为经过标准最小二乘法求得的待辨识参数估计值。

式(10)的推导过程如下：

最小二乘法的原理是根据得到的标准差分方程，将待辨识参数η的估计值带入， 使其方程的残差平方和最小，因此有：

$J=\underset{i=1}{\overset{N}{\Sigma }}{e}^2(n+i)=E^T\left(N\right)E\left(N\right)---(4-1)$

将E(N)＝Y(N)-Φ(N)η(N)带入式(4-1)，并使J对求偏导为0，即：

$\frac{\partial J}{\partial \eta }{|}_{\eta =\hat{\eta }}=0---(4-2)$

经整理后得到：

${\Phi }^T\left(N\right)\Phi \left(N\right)\hat{\eta }\left(N\right)={\Phi }^T\left(N\right)Y\left(N\right)---\left(4-3\right)$

进一步可得待辨识参数估计值的计算公式为式(10)。

(五)负荷聚合商提取k-2时刻的输入输出数据，并对输入输出数据进行预处理得 到k-1时刻的输入输出序列为[u(n+N+1),y(n+N+1)]。

(六)负荷聚合商通过k-2时刻的输入输出序列[u(n+N+1),y(n+N+1)]和 P(N)，对修正系数进行计算：

(七)根据计算更新后的待辨识参数估计值

根据更新式(1)，实时更新空调负荷模型。

以步骤(一)中列举的n＝1的情况进行说明，实时更新空调负荷模型的方法，根据 待辨识参数估计值可得参数a₁，b₀，b₁的辨识结果且存在如下关 系：

$\widehat{a_1}=-\widehat{ϵ}---\left(7-1\right)$

$\widehat{b_0}=1-\widehat{ϵ}---\left(7-2\right)$

$\widehat{b_1}=QR\left(1-\widehat{ϵ}\right)---\left(7-3\right)$

为了得到相对准确的辨识结果，用式(7-1)减式(7-2)，得到：

$\widehat{ϵ}=\frac{1-\left(\widehat{a_1}+\widehat{b_0}\right)}{2}---\left(7-4\right)$

将式(7-4)带入式(7-3)，可得：

$QR=\frac{2{\hat{b}}_1}{1+\left(\widehat{a_1}+\widehat{b_0}\right)}---(7-5)$

(八)计算k-2时刻的P(N+1)：

然后令P(N)＝P(N+1)，k＝k+1，返回步骤(5)。

所述步骤(二)和步骤(五)中，对采集的输入输出数据进行预处理，具体包括以 下内容：

(a)对于室内/室外温度，按照下述步骤进行预处理：

(a1)按照如下方法计算室内/室外温度的平均值和标准差S：

$\bar{x}=\frac{1}{5}(x_{i-2}+x_{i-1}+x_i+x_{i+1}+x_{i+2})---\left(15\right)$

$S=\sqrt{\frac{1}{5}\underset{i=-2}{\overset{2}{\Sigma }}{(x_i-\bar{x})}^2}---\left(16\right)$

式中：x_i表示连续采集的第i个室内/室外温度；

(a2)判断是否成立：若不成立，则说明x_i是正常数据；否则说明x_i是 坏数据，需要进行修正；

(a3)坏数据修正公式为：

$\widetilde{x_i}=\frac{\alpha }{4}(x_{i-2}+x_{i+2})+\frac{\beta }{4}(x_{i-1}+x_{i+1})---\left(17\right)$

式中：表示修正后的第i个室内/室外温度，α，β为自定义的修正系数，且 α+β＝2；

(b)对于空调负荷的开关状态量，按照下述步骤进行预处理：

设连续采集的三个开关状态量为{x_i-1,x_i,x_i+1}，对应采集的三个室内温度为 {y_i-1,y_i,y_i+1}，进行如下判断：

①若x_i-1·x_i·x_i+1＝1，则表示x_i是正常数据；

②若x_i-1·x_i+1＝1且x_i-1·x_i·x_i+1＝0，则表示x_i是坏数据，需要进行修正，令x_i＝x_i+1；

③若x_i-1·x_i+1＝0，则：

在(y_i-y_i-1)(y_i+1-y_i)>0时：若x_i-1·x_i＝1，则表示x_i是正常数据；若x_i-1·x_i＝0，则 表示x_i是坏数据，需要进行修正，令x_i＝x_i-1；

在(y_i-y_i-1)(y_i+1-y_i)<0时：若若x_i·x_i+1＝1，则表示x_i是正常数据；若x_i·x_i+1＝0， 则表示x_i是坏数据，需要进行修正，令x_i＝x_i+1；

式中：·表示逻辑与计算。

所述步骤(六)中的K(N)、步骤(七)中的步骤(八)中的P(N+1)推 导过程如下：

当增加一组输入输出序列为[u(n+N+1),y(n+N+1)]时，式(4)和式(7)可更 新为：

$Y(N+1)=\left(\begin{array}{c}Y\left(N\right)\\ y^T(N+1)\end{array}\right)---(8-1)$

则根据N+1组输入输出序列可得到更新后的待辨识参数估计值

$\hat{\eta }(N+1)={\left({\Phi }^T\right(N+1\left)\Phi \right(N+1\left)\right)}^{-1}{\Phi }^T(N+1)Y(N+1)---(8-3)$

将式(8-1)和(8-2)带入(8-3)可得：

令P(N)＝[Φ^T(N)Φ(N)]^-1，则：

根据矩阵的下列性质：

(A+BCD)^-1＝A^-1-A^-1B(C^-1+DA^-1B)^-1DA^-1(8-6)

可得：

将式(8-6)带入(8-4)可得递推公式(11)、(13)和(14)。

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