人们对建筑环境舒适性要求提高的同时建筑能耗也随之增加，2022年中国建筑能耗高达22.7×10⁸ t标准煤当量，占全国总能耗的45.5%^[1]。光伏建筑一体化(building integrated photovoltaic, BIPV)通过将光伏阵列安装在建筑围护结构外表面，为建筑提供清洁、绿色、环保的电能，成为当今建筑节能减排的重要方式^[2]。能源系统是以电力系统为核心，耦合光伏、热、冷等多种能源形成的“源-网-荷-储”一体化系统，可以有效提升能源利用率、促进光伏发电消纳^[3]。随着光伏发电装机比例的不断上升^[4-5]，光伏出力和负荷需求的不确定性为系统稳定运行带来巨大挑战^[6]，如何有效应对不确定性因素，实现系统安全经济运行成为能源系统优化调度热点问题。

鲁棒优化(robust optimization，RO)是优化理论中, 在不确定条件下使优化问题求解结果具有一定程度鲁棒性的方法^[7-8]。RO以不确定集合的形式描述系统中不确定变量的波动范围，充分考虑了系统建模过程中的不确定因素。相比其他优化方法，RO能够通过优化手段求得“最恶劣”场景下系统的调度方案，牺牲少量运行成本提高整个系统鲁棒性，因此受到国内外学者的广泛关注。Qiu等^[9]提出一种三级两阶段混合整数鲁棒优化模型，通过嵌套R&D算法求解得到正常和最恶劣突发场景下的最优调度方案。Chen等^[10]构建了一种电力市场考虑需求响应下的冷热电联产能源优化管理框架，通过设置不确定参数有效调节不确定系统的鲁棒性和经济性。刘一欣等^[11]建立了考虑微电网内可再生分布式电源和负荷不确定性的两阶段鲁棒优化经济调度模型，得到系统“最恶劣”场景下最小运行成本的调度方案。周特等^[12]针对综合能源系统低碳运行中碳-源-荷多个不确定参数带来的问题，构建基于两阶段鲁棒优化的综合能源系统低碳经济调度模型，通过算例验证了所提鲁棒优化方法的低碳经济性。不确定集是RO的核心，文献[9-12]均以不确定变量预测结果及其最大波动范围建立不确定集，忽视了预测精度等预测特征对调度结果可靠性的影响，割裂了预测和优化调度间的联系，针对此问题，Wang等^[13]在多能源参与的微电网电力市场中建立了预测-日前-实时多阶段最优能源管理框架, 通过长短期记忆神经网络对微电网中风电、光伏和电热负荷进行预测，降低不确定集保守性，将日前调度成本降至最低。Qin等^[14]提出一种结合能量回收和预测的微电网两阶段鲁棒优化调度框架，根据BP神经网络预测结果建立光伏、风电和负荷不确定集，最后通过仿真实验说明采用BP神经网络的预测误差在5%以内。文献[13-14]虽将预测阶段纳入优化调度体系中，但仅是对各不确定变量单独预测，忽略了可再生能源和不同负荷间的耦合关系，且未将预测结果中包含的概率信息融入调度过程中，影响预测精度的同时降低了系统运行经济型和鲁棒性。

本文所述BIPV能源系统内包含多种能源负荷，且各能源负荷互相耦合，若对各能源负荷逐一单独预测，割裂多元源荷之间的联系，忽视不同预测间的隐藏信息，对预测结果的准确度影响极大。多任务学习(multi-task learning，MTL)能够在具有一定相关性的多个不同任务之间共享信息，提高预测模型精度，多任务高斯过程^[15-16](muti-task Gaussian process，MTGP)是MTL领域中重要的研究方法，该方法考虑不同任务间相关性的同时深入挖掘历史数据的概率信息，得到预测结果的均值和方差，生成各预测任务在不同置信度下的置信区间。

综上所述，本文针对BIPV能源系统优化调度中光伏出力和电、热、冷负荷不确定性问题提出一种预测—配置—调度一体化策略。首先，充分考虑光伏出力和电、热、冷负荷间的耦合关系，根据其历史数据建立基于MTGP的源荷预测模型，通过各预测任务间相关性矩阵量化描述源荷耦合特性，得到精确的源荷预测均值和方差，建立正态概率分布形式的新型源荷不确定集。其次，针对中小型BIPV设计符合其建筑特性的能源系统，以系统年投资成本和年运营成本最小化为目标函数建立系统容量优化配置模型，求解得到系统经济最优场景下各设备容量配置^[17-18]，最后，在确定源荷不确定集和能源系统最优容量配置的基础上，建立系统两阶段鲁棒优化调度模型并采用列与约束生成算法和KKT(Karush-Kuhn-Tucker)条件迭代求解得到最恶劣场景下的经济最优调度方案。通过仿真实验对源荷预测结果、系统容量优化配置结果、鲁棒优化调度结果和不同优化方案对比结果进行分析。

1 BIPV能源系统建模

针对中小型BIPV光伏出力和电热冷负荷需求特性，本文设计图 1所示BIPV能源系统，该系统由光伏(photovoltaic，PV)、配电网(grid)、空气源热泵(air source heat pump，ASHP)、电制冷机(electric chiller，EC)、储能(energy storage，ES)、电、热、冷负荷多种能源和负荷组成。电能由配电网购买和光伏发电获取，多余电能可储存在储电(electricity energy storage，EES)设备中，除满足电负荷需求外，电能同时供应ASHP和EC用电需求。热能和冷能分别由ASHP和EC提供，在配电网电价较低时可以将多余的热、冷能以热水和冷冻水的形式储存于储热(heating energy storage，HES)、储冷(cooling energy storage，CES)设备中，满足热、冷负荷需求。

由图 1可知，电能可通过制热、制冷设备转换为热冷能，所以储电对储热和储冷是有替代作用的，但蓄电设备自身存在容量、单位时间内最大充放电功率等限制，且其投资成本高于储热、储冷设备，过度依赖储电设备不利于系统的经济运行。因此后文中综合考虑了储能设备间的耦合关系，建立系统容量优化配置模型求解得到系统设备最优容量配置。

1.1 光伏(PV)

PV利用半导体的光电效应，将太阳辐射能转化为电能，太阳能取之不尽，用之不竭，是最常见的可再生能源。PV满足以下约束：

