办公建筑在全生命周期中普遍存在能耗大、能效低等问题，导致严重的能源浪费。其节能潜力巨大，降耗幅度可以达到30%~50%^[1]。其中，办公建筑运行能耗占比最大。因此，研究办公建筑的运行能耗情况具有重要意义，对其进行实时准确地预测可以为优化运行效率提供数据决策，从而实现节能目标。

建筑能耗的预测方法主要分为两大类：正演模型和数据驱动模型。笔者根据获得的办公建筑历史能耗数据特点，利用数据驱动模型中的时间序列法进行能耗短期预测。机器学习预测方法是最常用的短期能耗预测方法，已经广泛应用于风速预测^[2]、电力需求预测^[3]、建筑能耗预测^[4]、建筑冷负荷预测^[5]等能源领域。目前国内外建筑能耗时间序列预测方法应用较为广泛的有：人工神经网络^[6](ANN, artificial neural network)、差分整合移动平均自回归^[7](ARIMA, autoregressive integrated moving average)、支持向量回归^[8](SVR, support vector regression)、多元线性回归^[9](MLR, multiple linear regression)等。办公建筑除了室外气象参数和围护结构等本体能耗外，还有机电设备、节假日、人员占用率等复杂的运行能耗，造成其能耗时间序列具有非线性特征。为了深入挖掘其时间序列内部的非线性规律，提高预测精度，研究者们将不同的算法与人工神经网络和支持向量机结合，建立了多种办公建筑混合预测模型。文献[9]根据天气、人员等多种因素利用MLR方法对办公建筑的逐日冷负荷进行了回归预测，与实际负荷相比平均绝对百分误差小于8%。文献[10]利用反馈神经网络(BPNN, back propagation neural network)，径向基函数神经网络(RBFNN, radial basis function neural network)，广义回归神经网络(GRNN, generalized regression neural network)，支持向量机(SVM, support vector machine)方法分别对办公建筑的逐时能耗进行分析，以中国广州某办公楼为实例，取得了较好的预测效果。文献[11]建立了小波支持向量机(wavelet-SVM)，小波偏最小二乘支持向量机(Wavelet-PLS-SVM, wavelet-partial least squares-support vector machine)模型对办公建筑的逐时能耗进行了分析，最终得出1 h，2 h，3 h和24 h的预测结果。在其他应用领域^[12]将平均速度、平均占有率和平均交通流时间序列利用贝叶斯理论集成到一个时间序列中，将混沌理论和SVR相结合用于交通流的时间序列预测，达到了较高的预测精度。上述研究都是基于多变量输入的，相较于单变量时间序列需要考虑的因素复杂，文献[13]使用历史时刻的办公建筑冷负荷根据贝叶斯理论映射到同一时间维度作为机器学习预测模型的单变量输入，利用Chaos-SVR和WD-SVR(wavelet decomposition-support vector regression)分别对建筑冷负荷时间序列进行预测，结果表明，Chaos-SVR效果更优。因此，选择恰当的预测模型不仅可以避免单变量输入的局限性，而且会取得较好的预测效果，但由于组合模型在迭代过程中存在误差累积的问题，导致不能很好地满足实际优化运营的需求。

为了解决以上问题，笔者探究了混沌现象背后隐藏的决定性规律使得将混沌理论应用于办公建筑运行能耗预测领域的可能性，展开了以下研究工作：

1) 分析办公建筑运行能耗的特点及构成，进行数据预处理和时间序列的混沌特性判别。

2) 将混沌理论(Chaos theory)和支持向量理论结合，构建Chaos-SVR办公建筑能耗预测模型。

3) 利用Markov链修正Chaos-SVR组合模型中由于参数传递存在的累积误差。

4) 以西安某办公建筑监测数据为例进行实例分析，并与其他预测方法进行比较，验证了模型的有效性。

文中主要研究办公建筑运行能耗的混沌特性，其他大型公共建筑比如商业建筑、医疗建筑等能耗是否具有混沌特性，模型的可推广性将是下一步的研究重点。

1 能耗数据特征分析