$ P_{\mathrm{PV}}(t)=\eta_{\mathrm{PV}} n_{\mathrm{PV}} r(t) $

(1)

$ n_{\mathrm{PV}} \leqslant n_{\text {arca }} $

(2)

式中：η_PV为光伏逆变器效率，n_PV为光伏板面积，r(t)为t时刻太阳辐射量，n_area为建筑可施工面积。

1.2 配电网(grid)

BIPV能源系统内光伏发电和EES放电无法满足电负荷需求时，需从配电网购电。光伏发电量富余时可向配电网出售多余电能，降低系统运行成本。系统和配电网的交互功率满足以下约束：

$ 0 \leqslant P_{\mathrm{G}}^{\text {buy }}(t) \leqslant U_{\mathrm{G}}(t) P_{\mathrm{G}}^{\max } $

(3)

$ 0 \leqslant P_{\mathrm{G}}^{\text {sell }}(t) \leqslant\left[1-U_{\mathrm{G}}(t)\right] P_{\mathrm{G}}^{\max } $

(4)

式中：P_G^buy(t)、P_G^sell(t)分别为t时刻系统向配电网购售电功率，U_G(t)为t时刻系统与配电网的购售电状态，取值为1时表示系统从配电网购电，取值为0时表示系统向配电网售电；P_G^max为系统与配电网交互功率的最大值。

1.3 空气源热泵(ASHP)

ASHP是一种利用电能使低位热源空气流向高位热源的节能装置，可以将空气中不能直接利用的低位热能转化为直接利用的高位热能。ASHP满足以下约束：

$ 0 \leqslant P_{\mathrm{ASHP}}(t) \leqslant \bar{P}_{\mathrm{ASHP}} $

(5)

$ Q_{\text {ASHP }}(t)=\operatorname{COP}_{\text {ASHP }} P_{\text {ASHP }}(t) $

(6)

式中：P_ASHP(t)为t时刻ASHP输出热功率，P_ASHP为ASHP额定功率，Q_ASHP(t)为t时刻ASHP输出热量，COP_ASHP为ASHP能效系数。

1.4 电制冷机(EC)

EC通过将电能转化为冷能，稳定调节室内温度，避免出现室内温度波动，并且无需特定环境条件和设施，因此适用各种场合。EC满足以下约束：

$ 0 \leqslant P_{\mathrm{EC}}(t) \leqslant \bar{P}_{\mathrm{EC}} $

(7)

$ Q_{\mathrm{EC}}(t)=\operatorname{COP}_{\mathrm{EC}} P_{\mathrm{EC}}(t) $

(8)

式中：P_EC(t)为t时刻EC输出冷功率，P_EC为EC额定功率，Q_EC(t)为t时刻EC输出冷量，COP_EC为EC能效系数。

1.5 储能设备(ES)

ES可以实现能量跨时间转移，当系统处于电价低或可再生能源出力较高时段时，ES进行能量存储，当电价较高或可再生能源出力较低时，ES释放能量满足负荷需求，实现系统的低成本、稳定运行。ES满足以下约束：

$ 0 \leqslant P_{\mathrm{ES}}^{\mathrm{dis}}(t) \leqslant U_{\mathrm{ES}}(t) P_{\mathrm{ES}}^{\max } $

(9)

$ 0 \leqslant P_{\mathrm{ES}}^{\mathrm{ch}}(t) \leqslant\left[1-U_{\mathrm{ES}}(t)\right] P_{\mathrm{ES}}^{\max } $

(10)

$ \eta_{\mathrm{ES}} \sum\limits_{t=1}^{N_{\mathrm{T}}}\left[P_{\mathrm{ES}}^{\mathrm{ch}}(t) \Delta t\right]-\frac{1}{\eta_{\mathrm{ES}}} \sum\limits_{t=1}^{N_{\mathrm{T}}}\left[P_{\mathrm{ES}}^{\mathrm{dis}}(t) \Delta t\right]=0 $

(11)

$ \begin{gathered} E_{\mathrm{ES}}^{\min } \leqslant E_{\mathrm{ES}}(0)+\eta_{\mathrm{ES}} \sum\limits_{t^{\prime}=1}^t\left[P_{\mathrm{ES}}^{\mathrm{ch}}(t) \Delta t\right]- \\ \frac{1}{\eta_{\mathrm{ES}}} \sum\limits_{t^{\prime}=1}^t\left[P_{\mathrm{ES}}^{\mathrm{dis}}(t) \Delta t\right] \leqslant E_{\mathrm{ES}}^{\max } \end{gathered} $

(12)

式中：P_ES^ch(t)、P_ES^dis(t)分别为t时刻ES充放功率，P_ES^max为ES单位时间内最大充放功率，U_ES(t)为t时刻ES充放状态，取值为1时表示ES放能，取值为0时表示ES充能，式(11)表示一个调度周期内初始和结束时刻ES剩余容量相等，保证调度的周期性和持续性；η_ES为ES充放效率，N_T为调度周期，取24 h；Δt为调度步长，取1 h；E_ES(0)为调度初始容量，E_ES^min、E_ES^max分别为ES调度过程中最大、最小剩余容量。式中ES包括EES、HES、CES设备。

1.6 系统功率平衡

系统功率平衡为：

$ \begin{gathered} P_{\mathrm{PV}}(t)+P_{\mathrm{G}}^{\text {buy }}(t)+P_{\mathrm{EES}}^{\mathrm{dis}}(t)=P_{\mathrm{G}}^{\mathrm{sell}}(t)+P_{\mathrm{EES}}^{\mathrm{ch}}(t)+ \\ P_{\mathrm{ASHP}}(t)+P_{\mathrm{EC}}(t)+P_{\mathrm{EL}}(t) \end{gathered} $

(13)

$ Q_{\mathrm{ASHP}}(t)+P_{\mathrm{HES}}^{\mathrm{dis}}(t)=P_{\mathrm{HES}}^{\mathrm{ch}}(t)+P_{\mathrm{HL}}(t) $

(14)

$ Q_{\mathrm{EC}}(t)+P_{\mathrm{CES}}^{\mathrm{dis}}(t)=P_{\mathrm{CES}}^{\mathrm{ch}}(t)+P_{\mathrm{CL}}(t) $

(15)