办公建筑能耗混沌时间序列预测研究的前提是确定该时间序列具有混沌特征^[14]。文中通过西安市建筑能耗监测平台，采集某办公建筑2019年8月1日0点到9月30日23点各项分类能耗数据，如表 1所示。

从表 1可以看出，办公建筑运行能耗包括水耗量、电耗量、燃气量、集中供冷/热量和其他能源使用量。为了对办公建筑运行能耗时间序列进行混沌特性分析，由于其组成结构复杂，需要对多项数据统一量纲。因此，将所有运行能耗的量纲转化为千克标准煤(kgce)的形式，从而得到该办公建筑运行能耗的逐时序列。然后对该逐时序列数据进行处理，利用最大Lyapunov指数法来判断其是否具有混沌特性。最大Lyapunov指数法是基于相空间重构的，因此对时间序列进行相空间重构是研究的前提。

1.1 相空间重构

相空间是一种工具，用于描述由确定自由度的确定性混沌系统生成的单变量或多变量时间序列重构的动态系统^[15]。相空间重构的主要原理是在拓扑等价意义下从高维空间中恢复混沌吸引子的动力学特性^[15]，对于总能耗的逐时时间序列x₁, x₂, x₃, …, x_n－1, x_n，适当选取嵌入维数m和延迟时间τ，得到重构相空间

$ Y_{m}(i)=[x(i), x(i+\tau), x(i+2 \tau), \cdots, x(i+(m-1) \tau)], $

(1)

式中，i=1, 2, …, n; N=n－(m－1)τ。N表示由时间序列构成m维后，由n个相点变成了N个相点。

嵌入维数和延迟时间的合理选取决定了办公建筑能耗时间序列相空间重构的质量，同时也保证了在实际预测中有足够高的预测精度。研究表明^[16]，将τ和m联合起来的嵌入窗宽的确定也会影响重构相空间的质量。C-C^[17]法通过序列的关联积分来构成统计量，通过统计量和延迟时间的关系图同时计算出延迟时间τ和最佳窗宽^[16]τ_w，再根据τ_ω=(m－1)τ求出嵌入维数。该方法既能有效地减少互信息量法的计算量，又能保持时间序列的非线性特征。故本文采用C-C法求取延迟时间和嵌入维数，具体步骤如下：

1) 将办公建筑能耗时间序列{x_i}, i=1, 2, …, n分成t个不相交的时间序列，当t=1时为时间序列本身，当t为一般的自然数时，时间序列长度l=Nt。构造办公建筑能耗时间序列的关联积分函数为

$ C(m, n, r, t)=\frac{2}{N(N-1)} \sum\limits_{1 \leqslant i \leqslant j \leqslant N} \theta\left(r-d_{i j}\right), r>0, $

(2)

式中，d_ij=‖X_i－X_j‖，X_i为相空间的点，θ(·)为阶跃函数：

$ \theta(\cdot)= \begin{cases}0, & x<0, \\ 1, & x \geqslant 0。\end{cases} $

(3)

2) 统计量代表了非线性时间序列的相关性。通过时间序列的关联积分来构成每个不相交的时间子序列的检验统计量：

$ S(m, n, r, t)=\frac{1}{t} \sum\limits_{s=1}^{t}\left[C_{s}(m, N / t, r, t)-C_{s}^{m}(1, N / t, r, t)\right], $

(4)

令n→∞，则：

$ S(m, r, t)=\frac{1}{t} \sum\limits_{s=1}^{t}\left[C_{s}(m, r, t)-C_{s}^{m}(1, r, t)\right],(m=2,3, \cdots) 。$

(5)

3) 如果时间序列是独立同分布的，对于固定的m和t，当N→∞时，r取任意值时，S(m, r, t)恒等于零。但实际的时间序列是有限的，因此一般S(m, r, t)不等于零。那么，局部最大时间间隔可以取S(m, r, t)的零点或对所有的半径r相互差别最小的时间点，这意味着这些点几乎是均匀分布的。选择对应值最大和最小的2个半径r，定义关于r的最大偏差：

$ \Delta S(m, t)=\max \left\{S\left(m, r_{j}, t\right)\right\}-\min \left\{S\left(m, r_{j}, t\right)\right\}{\text {。}} $

(6)

局部最大时间t应该是S(m, r, t)的零点和ΔS(m, t)的最小值。但是S(m, r, t)的零点对所有m, r应几乎相等；ΔS(m, t)的最小值对所有m应几乎相等。延迟时间τ对应着这些局部最大时间t中的第一个。

4) 通过统计量和延迟时间的关系图来确定τ和嵌入窗宽τ_ω，从而确定出嵌入维数^[16]。

几种重要渐进分布的数学统计结果表明：当2≤m≤5, σ/2≤r≤2σ, N≥500时(σ指时间序列的均方差或标准差)，渐近分布可以通过有限序列很好的近似，S(m, n, r, 1)能代表序列的相关性。具体计算时，根据办公建筑能耗时间序列特性，取m=2, 3, 4, 5, r_i=iσ/2, i=1, 2, 3, 4, 相点数n=3 000，计算下列3个检验统计量^[17]：

$ \bar{S}(t)=\frac{1}{16} \sum\limits_{m=2}^{5} \sum\limits_{j=1}^{4} S\left(m, r_{j}, t\right) , $

(7)

$ \Delta \bar{S}(t)=\frac{1}{4} \sum\limits_{m=2}^{5} \Delta S(m, t), $

(8)

$ S_{\text {cor }}(t)=\Delta \bar{S}(t)+|\bar{S}(t)|。$

(9)

根据上述计算结果在MATLAB中绘制图像，得到C-C法重构办公建筑能耗时间序列统计量曲线如图 1所示。

分析得到，S(t)第一个零点或ΔS(t)第一次取得的极小值，时间序列独立的第一个局部最大值的对应的延迟时间即为最佳延迟时间τ=2。S_cor(t)取得最小值，时间序列独立的第一个整体最大时间窗口t为最佳窗宽τ_w=2，根据嵌入时间窗宽公式τ_w=(m－1)τ可确定嵌入维数m=2。

因此，该办公建筑总能耗时间序列重构的相空间为Y₂(i)=[x(i), x(i+2)]。

1.2 混沌特性判别

一个系统的Lyapunov指数谱中有一个或一个以上的正值就可以肯定其存在混沌特性^[18]。在相空间重构后，相空间中的两条相邻轨线将随时间逐渐发散或聚合，而Lyapunov指数是其轨线的收敛或发散率，正的最大Lyapunov指数反映时间序列具有混纯特征，负的最大Lyapunov指数表明时间序列具有随机性或周期性。最大Lyapunov指数越大，时间序列的非线性越强，对初值越敏感^[18]。

小数据量法^[19]可以充分利用所有能够利用的数据，对小数据组来说比较可靠，计算量不大，相对容易操作，计算结果精度高。因此采用小数据量方法计算最大Lyapunov指数。小数据量法求取最大Lyapunov指数的步骤如下：

1) 平均周期是从相空间的意义上度量一个吸引子全局属性的最小单元，利用快速傅里叶变换得到频谱信息，计算出能耗的平均频率根据其倒数估算出时间序列的平均周期τ_p。

经过傅里叶变换后，平均周期可表示为

$ \tau_{p}=\sum\limits_{i=1}^{L} P_{i}\left(\sum\limits_{i=1}^{L} f_{i} P_{i}\right)-1, $

(10)

式中：f_i(i=1, 2, …L)是0 Hz和奈奎斯特频率之间均匀分布的频率点；P_i是在频率点f_i处的幅值。用C-C法计算出延迟时间τ和嵌入维数m。对于办公建筑能耗时间序列{x_i}，重构后的相空间为式(1)。

2) 寻找重构相空间中每个点Y_m(t)的最近邻点Y_m($\hat t$)，并限制短暂分离：

$ L_{t}(0)=\min _{\hat{t}}\left\|Y_{m}(t)-Y_{m}(\hat{t})\right\|,|t-\hat{t}|>\tau_{p}, $