式中：式(13)~(15)分别表示t时刻电、热、冷平衡约束，P_EL(t)、P_HL(t)、P_CL(t)分别为t时刻电、热、冷负荷。

2 BIPV能源系统鲁棒优化调度模型

本文提出的预测—配置—调度一体化BIPV能源系统鲁棒优化调度模型中，首先建立了基于MTGP的源荷不确定集，然后对能源系统各设备进行容量最优配置，建立以年度总成本最小化为目标函数的容量优化配置模型，最后提出日前两阶段鲁棒优化调度模型。

2.1 源荷不确定集模型

鲁棒优化调度中，不确定集直接影响系统的优化调度结果，本文通过MTGP建立源荷不确定集描述BIPV能源系统中光伏出力和电、热、冷负荷波动引起的不确定性。

2.1.1 基于MTGP的日前预测模型

高斯过程(Gaussian process，GP)可以看作一个随机过程。假设训练集输入向量P ={p_i|i=1, 2, …, N}，训练集输出向量Q ={q_i|i=1, 2, …, N}，高斯过程f(p)由一个均值函数m(p)和一个协方差函数k(p, p′)组成，可表示为

$ f(p) \sim \operatorname{GP}\left(m(p), k\left(p, p^{\prime}\right)\right) $

(16)

式中：m(p)为f(p)在点p处的期望，其结果并不影响预测结果，为简便计算可假设m(p)=0；k(p, p′)为核函数。

输入p_i和输出q_i间的映射关系f(p)不能直接获得，可表示为

$ q=f(p)+\varepsilon $

(17)

式中ε为观测噪声，服从均值为0，方差为σ_n²的高斯分布，记作ε~N(0, σ_n²)。

由于f(p)服从多元高斯分布，q可表示为

$ q \sim N\left(0, k\left(p, p^{\prime}\right)+\sigma_n^2 \delta_{i j}\right) $

(18)

式中δ_ij为克罗内克函数，当i=j时取1，其余情况取0。

假设测试集输入P_*与训练集输入P独立同分布，那么测试集输出Q_*与训练集输出Q的联合高斯分布可表示为

$ \left[\begin{array}{l} \boldsymbol{Q} \\ \boldsymbol{Q}_* \end{array}\right] \sim N\left(0, \left[\begin{array}{cc} \boldsymbol{K}(\boldsymbol{P}, \boldsymbol{P})+\sigma_n^2 \boldsymbol{I} & \boldsymbol{K}\left(\boldsymbol{P}, \boldsymbol{P}_*\right) \\ \boldsymbol{K}\left(\boldsymbol{P}_*, \boldsymbol{P}\right) & \boldsymbol{K}\left(\boldsymbol{P}_*, \boldsymbol{P}_*\right) \end{array}\right]\right) $

(19)

式中：I为单位矩阵，K(P, P)、K(P, P_*)、K(P_*, P_*) 分别为训练集与训练集、训练集与测试集、测试集与测试集输入之间的协方差矩阵。

根据多元高斯分布的性质可知：

$ \boldsymbol Q_* \mid \boldsymbol{Q} \sim N\left(\boldsymbol{\mu}_*, \boldsymbol{\sigma}_*\right) $

(20)

其中