(11)

式中：$\hat t$=1, 2, …, N，t≠$\hat t$；τ_p为办公建筑能耗时间序列的平均周期；L_t(0)是初始时刻一对最近邻点之间的距离。

3) 根据相空间中每个点的最邻近点的平均发散速率估计最大Lyapunov指数为

$ \lambda(i)=\frac{1}{i \varDelta t} \frac{1}{(N-i)} \sum\limits_{j=1}^{N-i} \ln \frac{L_{t}(i)}{L_{t}(0)}, $

(12)

式中：Δt为样本周期。对相空间中的每个点Y_m(t)，计算出该邻域点经过i个离散事件步长后的距离

$ L_{t}(i)=\left\|Y_{m}(t+i)-Y_{m}(\hat{t}+i)\right\|, i=1,2, \cdots, \min (N-t, N-\hat{t}), $

(13)

估计得到：

$ L_{t}(i)=C_{j} \mathrm{e}^{\lambda(i \Delta t)}, C_{j}=L_{t}(0) \text { 。} $

(14)

最大Lyapunov指数的几何意义是量化初始闭轨道的指数发散和估计系统的总体混沌水平的量^[11]，对式(14)两边取对数得到

$ \ln L_{t}(i)=\ln C_{j}+\lambda(i \Delta t)(i=1,2, \cdots N)。$

(15)

4) 对每个i求出所有t的lnL_t(i)平均值y(i)，即：

$ y(i)=\frac{1}{q \Delta t} \sum\limits_{j=1}^{q} \ln L_{t}(i)。$

(16)

式中：q是非零L_t(i)的数目。最大Lyapunov指数相当于式(15)直线的斜率，可以通过最小二乘法逼近这条直线得到。

根据相空间重构C-C法求出的最佳延迟时间τ=2及嵌入维数m=2，利用小数据量法计算得出如图 2所示办公建筑能耗时间序列的最小二乘法拟合直线。

根据直线斜率可以得到最大Lyapunov指数值为0.007 7，此为正值，由此可得出办公建筑总能耗量时间序列具有混沌特性。

2 Chaos-SVR组合模型预测

办公建筑的Chaos-SVR的预测方法主要运用混沌理论，构造非线性映射利用相空间重构技术将原始时间序列数据所隐含的信息进行还原重构，并将其作为SVR模型的输入，进行非线性训练和预测。基于办公建筑能耗的Chaos-SVR时间序列模型建立及预测过程如下：

1) 构建预测模型结构，主要包括3部分：输入向量、支持向量的非线性变化和输出。西安某办公建筑运行能耗的时间序列进行重构后的相空间为Y₂(i)=[x(i), x(i+2)]，相点数为N=n-2，每个相点包含了每个序列的主要特征，能够逼近建筑能耗的真实状态。向量x_i=(x_i, l)是模型的输入向量，β=(β₁, β₂, …β_l)是输出的权重向量，y_i=f(x_i)是模型的输出。Chaos-SVR模型结构如图 3所示。

2) 利用支持向量机中的ε不敏感损失函数，得到ε-SVR模型，采用高斯径向基函数f(x_i, x)=exp(－γ‖x_i－x‖²)作为ε-SVR的核。ε-SVR估计函数的形式为：${y_i} = f\left( {{x_i}} \right) = \sum\limits_{i = 1}^l {\left( {{\alpha _i} - {\alpha _i}^*} \right)} f\left( {{x_i}, x} \right) + b$，其中α_i, α_i^*是lagrange乘子。对于ε-SVR，一个重要的问题是模型参数的选择，采用遗传算法(genetic algorithm，GA)来优化惩罚因子C和核函数参数γ。

3) 进行网络循环学习，修正SVR的输出及权值，直到误差控制在允许范围或迭代次数达到上限，训练结束。

4) 取样本前n个数据之后连续的Q个数据点作为预测数据。对模型进行预测，网络的输出值即为预测值，将预测值进行相应的反归一化处理，得到预测的实际值。

为了验证模型的有效性，选取2019年9月4日~9月24日共计20天480组数据作为训练数据，以9月25日和26日2天共计48组数据作为预测数据来进行模型验证。并与目前在时间序列预测领域应用广泛的2种机器学习方法——非线性自回归神经网络(nonlinear auto regressive neural network)和支持向量回归(SVR)方法的预测值和实际值进行比较，利用均方根误差RMSE(root mean squared error)平均绝对百分比误差MAPE(mean absolute percent error)作为评价指标。预测对比结果如图 4所示。

预测结果误差对比如图 5所示。

预测结果误差如表 2所示。

从图 4中可以看出，Narnet预测曲线波动最大，SVR预测曲线较为平缓，Chaos-SVR预测曲线距离真实值曲线最接近，部分值几乎完全拟合。由图 5得，Narnet预测方法误差最大，SVR预测方法误差较大，Chaos-SVR预测方法误差在0值附近小范围浮动。由表 2可以看出，Chaos-SVR预测模型的RMSE为6.326 1，MAPE为0.287 5，相较于其他2种方法均最小。实验是在同样的环境和同一组数据的基础上进行的，Chaos-SVR预测方法相较于Narnet和SVR方法，预测效果最好，预测精度明显提高。但Chaos-SVR组合模型在预测的过程中，还有部分值存在较大的相对误差，为了更加精确地进行办公建筑运行能耗预测，需要对这部分的误差进行修正。导致这种误差存在的原因可能是由于组合模型在迭代过程中，由于参数传递而出现的过程误差积累，这样的误差积累会造成预测效果存在一定的缺陷。为了进一步提高预测精度，文中利用Markov链进行修正。

3 Markov链修正误差

Markov链是一种具有无后效性的随机过程^[20]，Markov链适用于对存在波动性的预测问题进行修正描述^[21]。它可以根据某一已知条件的时刻条件下推求出下一时刻的概率分布，得到下一时刻的状态而与其他时刻无关。因此，在传统Chaos-SVR预测模型基础之上利用Markov链对累积误差进行修正，从而实现办公建筑能耗时间序列的精确预测。

1) 误差状态区间^[22]划分。

误差状态区间表示未来时刻预测值的误差状态，用E₁, E₂, …, E_m分别来表示，并且状态转移只在t₁, t₂, …, t_m等可数时刻发生。定义误差均值X和标准差S分别为

$ \overline{X}=\frac{1}{n} \sum\limits_{k=1}^{n} \overline{\boldsymbol{A}}(t), $

(17)

$ S=\sqrt{\frac{1}{n-1} \sum\limits_{k=1}^{n}(\overline{\boldsymbol{A}}(t)-\overline{X})}, $

(18)

式中：A(t)为Chaos-SVR预测模型所求出的预测值矩阵，A(t)为真实值矩阵，根据中心极限原理^[15]，将误差状态区间划分为5组，即(－∞, X－S)，[X－S, X－0.5S)，[X－$\overline {0.5S} $, X+0.5S)，[X+0.5S, X+S)，[X+S, ∞]。

2) 状态转移概率矩阵。

Markov链由状态E_i经过k步转移到状态E_j的转移概率用P_ij^(k)表示，其含义是指未来时刻会处于哪一个误差状态区间，为

$ P_{i j}^{(k)}=\frac{n_{i j}^{(k)}}{N_{i}}, $

(19)

式中：N_i表示状态E_i出现的总次数；n_ij^(k)表示状态E_i误差经过k步转移到状态E_j的次数；n为划分的状态区间数目。那么一步状态转移概率矩阵为

$ \boldsymbol{P}^{(1)}=\left[\begin{array}{cccc} P_{11} & P_{12} & \cdots & P_{1 m} \\ p_{21} & p_{22} & \cdots & p_{2 m} \\ \vdots & \vdots & & \vdots \\ p_{m 1} & p_{m 2} & \cdots & p_{m m} \end{array}\right], $

(20)

使用C-K方程(切普曼柯尔莫哥洛夫方程)计算k步的转移概率矩阵为