$ \left\{\begin{array}{l} \boldsymbol{\mu}_*=\boldsymbol{K}\left(\boldsymbol{P}, \boldsymbol{P}_*\right)^{\mathrm{T}}\left(\boldsymbol{K}(\boldsymbol{P}, \boldsymbol{P})+\boldsymbol\sigma_n^2 \boldsymbol I\right)^{-1} \boldsymbol Q \\ \boldsymbol{\sigma}_*=\boldsymbol{K}\left(\boldsymbol{P}_*, \boldsymbol{P}_*\right)-\boldsymbol{K}\left(\boldsymbol{P}_{, } \boldsymbol{P}_*\right)^{\mathrm{T}}\left(\boldsymbol{K}(\boldsymbol{P}, \boldsymbol{P})+\boldsymbol\sigma_n^2 \boldsymbol I\right)^{-1} \boldsymbol{K}\left(\boldsymbol{P}, \boldsymbol{P}_*\right) \end{array}\right. $

式中，μ_*、σ_*分别为预测值的期望和方差。

MTGP是由GP扩展得到的一种多任务学习模型，假设训练集输入向量P ={p_i^j|j=1, 2, …, M; i=1, 2, …, N^j}，训练集输出向量Q ={q_i^j|j=1, 2, …, M; i=1, 2, …, N^j}，其中：M为任务个数，N^j为任务j输入向量的长度。为表明p_i和q_i与任务j的隶属关系，引入附加标签L到输入中，其中L ={l_i^j|l_i^j=j, j=1, 2, …, M; i=1, 2, …, N^j}，训练集D ={(l_i^j, p_i^j, q_i^j)|j=1, 2, …, M; i=1, 2, …, N^j}，新的协方差函数可表示为

$ k_{\mathrm{MTGP}}\left(l, l^{\prime}, p, p^{\prime}\right)=k_l\left(l, l^{\prime}\right) \times k_x\left(p, p^{\prime}\right) $

(21)

任务j中输入(l_i^j, p_i^j)和输出q_i^j间的映射关系可表示为

$ q^j=f^j\left(l^j, p^j\right)+\varepsilon^j $

(22)

式中ε^j为任务j的观测噪声，服从均值为0，方差为σ_n^j2的高斯分布，记作ε^j~N(0, σ_n^j2)。

为简化计算，假设N^j=N(j=1, 2, …, M)，所有训练集任务间的协方差矩阵可表示为

$ \boldsymbol{K}_{\mathrm{MTGP}}\left(\boldsymbol{L}, \boldsymbol{P}, \boldsymbol{\theta}_l, \boldsymbol{\theta}_x\right)=\boldsymbol{K}_l\left(\boldsymbol{L}, \boldsymbol{\theta}_l\right) \otimes \boldsymbol{K}_p\left(\boldsymbol{P}, \boldsymbol{\theta}_p\right) $

(23)

式中：$\otimes$为克罗内克积，θ_l、θ_p分别为任务相关性矩阵K_l和GP协方差矩阵K_p中的超参数。根据下式可以看出，如果K_l为对角矩阵，各任务间没有相关性，式(23)退化为M个相互独立的GP核函数。K_l的常用形式可表示为

$ \boldsymbol{K}_l=\left(\begin{array}{cccc} \theta_1 & 0 & \cdots & 0 \\ \theta_2 & \theta_3 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ \theta_{k-m+1} & \theta_{k-m+2} & \cdots & \theta_k \end{array}\right) \cdot\left(\begin{array}{cccc} \theta_1 & 0 & \cdots & 0 \\ \theta_2 & \theta_3 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ \theta_{k-m+1} & \theta_{k-m+2} & \cdots & \theta_k \end{array}\right)^{\mathrm{T}} $

(24)

参考GP中对输出期望和方差的推断(20)，结合式(23)，MTGP对于给定测试集(x_*, l_*)，其输出的期望和方差可表示为

$ \left\{\begin{array}{l} \boldsymbol{\mu}_*=\left(\boldsymbol{k}_{l_*} \otimes \boldsymbol{k}_{p_*}\right)^{\mathrm{T}}\left(\boldsymbol{K}_{\mathrm{MTGP}}+\boldsymbol{R} \otimes \boldsymbol{I}\right)^{-1} \boldsymbol{Q} \\ \boldsymbol{\sigma}_*=\boldsymbol{k}_{l_{**}} \otimes \boldsymbol{k}_{p_{**}}-\left(\boldsymbol{k}_{l_*} \otimes \boldsymbol{k}_{p_*}\right)^{\mathrm{T}}\left(\boldsymbol{K}_{\mathrm{MTGP}}+\boldsymbol{R} \otimes \boldsymbol{I}\right)^{-1}\left(\boldsymbol{k}_{l_*} \otimes \boldsymbol{k}_{p_*}\right) \end{array}\right. $

(25)

式中：k_l*为K_l的第l_*列，k_p*、k_p**分别为训练集与测试集、测试集与测试集输入之间的协方差向量，R为对角线元素为σ_n^j2的对角矩阵。

2.1.2 基于MTGP的源荷不确定集

传统盒式不确定集^{[11, 19]}通过预测结果的等比例缩放描述源荷端波动的不确定性，然而该不确定集并未考虑预测结果中均值和方差等概率特性，导致优化调度模型忽视了日前预测和日前调度之间的对应关系。合理利用源荷预测得到的概率特性能有效降低优化调度决策结果的保守性，为源荷预测与优化调度建立紧密联系。由MTGP的日前预测模型可知，利用MTGP在源荷预测阶段可以充分考虑光伏出力和电热冷负荷间的耦合关系，建立各预测任务间相关性矩阵K_l量化描述光伏出力和电、热、冷负荷间的相关性，得到更精确的预测均值和方差。光伏出力和电、热、冷负荷预测结果均可由正态分布概率模型描述，因此本文提出以下基于MTGP的新型源荷不确定集：

$ \boldsymbol{U}=\left\{\begin{array}{l} u_{\mathrm{PV}}(t) \in\left[\hat{u}_{\mathrm{PV}}(t)-z_{\alpha / 2} \sigma_{\mathrm{PV}}(t), \hat{u}_{\mathrm{PV}}(t)+z_{\alpha / 2} \sigma_{\mathrm{PV}}(t)\right] \\ u_{\mathrm{EL}}(t) \in\left[\hat{u}_{\mathrm{EL}}(t)-z_{\alpha / 2} \sigma_{\mathrm{EL}}(t), \hat{u}_{\mathrm{EL}}(t)+z_{\alpha / 2} \sigma_{\mathrm{EL}}(t)\right] \\ u_{\mathrm{HL}}(t) \in\left[\hat{u}_{\mathrm{HL}}(t)-z_{\alpha / 2} \sigma_{\mathrm{HL}}(t), \hat{u}_{\mathrm{HL}}(t)+z_{\alpha / 2} \sigma_{\mathrm{HL}}(t)\right] \\ u_{\mathrm{CL}}(t) \in\left[\hat{u}_{\mathrm{CL}}(t)-z_{\alpha / 2} \sigma_{\mathrm{CL}}(t), \hat{u}_{\mathrm{CL}}(t)+z_{\alpha / 2} \sigma_{\mathrm{CL}}(t)\right] \end{array}\right. $

(26)

式中：u =[u_PV(t), u_EL(t), u_HL(t), u_CL(t)]∈ R ^NT×4，其中u_PV、u_EL、u_HL、u_CL为考虑不确定性后光伏出力和电、热、冷负荷的不确定变量；$\hat{u}_{\mathrm{PV}}、\hat{u}_{\mathrm{EL}}、\hat{u}_{\mathrm{HL}}、\hat{u}_{\mathrm{CL}}$为预测期望，σ_PV²、σ_EL²、σ_HL²、σ_CL²为预测方差，z_α/2为标准正态分布的双侧α分位数，其中α为显著性水平，可以通过α调节不确定集的上、下边界。

2.2 BIPV能源系统容量优化配置模型

本文建立式(27)所示系统容量优化配置模型，其目的在于确定系统年度总成本最优情况下各设备容量配置值，年度总成本包括年投资成本(capital expenditure，CAPEX)和年运营成本(operating expenses，OPEX)。

$ \min C=\min \left(C_{\text {CAPEX }}+C_{\mathrm{OPEX}}\right) $

(27)

其中

$ \left\{\begin{array}{l} C_{\mathrm{CAPEX}}=C_{\mathrm{CAPEX}}^{\mathrm{PV}}+C_{\mathrm{CAPEX}}^{\mathrm{ASHP}}+C_{\mathrm{CAPEX}}^{\mathrm{EC}}+C_{\mathrm{CAPEX}}^{\mathrm{ES}} \\ C_{\mathrm{OPEX}}=C_{\mathrm{OPEX}}^{\mathrm{G}}+C_{\mathrm{OPEX}}^{\mathrm{ES}} \end{array}\right. $

C_CAPEX^PV、C_CAPEX^ASHP、C_CAPEX^EC、C_CAPEX^ES分别为PV、ASHP、EC、ES设备年投资成本，可表示为：

$ C_{\mathrm{CAPEX}}^{\mathrm{PV}}=c_{\mathrm{PV}} n_{\mathrm{PV}} \mathrm{CRF}_{\mathrm{pv}} $

(28)

$ C_{\mathrm{CAPEX}}^{\mathrm{ASHP}}=c_{\mathrm{ASHP}} \bar{P}_{\mathrm{ASHP}} \mathrm{CRF}_{\mathrm{ASHP}} $

(29)

$ C_{\text {CAPEX }}^{\mathrm{EC}}=c_{\mathrm{EC}} \bar{P}_{\mathrm{EC}} \mathrm{CRFF}_{\mathrm{EC}} $

(30)

$ \begin{aligned} C_{\mathrm{CAPEX}}^{\mathrm{ES}}= & c_{\mathrm{EES}} E_{\mathrm{EES}}^{\max } \mathrm{CRF}_{\mathrm{EES}}+c_{\mathrm{HES}} E_{\mathrm{HES}}^{\max } \mathrm{CRF}_{\mathrm{HES}}+ c_{\mathrm{IES}} E_{\mathrm{CES}}^{\max } \mathrm{CRF}_{\mathrm{IES}} \end{aligned} $

(31)

式中：n_pv为光伏板面积，P_ASHP为ASHP额定功率，P_EC为EC额定功率，E_EES^max、E_HES^max、E_CES^max分别为EES、HES、CES容量，c_j为设备j的投资成本，CRF为资金回复系数，可表示为

$ \text { CRF }=\frac{i(1+i)^r}{(1+i)^r-1} $

(32)

式中：i为折现率，r为设备使用寿命。

C_OPEX^G、C_OPEX^ES分别为配电网交互、ES设备年运营成本，可表示为

$ C_{\mathrm{OPEx}}^{\mathrm{G}}=\sum\limits_{t=1}^{8760}\left[k_{\mathrm{G}}^{\text {buy }}(t) P_{\mathrm{G}}^{\text {buy }}(t)-k_{\mathrm{G}}^{\text {sell }}(t) P_{\mathrm{G}}^{\text {sell }}(t)\right] $

(33)

$ \begin{aligned} C_{\mathrm{OPEX}}^{\mathrm{ES}}= & k_{\mathrm{EES}} \sum\limits_{t=1}^{8760}\left[P_{\mathrm{EES}}^{\mathrm{dis}}(t) / \eta_{\mathrm{EES}}+P_{\mathrm{EES}}^{\mathrm{ch}}(t) \eta_{\mathrm{EES}}\right]+ \\ & k_{\mathrm{HES}} \sum\limits_{t=1}^{8760}\left[P_{\mathrm{HES}}^{\mathrm{dis}}(t) / \eta_{\mathrm{HES}}+P_{\mathrm{HES}}^{\mathrm{dh}}(t) \eta_{\mathrm{HES}}\right]+ \\ & k_{\mathrm{CES}} \sum\limits_{t=1}^{8760}\left[P_{\mathrm{CES}}^{\mathrm{dis}}(t) / \eta_{\mathrm{CES}}+P_{\mathrm{CES}}^{\mathrm{ch}}(t) \eta_{\mathrm{CES}}\right] \end{aligned} $

(34)

式中：k_G^buy(t)、k_G^sell(t)分别为t时刻配电网交易电价，k_EES、k_HES、k_CES分别为折算后的单位充放电、热、冷能成本。

式(27)所示系统容量优化配置模型决策变量W可表示为

$ \begin{aligned} W= & {\left[n_{\mathrm{pv}}, \bar{P}_{\mathrm{ASHP}}, \bar{P}_{\mathrm{EC}}, E_{\mathrm{EES}}^{\max }, E_{\mathrm{HES}}^{\max }, E_{\mathrm{CES}}^{\max }, U_{\mathrm{G}}(t), \right.} \\ & U_{\mathrm{EES}}(t), U_{\mathrm{HES}}(t), U_{\mathrm{CES}}(t), P_{\mathrm{G}}^{\mathrm{buy}}(t), P_{\mathrm{G}}^{\mathrm{sell}}(t), \\ & P_{\mathrm{EES}}^{\mathrm{ch}}(t), P_{\mathrm{EES}}^{\mathrm{dis}}(t), P_{\mathrm{HES}}^{\mathrm{ch}}(t), P_{\mathrm{HES}}^{\mathrm{dis}}(t), P_{\mathrm{ASHP}}(t), \\ & \left.P_{\mathrm{CES}}^{\mathrm{ch}}(t), P_{\mathrm{CES}}^{\mathrm{dis}}(t), P_{\mathrm{EC}}(t)\right]^{\mathrm{T}} \end{aligned} $

(35)

式中，t=1, 2, …, 8 760。

2.3 BIPV能源系统两阶段鲁棒优化调度模型

本文建立式(36)所示两阶段鲁棒优化调度模型，其目的在于找到不确定变量u在不确定集U中最恶劣场景下的经济性最优调度方案，并满足式(1)~(15)所述约束条件:

$ \min _{\boldsymbol{x}}\left\{\left(\max _{\boldsymbol{u} \in \boldsymbol{U}} \min _{\boldsymbol{y} \in \boldsymbol{\Omega ( x , \boldsymbol { u } )}} C^{\mathrm{G}}+C^{\mathrm{ES}}\right)\right\} $

(36)

式中，C^G、C^ES分别为配电网交互成本和储能成本。

C^G、C^ES可表示为：

$ C^{\mathrm{G}}=\sum\limits_{t=1}^{N_{\mathrm{T}}}\left[k_{\mathrm{G}}^{\text {buy }}(t) P_{\mathrm{G}}^{\text {buy }}(t)-k_{\mathrm{G}}^{\text {sell }}(t) P_{\mathrm{G}}^{\text {sell }}(t)\right] $

(37)

$ \begin{aligned} C^{\mathrm{ES}}= & k_{\mathrm{EES}} \sum\limits_{t=1}^{N_{\mathrm{T}}}\left[P_{\mathrm{EES}}^{\mathrm{dis}}(t) / \eta_{\mathrm{EES}}+P_{\mathrm{EES}}^{\mathrm{ch}}(t) \eta_{\mathrm{EES}}\right]+ \\ & k_{\mathrm{HES}} \sum\limits_{t=1}^{N_{\mathrm{T}}}\left[P_{\mathrm{HES}}^{\mathrm{dis}}(t) / \eta_{\mathrm{HES}}+P_{\mathrm{HES}}^{\mathrm{ch}}(t) \eta_{\mathrm{HES}}\right]+ \\ & k_{\mathrm{CES}} \sum\limits_{t=1}^{N_{\mathrm{T}}}\left[P_{\mathrm{CES}}^{\mathrm{dis}}(t) / \eta_{\mathrm{CES}}+P_{\mathrm{CES}}^{\mathrm{ch}}(t) \eta_{\mathrm{CES}}\right] \end{aligned} $

(38)

紧凑形式可表示为

$ \left\{\begin{aligned} \min _x & \left\{\max _{\boldsymbol{u} \in \boldsymbol{U}} \min _{\boldsymbol{y \in \Omega ( \boldsymbol x ,\boldsymbol u )}} \boldsymbol{c}^{\mathrm{T}} \boldsymbol{y}\right\} \\ \text { s. t. } & \boldsymbol{A y} \geqslant \boldsymbol{a} \\ & \boldsymbol{B y}=\boldsymbol{b} \\ & \boldsymbol{C x}+\boldsymbol{E} \boldsymbol{y} \geqslant \boldsymbol{e} \\ & \boldsymbol{F}_u \boldsymbol{y}=\hat{\boldsymbol{u}} \end{aligned}\right. $

(39)

其中

$ \left\{\begin{aligned} \boldsymbol{x}= & {\left[U_{\mathrm{G}}(t), U_{\mathrm{EES}}(t), U_{\mathrm{HES}}(t), U_{\mathrm{CES}}(t)\right]^{\mathrm{T}} } \\ \boldsymbol{y}= & {\left[P_{\mathrm{G}}^{\mathrm{buy}}(t), P_{\mathrm{G}}^{\mathrm{elll}}(t), P_{\mathrm{EES}}^{\mathrm{ch}}(t), P_{\mathrm{EES}}^{\mathrm{dis}}(t), P_{\mathrm{HES}}^{\mathrm{ch}}(t), P_{\mathrm{HES}}^{\mathrm{dis}}(t), \right.} \\ & P_{\mathrm{ASHP}}(t), Q_{\mathrm{ASHP}}(t), P_{\mathrm{CES}}^{\mathrm{ch}}(t), P_{\mathrm{CES}}^{\mathrm{dis}}(t), P_{\mathrm{EC}}(t), Q_{\mathrm{EC}}(t), \\ & \left.P_{\mathrm{PV}}(t), P_{\mathrm{EL}}(t), P_{\mathrm{HL}}(t), P_{\mathrm{CL}}(t)\right]^{\mathrm{T}} \end{aligned}\right. $

式中：U为源荷不确定集，Ω (x, u) 为给定(x, u) 时y的可行域，c为目标函数对应的系数向量，A、B、C、E、F_u为对应约束的系数矩阵，a、b、e为常数向量，$\hat{\boldsymbol{u}}$为对应场景下不确定变量的取值，x、y为决策变量, t=1, 2, …, N_T。

两阶段鲁棒优化调度模型(39)中，第1阶段决策变量x是对配电网购售电和储能充放能做出决策，第2阶段优化变量u和y是求解最恶劣场景下的最优解，式中max结构是为了找到运营成本最大时的最恶劣场景。

3 两阶段鲁棒优化调度模型求解

本文建立的BIPV能源系统两阶段鲁棒优化调度模型为min-max-min结构的两阶段3层优化模型，其结构复杂，无法同时求解，因此采用列约束生成(column and constraint generation，C&CG)算法将式(39)分解为主问题(master problem，MP)和子问题(subproblem，SP)反复迭代求解直至运营成本上、下界收敛。

分解后得到主问题可表示为

$ \left\{\begin{aligned} \min _x \theta & \\ \text { s.t. } & \theta \geqslant \boldsymbol{c}^{\mathrm{T}} \boldsymbol{y}_i \\ & \boldsymbol{A} \boldsymbol{y}_i \geqslant \boldsymbol{a} \\ & \boldsymbol{B} \boldsymbol{y}_i=\boldsymbol{b} \\ & \boldsymbol{C x}+\boldsymbol{E} \boldsymbol{y}_i \geqslant \boldsymbol{e} \\ & \boldsymbol{F}_u \boldsymbol{y}_i=\boldsymbol{u}_i^* \\ & i \leqslant z \end{aligned}\right. $

(40)

式中：θ为辅助变量，z为当前迭代次数，y_i为i-1次迭代中子问题的解，u_i^*为i-1次迭代中得到最恶劣场景u的取值。

分解后得到子问题可表示为

$ \max\limits_{\boldsymbol{u} \in \boldsymbol{U}} \min\limits_{\boldsymbol{y} \in \Omega(\boldsymbol{x}, \boldsymbol{u})} \boldsymbol{c}^{\mathrm{T}} \boldsymbol{y} $

(41)

式中双层max-min问题难以直接求解，使用KKT条件^[20]将内层min问题转化为无目标的方程组，与外层max问题合并，得到以下形式：

$ \begin{cases}\max\limits_\limits{\boldsymbol{u} \in \boldsymbol{U}} & \boldsymbol{c}^{\mathrm{T}} \boldsymbol{y} \\ \text { s.t. } & \boldsymbol{A} \boldsymbol{y} \geqslant \boldsymbol{a} \\ & \boldsymbol{B} \boldsymbol{y}=\boldsymbol{b} \\ & \boldsymbol{C x}+\boldsymbol{E} \boldsymbol{y} \geqslant \boldsymbol{e} \\ & \boldsymbol{F}_{\boldsymbol{u}} \boldsymbol{y}=\boldsymbol{u} \\ & \lambda_1, \lambda_2 \geqslant 0 \\ & \lambda_1(-\boldsymbol{A} \boldsymbol{y}+\boldsymbol{a})=0 \\ & \lambda_2(-\boldsymbol{C} \boldsymbol{x}-\boldsymbol{E} \boldsymbol{y}+\boldsymbol{e})=0 \\ & \boldsymbol{c}-\boldsymbol{A}^{\mathrm{T}} \boldsymbol{\lambda}_1-\boldsymbol{E}^{\mathrm{T}} \boldsymbol{\lambda}_2+\boldsymbol{B}^{\mathrm{T}} \omega_1+\boldsymbol{F}_\boldsymbol{u}^{\mathbf{T}} \boldsymbol\omega_2=0\end{cases} $

(42)

式中：λ₁、λ₂为不等式约束引入的辅助变量，ω₁、ω₂为等式约束引入的辅助变量。对于约束条件中的二次项，可通过big-M^[20]法进行线性化。

子问题求得的最优解中，源荷不确定变量为不确定集的一个极点，其取值为不确定集边界值。本文研究的BIPV能源系统中，光伏出力取到最小边界、负荷取到最大边界时系统运营成本更高，符合“最恶劣”场景定义，因此式(26)可改写为

$ \left\{\begin{array}{l} u_{\mathrm{PV}}(t)=\hat{u}_{\mathrm{PV}}(t)-\tau_{\mathrm{PV}}(t) z_{\alpha / 2} \sigma_{\mathrm{PV}}(t), \sum\limits_{t=1}^{N_{\mathrm{T}}} \tau_{\mathrm{PV}}(t) \leqslant \mathit\Gamma_{\mathrm{PV}} \\ u_{\mathrm{EL}}(t)=\hat{u}_{\mathrm{EL}}(t)+\tau_{\mathrm{EL}}(t) z_{\alpha / 2} \sigma_{\mathrm{EL}}(t), \sum\limits_{t=1}^{N_{\mathrm{T}}} \tau_{\mathrm{EL}}(t) \leqslant \mathit\Gamma_{\mathrm{EL}} \\ u_{\mathrm{HL}}(t)=\hat{u}_{\mathrm{HL}}(t)+\tau_{\mathrm{HL}}(t) z_{\alpha / 2} \sigma_{\mathrm{HL}}(t), \sum\limits_{t=1}^{N_{\mathrm{T}}} \tau_{\mathrm{HL}}(t) \leqslant \mathit\Gamma_{\mathrm{HL}} \\ u_{\mathrm{CL}}(t)=\hat{u}_{\mathrm{CL}}(t)+\tau_{\mathrm{CL}}(t) z_{\alpha / 2} \sigma_{\mathrm{CL}}(t), \sum\limits_{t=1}^{N_{\mathrm{T}}} \tau_{\mathrm{CL}}(t) \leqslant \mathit\Gamma_{\mathrm{CL}} \end{array}\right. $

(43)

式中: τ=[τ_PV, τ_EL, τ_HL, τ_CL]^T为二进制变量，取1时表示对应时刻不确定变量取值为不确定集边界值，取0时表示对应时刻不确定变量取值为预测结果；Γ=[Γ_PV, Γ_EL, Γ_HL, Γ_CL]^T为各不确定变量的调节参数，取值为0~N_T范围内的正整数，表示一个调度周期内各不确定变量取到不确定集边界值的总时间，其取值越大，调度方案越保守。

经过上述推导，本文建立的两阶段鲁棒优化调度模型已转化为混合整数线性规划问题，通过C&CG算法流程见图 2。

4 算例分析

本文以美国亚利桑那州州立大学某办公楼为研究对象，针对其建筑类型和负荷需求特性进行能源系统设计及鲁棒优化调度研究。美国亚利桑那州地处沙漠地带，气候干燥，平均温度较高，故冷负荷需求远高于热负荷需求，为此设计如图 1所示带有冷、热、电分布式储能的光伏一体化系统，在充分利用当地丰富太阳能资源的基础上满足建筑的电、热、冷负荷需求，全年冷、热、电负荷曲线见图 3。系统各设备运行参数见表 1，日前分时电价见表 2，所建立的优化模型采用matlabR2023a结合yalmip插件调用商业求解器Gurobi求解器进行求解。

4.1 MTGP能源荷预测结果分析

选取8—10月光伏出力和电、热、冷负荷的历史数据为训练集样本，11月1日的数据作为测试集样本，根据MTGP的日前预测模型中关于MTGP的推导可得不同置信度下光伏出力和电、热、冷负荷的置信区间见图 4。

由图 4可知，随着预测结果置信度的降低，对应置信区间边界逐渐变窄，与实际情况相符。本文选用区间覆盖率(coverage probability，CP)^[21]描述置信区间的可信度，其值越大，表示对应预测置信区间内包含真实数据的点越多，系统鲁棒性越强。设定光伏出力不确定调节参数Γ_PV=6，电、热、冷负荷不确定调节参数Γ_EL=12、Γ_HL=12、Γ_CL=12，此时不同置信度下CP和运营成本见表 3。

由表 3可以看出，相同置信度下热负荷的CP比光伏出力和电冷负荷小很多，这是因为数据来源地区温度较高，热负荷与其他负荷相比需求较小且随机性高，导致热负荷预测误差相对较大。置信度增加的同时光伏出力和电、热、冷负荷的CP随之变大，但对应系统运营成本也相应增加，综合考虑系统CP和经济性，后续仿真试验源荷不确定集均选用95.0%置信区间进行。在实际运营中调度人员可根据实际需求灵活调节不确定变量置信度，综合提升系统经济性和鲁棒性。

4.2 容量优化配置结果分析

考虑系统整体经济性，以年度总成本最小化为目标时，系统容量最优配置结果见表 4。

由表 4可以看出，虽然单位面积光伏板投资成本较高，但优化结果中n_PV取其约束上限，等于建筑可施工面积，这是由于光伏发电作为可再生能源清洁环保且无需额外消耗燃料，因此需尽可能多的铺设光伏板保证光伏最大化利用，符合如今大力发展可再生能源的趋势。优化结果中ASHP额定功率和储热容量远小于EC额定功率和储冷容量，与研究对象所处地域环境温度较高、热负荷远低于冷负荷的实际情况相符。系统中制热、制冷设备均由电能驱动，因此在负荷高峰期储电量不仅用于满足电负荷需求，还可转化为热冷量供应热、冷负荷需求，所以优化结果中储电容量相对较大。由此可以看出该容量优化配置结果实现系统年度总成本最小化的同时符合研究对象实际情况，保证系统能够长期经济运行。

4.3 鲁棒优化调度结果分析

设定源荷不确定集为上述源荷预测95%置信区间，不确定调节参数Γ_PV=6，Γ_EL=12、Γ_HL=12、Γ_CL=12, 系统两阶段鲁棒优化调度结果见图 5、6。

图 5(a)~5(d)分别为系统光伏出力、电、热、冷负荷的最恶劣场景。由图 5(a)可知，最恶劣场景下光伏出力在11:00―16:00这6个时段取到不确定集下界，这是由于该时段光照强度达到峰值，系统内各设备能耗均由光伏发电供应，若此时发生阴影遮挡光伏组件等突发事件导致光伏发电量减少，会对系统内设备稳定运行造成影响，为维持系统功率平衡，临时向电网购电也将增加补偿成本，因此根据优化结果中光伏出力最恶劣场景进行调度，提前预估高峰期因意外事故造成光伏减产可以有效提高系统运行稳定性和经济性。由图 5(b)~5(d)可知，电负荷在11:00―20:00和22:00―23:00这12个时段取到不确定集上界，热负荷在11:00―13:00和15:00―23:00这12个时段取到不确定集上界，冷负荷在11:00―18:00和20:00―23:00这12个时段取到不确定集上界，由此可以看出，优化结果中各负荷最恶劣场景均取在负荷高峰期时段，该时段系统负荷需求大、随机性强，临时增加负荷系统将高价从电网购电，以不确定集上边界调度可以提前对高峰时段临时增加负荷进行规划，减小系统补偿成本。综上所述，本文使用的两阶段鲁棒优化调度方法可以基于源荷不确定集确定系统运行最恶劣场景，进一步降低源荷不确定性对系统运行带来的不利影响，提高系统鲁棒性。

图 6(a)~6(c)分别为电能、热能、冷能优化调度结果。从图 6(a)中可以看出，由于分时电价机制，在01:00―06:00时段，光伏出力为0且电价处于低谷期，系统通过从配电网购电维持运行，由于电价低且负荷较小，EES在此时段充电；在07:00―15:00时段，由于光照强度逐渐增加且电价处于高峰期，系统主要通过光伏出力维持运行，在此期间由于电价过高，系统尽量避免从电网购电，光伏发电不足时优先调用EES供应负荷需求；在16:00―18:00时段，光照强度逐渐减少且电价处于平段期，仅依靠光伏发电远远不能满足负荷需求，此时系统通过从配电网购电维持运行；在19:00―21:00时段，光伏出力为0但电价处于高峰期，此时负荷需求相对较小，系统调用EES维持运行；在22:00―23:00时段，光伏出力为0且电价处于平段期，此时EES剩余容量不满足不满足负荷需求，系统通过从配电网购电维持运行；在24:00时段，光伏出力为0且电价处于低谷期，系统通过从配电网购电维持运行，为了保证调度周期内初始和结束时刻EES剩余容量相等，EES在此时充电，以保证其循环利用。由图 6(b)可以看出，系统在电价高峰期的05:00―07:00时段将多余的热能储存在HES中，等到电价较高的19:00和22:00时段将HES中的热能供应给热负荷需求，降低系统运行成本。图 6(c)中冷能同理。综上所述，系统在电价低谷期向配电网购电维持运行并向ES充能，在电价平段期通过光伏发电和配电网购电协调维持运行，在电价高峰期通过光伏发电和ES放能持运行，利用ES避免了高价购电，实现了系统的经济运行。

4.4 不同优化方案对比分析

为验证本文所提鲁棒优化调度模型的可行性，设置3组调度方案分别进行仿真试验，将常见优化调度模型和本文提出的优化调度模型进行比较。3种方案优化调度的运行成本见表 5。

1) 方案1，不考虑源荷不确定性的确定性优化。

2) 方案2，考虑源荷不确定性的两阶段鲁棒优化调度模型，建立传统盒式不确定集，光伏出力不确定集的上、下边界为预测值的10%，电、热、冷负荷不确定集的上、下边界为预测值的20%。

3) 方案3，考虑源荷不确定性的两阶段鲁棒优化调度模型，建立基于MTGP的源荷不确定集，即本文所提策略。

由表 5可以看出，方案1所提的确定性优化运行成本远少于方案2、3，但这并不说明其调度方案更加优越，方案1以光伏出力和电、热、冷负荷的预测值进行日前优化调度，未考虑源荷预测误差对系统运行带来的影响，到调度日阶段由源荷预测误差带来的能源缺口均需系统从实时市补偿，而实时市场的能源价格一般远高于日前市场。假设实时市场中配电网售电价格为日前市场对应时段的1.5倍^[11]，为满足方案1实时负荷需求，系统需高价购买电量补偿能源缺口，其总成本(￥535.29)远超方案2、3所用鲁棒优化调度。方案2、3均考虑了源荷不确定性对系统运行带来的影响，但方案3建立源荷不确定集过程中充分考虑了预测结果中蕴含的概分布率信息，方案2仅仅利用预测值的等比例缩放，导致其不确定集边界过于宽泛，保守度过高，最终方案3经济性明显优于方案2。

5 结论

本文提出一种计及源荷不确定性的BIPV能源系统预测—配置—调度一体化的鲁棒优化调度策略，建立了基于MTGP的源荷不确定集、系统容量优化配置模型和两阶段鲁棒优化调度模型，通过仿真试验可得以下结论：

1) 在源荷预测阶段充分考虑光伏出力和电、热、冷负荷间耦合关系，基于多任务学习思想确定各预测任务间相关性矩阵K_l量化描述光伏出力和电、热、冷负荷间的相关性，建立基于MTGP的源荷不确定集。实际工况中该不确定集可根据需求选择不同置信度下的置信区间，调整调度方案的保守性，提高不同需求下优化调度模型的灵活性。

2) 根据美国亚利桑那州州立大学某办公楼建筑负荷特性和所处地理环境位置，设计符合其用能需求的BIPV能源系统结构，建立容量优化配置模型并求解得到经济最优情况下系统设备和ES容量。在考虑分时电价机制下，PV和ES充分发挥削峰填谷作用，电价高峰时期利用PV出力和ES供能维持系统运行，避免高价购电，降低了系统运营成本。

3) 本文建立的两阶段鲁棒优化调度模型与确定性优化相比，充分考虑了源荷不确定性对系统稳定运行带来的影响，虽然运行成本有所增加，但补偿成本远小于确定性优化。与盒式不确定集鲁棒优化调度模型相比，本文所提优化策略深度挖掘源荷历史数据中包含的概率信息，在保证系统稳定运行的同时降低了运营成本，在提高经济性的同时考虑鲁棒性。