$ \boldsymbol{P}(k)=(\boldsymbol{P}(1))^{k}。$

(21)

3) 预测模型。

设初始向量为P(0)，经k步转移后的状态向量为

$ \boldsymbol{P}(k)=\boldsymbol{P}(0) \times \boldsymbol{P}(k)=\boldsymbol{P}(0) \times(\boldsymbol{P}(1))^{k}，$

(22)

据此可以判断k步所处的状态空间。

4) 根据步骤2)所得概率转移矩阵和步骤3)所得状态向量，得出误差修正后的值。

以该办公建筑9月25日的预测数据为例进行Markov链参数的设置。首先通过预测能耗值与实际能耗值的误差均值X=0.008和标准差S=0.088确定误差状态空间为：E₁=(-∞, -0.08)，E₂=[-0.08, -0.036)，E₃=[-0.036, 0.052)，E₄=[0.052, 0.096)，E₅=[0.096, +∞)。

根据划分的状态区间对Chaos-SVR模型预测值的误差值进行状态划分，划分结果如表 3所示。

根据上述误差状态区间的确定可以得到Chaos-SVR预测结果的Markov链的状态转移情况，进一步计算出状态转移概率矩阵P为

$ \boldsymbol{P}=\left[\begin{array}{ccccc} 0 & 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{1} & 0 \\ \frac{2}{16} & \frac{1}{16} & \frac{8}{16} & \frac{1}{16} & \frac{2}{16} \\ 0 & 0 & \frac{3}{3} & 0 & 0 \\ 0 & 0 & \frac{2}{2} & 0 & 0 \end{array}\right] \text { 。} $

(23)

根据Markov链预测模型得到办公建筑9月25日的能耗预测状态向量，针对Chao-SVR预测模型存在的累积误差进行修正，得到该办公建筑能耗的最终预测值。为了验证基于Markov链修正后的Chaos-SVR预测效果，对修正前后的预测结果进行比较。对比结果如图 6所示。

Markov链修正前后预测结果误差对比如图 7所示。

Markov链修正前后的误差计算结果如表 4所示。

由图 6得，经Markov修正后的模型预测曲线和真实值曲线拟合度很高，直观地可以看到主要对2个时间点的运行能耗值进行了修正。从图 7预测结果误差对比中可以得到证实，经过Markov链修正之后得误差非常小，不仅对部分累积误差进行修正接近于0值，其他误差也进行了相应修正，总体预测精度显著提高。由表 4可以看到，经Markov修正之后的RMSE和MAPE均降低，这说明经过Markov链对组合模型中参数传递存在的累积误差进行了有效改进，它的修正过程适用于办公建筑运行能耗这一非线性系统。Markov链的无效性特征在组合模型的误差修正中得到了充分体现，只对最后的预测结果做处理而不影响传统预测模型的运算过程，更符合实际办公建筑运行能耗的变化规律。因此，与Chaos-SVR时间序列预测模型相比，基于Markov链修正的Chaos-SVR时间序列预测模型更优良。

4 结论

针对办公建筑运行能耗预测问题，本文根据混沌理论构造非线性映射对办公建筑运行能耗时间序列进行相空间重构，提出Chaos-SVR的组合预测模型，并利用Markov链对组合模型中累积误差进行修正的方法。并将预测模型应用于西安某实际的办公建筑，得出以下结论：

1) 用C-C法求解办公建筑总能耗时间序列嵌入维数和延迟时间，对其进行了相空间重构。同时根据小数据量法求解最大Lyapunov指数值为0.007 7，由此判定办公建筑能耗时间序列存在混沌特性。

2) 将常用的预测方法比如Narnet, SVR与本文所建立方法Chaos-SVR预测结果进行比较，本文所建立方法预测误差最小，RMSE为6.326 1，MAPE为0.287 5，利用Markov链对组合模型累积误差修正之后的RMSE和MAPE分别降为0.565 5和0.044 1。结果表明，Markov修正之后的预测精度显著提高，为办公建筑的优化运营与节能管理提供充分的决策依据，符合工程实际应用。