基于扩散模型的数据中心虚拟电厂分布鲁棒优化调度策略

Robust Optimal Scheduling Strategy for Data Center Virtual Power Plant Distribution Based on Diffusion Model

摘要: 具有强不确定性的可再生能源并网给包含数据中心的虚拟电厂的安全运行带来了巨大的挑战。为此，本文提出一种基于扩散模型的数据中心虚拟电厂分布鲁棒优化调度模型。首先，为了更准确地刻画光伏的不确定性，采用扩散模型对辐照度–温度联合场景数据进行泛化处理，提高光伏出力数据的准确性和多样性。然后，通过高斯混合聚类模型对光伏出力数据进行削减，并将得到的场景作为分布鲁棒优化集合的初始场景。再次，为进一步挖掘数据中心的需求响应潜力，建立了包含可转移批处理负载的数据中心虚拟电厂的数学模型。最后，为了求解min-max-min分布鲁棒优化问题，通过列和约束生成算法进行求解。仿真结果表明，本文所提出的基于扩散模型的光伏出力场景生成方法能更准确、更有效地刻画光伏的不确定性。此外，本文所提的分布鲁棒优化模型能平衡经济运行和鲁棒性之间的关系。

Abstract: The integration of highly uncertain renewable energy into the grid poses significant challenges to the secure operation of virtual power plants (VPPs) that include data centers. To address this issue, this paper proposes a distributionally robust optimization (DRO) scheduling model for data center VPPs based on a diffusion model. First, to more accurately characterize the uncertainty of photovoltaic (PV) power generation, a diffusion model is employed to generalize irradiance-temperature joint scenario data, enhancing both the accuracy and diversity of PV output data. Then, a Gaussian mixture clustering model is applied to reduce the dimensionality of PV output data, and the resulting scenarios are used as the initial set for DRO. Furthermore, to further explore the demand response potential of data centers, a mathematical model for data center VPPs incorporating transferable batch processing loads is developed. Finally, the min-max-min DRO problem is solved using the column-and-constraint generation algorithm. Simulation results demonstrate that the proposed diffusion model-based PV output scenario generation method in this paper can more accurately and effectively capture the uncertainty of PV. Additionally, the proposed DRO model in this paper can effectively balance the relationship between economic operation and robustness.

1. 引言

当前，人类正面临能源短缺、环境污染和气候变化等重大挑战。为满足社会对高效、低碳、环保、安全且可靠能源的日益增长需求，分布式新能源，如风能和太阳能，将成为未来能源供应的主要来源[1]。虚拟电厂(virtual power plant, VPP)是一种以电力系统为框架的协调管理系统，它集成传统能源、可再生能源、储能系统、灵活负载，并通过先进的信息通信技术实现分布式能源的聚合和调度[2]。然而，随着可再生源能在电力系统中的渗透率越来越高，其固有的不确定性严重影响VPP的稳定协同运行和能源消纳率。因此，研究准确描述可再生能源不确定性的调度控制策略，对VPP优化运行至关重要。

目前，随机优化[3]和鲁棒优化[4]是VPP中处理不确定性的两种常用方法。然而，随机优化需要获取准确的不确定性概率分布，并对大量的场景数据进行采样，保证概率分布的可信度。但所获得的概率分布往往无法准确描述实际可再生能源的真实分布，且随着场景数据的增多，导致计算时间较长[5] [6]。鲁棒优化通过利用不确定性集来描述参数的波动区间，但由于过于强调最坏环境，往往导致过于保守的求解方案[7]。新兴的分布鲁棒优化(distribution robust optimization, DRO)结合了这两种方法的优点。不需要设置概率分布的类型和参数，通过数据驱动方式建立分布模糊集，并针对最坏的概率分布进行决策。一般的DRO通常采用基于概率距离和矩信息构建模糊集，但是由此产生的NP-hard问题难以解决[8] [9]。

近年来，人工智能技术广泛应用于可再生能源场景生成方面。与传统的概率模型场景生成相比，深度学习表现出更强的泛化能力和数据表达能力[10]。文献[11]提出一种基于Wasserstein距离的生成对抗神经网络(generate adversarial networks, GANs)，用于风光出力不确定性场景的生成，相较于传统场景生成方法表现出了随机变量的时空特性。

作为体量巨大且快速增长的新兴负荷之一，全面发掘数据中心负荷调节潜力具有重要意义，数据中心虚拟电厂(data center virtual power plant, DCVPP)能够通过负荷转移实现太阳能产电量的完全消纳与分时电价的合理响应，从而提高系统各方面性能[12]。

针对上述问题，本文提出一种基于去噪扩散概率模型(denoising diffusion probabilistic models, DDPM)的数据中心虚拟电厂分布鲁棒优化运行策略。首先，为解决光伏发电的不确定性，创新性地引入DDPM来生成辐照度–温度联合场景，用于生成DRO集合的初始概率分布场景，提高数据的准确性和多样性。其次，为了进一步挖掘数据中心需求响应的潜力，本文考虑数据中心批处理负载进行转移的时间限制，建立了有效促进光伏消纳和响应分时电价的需求响应机制。最后，为提高求解效率，采用列和约束生成算法(column and constraint generation, C&CG)求解min-max-min分布鲁棒优化问题。

2. 辐照度–温度联合场景的生成和缩减

2.1. 影响光伏出力的气象特征的选取

太阳辐照是光伏电源出力的基础，光伏电源出力与太阳辐照的关系如式(1)所示：

$P_t^{PV}=P_{\max }^{PV}G\frac{1+\alpha \left(T-T_{\max }\right)}{G_{\max }}$ (1)

式中，$P_t^{PV}$ 为光伏电源功率；$G$ 为辐照度$(\text{W}/{\text{m}}^2)$ ；$P_{\max }^{PV}$ 为光伏发电的额定功率；T为光伏表面温度，取环境温度值；$\alpha$ 为功率温度系数，通常每摄氏度取−0.002~−0.005 (1/℃)；$G_{\max }$ 为额定辐照度取$1\text{kW}/{\text{m}}^2$ ；$T_{\max }$ 为额定功率下的光伏板表面温度25℃。

由式(1)得出，辐照度和温度的高低直接影响光伏出力。因此，本文选取辐照度和温度的不确定性场景集来表征光伏出力的不确定性。

2.2. 基于扩散模型的辐照度–温度场景生成

DDPM是一个基于最大似然训练范式的生成模型，能够生成较高质量数据，其核心包括正向扩散、反向去噪和模型优化[13]。基于DDPM的辐照度–温度随机场景生成结构图如图1所示。

Figure 1. DDPM-based framework diagram for irradiance-temperature stochastic scene generation

图1. 基于DDPM的辐照度–温度随机场景生成结构图

2.2.1. 正向扩散

正向扩散是逐步向原始辐照度–温度数据$x_0$ 中添加高斯噪声，使其变成标准高斯噪声$x_T=\mathcal{N}\left(0,I\right)$ 的过程。整个过程进行$T$ 步扩散，通过将高斯噪声添加到前一步数据中，以此得到每步扩散后的数据$x_t$ ，如式(2)所示：

$q\left(x_t|x_{t-1}\right)=\mathcal{N}\left(x_t;\sqrt{1-{\beta }_t}{x}_{t-1},{\beta }_{t}I\right)$ (2)

式中，$t$ 为加入噪声的次数；$x_0$ 表示未加入噪声的原始数据；$x_t$ 表示加入了t次噪声的数据；$q(\cdot )$ 表示生成数据的概率；${\left\{{\beta }_t\right\}}_{t=1}^T$ 为每一步所采用的方差，它介于0~1之间；$I$ 为单位矩阵。

扩散过程的每一步都生成一个带噪音的数据，整个扩散过程也就是一个马尔可夫链，即：

$q\left(x_{1:T}|x_0\right)={\displaystyle \prod _{t=1}^{T}q}\left(x_t|x_{t-1}\right)$ (3)

在正向扩散中，可以基于原始数据$x_0$ 在任意时间步骤$t$ 对$x_t$ 进行采样。

$\begin{array}{l}{x}_t=\sqrt{{\alpha }_t}{x}_{t-1}+\sqrt{1-{\alpha }_t}{\epsilon }_{t-1}\\ =\sqrt{{\alpha }_t}\left(\sqrt{{\alpha }_{t-1}}{x}_{t-2}+\sqrt{1-{\alpha }_{t-1}}{\epsilon }_{t-2}\right)+\sqrt{1-{\alpha }_t}{\epsilon }_{t-1}\\ =\sqrt{{\alpha }_t{\alpha }_{t-1}}{x}_{t-2}+\sqrt{1-{\alpha }_t{\alpha }_{t-1}}{\epsilon }_{t-2}\\ =\dots \\ =\sqrt{{\overline{\alpha }}_t}{x}_0+\sqrt{1-{\overline{\alpha }}_t}\epsilon \end{array}$ (4)

式中，${\overline{\alpha }}_t={\displaystyle \prod _{i=1}^t{\alpha }_i},{\alpha }_i=\sqrt{1-{\beta }_t};\epsilon =\mathcal{N}\left(0,I\right)$ 为经过重新参数化后得到的噪声。

2.2.2. 反向去噪

反向去噪过程中，去噪分布可以表示为：

$p_{\theta }\left(x_t|x_{t-1}\right)=\mathcal{N}\left(x_{t-1};{\mu }_{\theta }\left(x_t,t\right),{\sum }_{\theta }\left(x_t,t\right)\right)$ (5)

$p_{\theta }\left(x_{0:T}\right)=p\left(x_T\right){\displaystyle \prod _{t=1}^T{p}_{\theta }}\left(x_{t-1}|x_t\right)$ (6)

式中，$\theta$ 是神经网络的可学习参数。$p\left(x_T\right)$ 表示随机高斯噪声；$p_{\theta }\left(x_t|x_{t-1}\right)$ 是一个参数化高斯分布，需要训练网络来计算均值${\mu }_{\theta }\left(x_t,t\right)$ 和方差${\sum }_{\theta }\left(x_t,t\right)$ 。

通过利用神经网络使$p_{\theta }\left(x_t|x_{t-1}\right)$ 近似于真实后验分布$q\left(x_{t-1}|x_t,x_0\right)$ ，真实后验分布表示如下：

$q\left(x_{t-1}|x_t,x_0\right)=\mathcal{N}\left(x_{t-1};\tilde{\mu }\left(x_t,x_0\right),{\tilde{\beta }}_{t}I\right)$ (7)

由贝叶斯公式得：

$q\left(x_{t-1}|x_t,x_0\right)=q\left(x_t|x_{t-1},x_0\right)\frac{q\left(x_{t-1}|x_0\right)}{q\left(x_t|x_0\right)}$ (8)

由正向扩散可知：

$\{\begin{array}{c}q\left(x_t|x_{t-1},x_0\right)=\mathcal{N}\left(x_t;\sqrt{{\alpha }_t}{x}_{t-1},\left(1-{\alpha }_t\right)I\right)\\ q\left(x_t|x_0\right)=\mathcal{N}\left(x_t;\sqrt{{\overline{\alpha }}_t}{x}_0,\left(1-{\overline{\alpha }}_t\right)I\right)\\ q\left(x_{t-1}|x_0\right)=\mathcal{N}\left(x_{t-1};\sqrt{{\overline{\alpha }}_{t-1}}{x}_0,\left(1-{\overline{\alpha }}_{t-1}\right)I\right)\end{array}$ (9)

经计算化简，得到均值和方差：

$\{\begin{array}{c}{\tilde{\mu }}_t=\frac{1}{\sqrt{{\alpha }_t}}\left(x_t-\frac{1-{\alpha }_t}{\sqrt{1-{\overline{\alpha }}_t}}{\epsilon }_t\right)\\ {\tilde{\beta }}_t=\frac{1-{\overline{\alpha }}_{t-1}}{1-{\overline{\alpha }}_t}{\beta }_t\end{array}$ (10)

2.2.3. 模型优化

根据文献[14]，上述模型的训练目标可由创建变分下限(variational lower bound, VLB)表示，如式(11)所示：

$L_{\text{VLB}}=E\left[D_{KL}\left(q\left(x_T|x_0\right)\Vert p\left(x_T\right)\right)+{\displaystyle \sum _{t=2}^T{D}_{KL}}\left(q\left(x_{t-1}|x_t,x_0\right)\Vert p_{\theta }\left(x_{t-1}|x_t\right)\right)-log{p}_{\theta }\left(x_0|x_1\right)\right]$ (11)

式中，$E(\cdot )$ 为数学期望；$D_{KL}\left(q(\cdot )\Vert p(\cdot )\right)$ 为2个高斯分布的差异。

由式(11)第2项表示神经网络的训练目标为：保证每个时间步骤$t$ 的去噪分布$p_{\theta }\left(x_{t-1}|x_t\right)$ 尽可能接近正向扩散中的真实后验分布。其中，采用了KL散度(Kullback-Leibler divergence, KLD)衡量2个高斯分布的差异。

为了更好地提高样本质量和简化步骤，根据文献[14]可得出下式：

$L={\text{E}}_{t,x_0,\epsilon ~\mathcal{N}\left(0,I\right)}\left[{\Vert {\epsilon }_t-{\epsilon }_{\theta }\left(x_t,t\right)\Vert }^2\right]$ (12)

由式(12)可知，模型优化为通过对比U-Net网络生成的预测噪声与原始随机噪声之间的损失来不断学习到数据的特征，进而提高生成样本的真实性。

2.3. 基于高斯混合聚类的辐照度–温度场景削减

高斯混合模型(Gaussian mixture model, GMM)是高斯模型的扩展形式，它假设所有的样本数据均服从混合高斯分布，并通过多个高斯函数的组合来有效地逼近数据集中任意连续的概率分布，利用期望最大化算法迭代估计各高斯成分的权重、均值和协方差参数。以每个高斯分布为一个簇，将样本数据划分到距离接近的簇中，其概率分布函数如下[15]：

$\{\begin{array}{l}{p}_M\left(x\right)={\displaystyle \sum _{i=1}^k{\alpha }_i}\cdot p\left(x|{\mu }_i,{\Sigma }_i\right)\\ p\left(x\right)=\frac{1}{{\left(2\pi \right)}^{\frac{n}{2}}{\left|\Sigma \right|}^{\frac{1}{2}}}{\text{e}}^{-\frac{1}{2}{\left(x-\mu \right)}^{\text{T}}{\Sigma }^{-1}\left(x-\mu \right)}\\ {\displaystyle \sum _{i=1}^k{\alpha }_i}=1\end{array}$ (13)

式中，${\alpha }_i$ 代表第i个高斯分模型的权重，一共k个高斯分模型；$\mu$ 是n维均值向量；$\Sigma$ 是n × n维协方差矩阵。

3. DCVPP的物理模型

3.1. DCVPP总体结构

DCVPP主要由微型燃气轮机、光伏发电机组等分布式电源，电制冷机、储能装置以及数据中心负荷组成，DCVPP优化调度框架如图2所示。

3.2. 数据中心模型

1) 数据中心服务器负载分类及模型

数据中心的核心功能是处理用户提交的数据计算和存储等任务，并根据服务内容是否允许延迟处理，将任务分为交互性工作负载和批处理工作负载。批处理工作负载可以在预设时限内处理，其电能消耗视为可延迟负荷；而交互式工作负载需立即处理，电能消耗为刚性负荷，不可延迟[16]。

考虑到实际运行过程中，数据中心的服务器数量通常为一个固定值，因此电功率在某一时段数据中心的电功率有一个上限值，可表示为式(14)的形式：

$0\le P_{\text{IT,}t}\le P_{\text{IT}}^{\text{max}}$ (14)

式中，$P_{\text{IT}}^{\text{max}}$ 为数据中心最大的数据处理负荷承受值。

Figure 2. DCVPP operation framework diagram

图2. DCVPP运行框架图

数据中心通常包括IT机房和办公区两部分，IT机房的主要耗能单元以服务器为主，数据中心在t时刻的运行成本可表示为：

$F_{\text{DC}}=F_{\text{IT}}+F_{of}$ (15)

$F_{\text{IT}}=C_{\text{IT1}}+C_{\text{IT2}}$ (16)

$C_{\text{IT1}}=\left(c_{\text{grid}}+c_{\text{ma}}\right){\displaystyle {\sum }_{t=1}^T{P}_{\text{IT,t}}}$ (17)

$C_{\text{IT2}}={\displaystyle \sum _{t=1}^T{\eta }^{\text{CW}}{P}_{c,t}^{\text{CW}}}$ (18)

$P_{c,t}^{\text{CW}}={\mu }_c^{\text{CW}}{P}_{e,t}^{\text{CW}}$ (19)

式中，$F_{\text{DC}}$ 表示数据中心的运行成本；$F_{\text{IT}}$ 表示IT机房的运行成本；$F_{of}$ 表示办公区的耗电成本；$C_{\text{IT1}}$ 表示服务器的运行成本；$C_{\text{IT2}}$ 表示数据中心的冷却成本；$c_{\text{com}}$ 、$c_{\text{ma}}$ 表示VPP运行时的提供给用户的补偿费用、维护费用；$P_{\text{IT,t}}$ 为数据中心在t时刻的服务器功率。${\mu }_c^{\text{CW}}$ 表示中央空调制冷成本系数；$P_{c,t}^{\text{CW}}$ 表示中央空调制冷功率；$P_{e,t}^{\text{CW}}$ 表示中央空调吸收的电功率。

2) 数据中心中央空调模型

数据中心制冷优先考虑吸收式制冷机，若不足以满足需求，通过内部进行供冷，功率平衡如式(20)所示：

$P_{c,t}^{\text{CW}}={\kappa }_{\text{CW}}{P}_{\text{IT},t}-P_{c,t}^{\text{AC}}$ (20)

式中，${\kappa }_{\text{CW}}$ 为服务器保持正常工作状态的冷功率需求常数。

3.3. 燃气轮机模型

燃气汽轮机的运行成本可用二次曲线表示：

$C_{\text{GT}}={\displaystyle {\sum }_{t=1}^T\left[a{\left(P_{e,t}^{\text{GT}}\right)}^2+b{P}_{e,t}^{\text{GT}}+c+{\lambda }_{\text{gas,}t}{Q}_t^{\text{GT}}\right]}$ (21)

式中，$C_{\text{GT}}$ 表示燃气轮机的运行成本；$a、b、c$ 表示成本系数，该系数由燃气轮机的型号特性与气价决定；$P_{e,t}^{\text{GT}}$ 表示燃气轮机在t时刻的电输出功率；${\lambda }_{\text{gas,}t}$ 表示天然气价格。

燃气轮机电功率、热功率与燃料之间的数学模型可表述为：

$P_{e,t}^{\text{GT}}={\mu }_e^{\text{GT}}{Q}_t^{\text{GT}}$ (22)

$P_{h,t}^{\text{GT}}={\mu }_h^{\text{GT}}{Q}_t^{\text{GT}}$ (23)

式中，${\mu }_e^{\text{GT}}、{\mu }_h^{\text{GT}}$ 表示燃气锅炉产电和产热效率；$Q_t^{\text{GT}}$ 表示燃气轮机耗气量。

燃气轮机在运行中受到输出功率和爬坡约束的限制：

$P_{e,t,\min }^{\text{GT}}\le P_{e,t}^{\text{GT}}\le P_{e,t,\max }^{\text{GT}}$ (24)

$-R_{i,e}^{\text{GT}}\le P_{i,e,t}^{\text{GT}}-P_{i,e,t-1}^{\text{GT}}\le R_{i,e}^{\text{GT}}$ (25)

式中，$P_{e,t,\max }^{\text{GT}}、P_{e,t,\min }^{\text{GT}}$ 表示燃气轮机输出功率上限与下限；$R_e^{\text{GT}}$ 表示燃汽轮机的爬坡速率。

3.4. 储能设备模型

储能系统的运行成本可表示为：

$C_{\text{ESS}}={\displaystyle {\sum }_{t=1}^T\left[{\kappa }_{\text{ESS}}\left(P_{\text{ch},t}^{\text{ESS}}+P_{\text{dis},t}^{\text{ESS}}\right)\right]}$ (26)

式中，$C_{\text{ESS}}$ 为储能系统在t时刻的运行成本；$P_{\text{ch},t}^{\text{ESS}}、P_{\text{dis},t}^{\text{ESS}}$ 为储能系统在t时刻的充放电功率；${\kappa }_{\text{ESS}}$ 为储能系统折算后的单位成本系数。

在t时刻，储能系统的荷电状态的数学模型如下：

$S_{t+1}^{\text{ESS}}=S_t^{\text{ESS}}+{\eta }_{\text{ch}}\left(P_{\text{ch},t}^{\text{ESS}}\times \Delta t\right)-\frac{1}{{\eta }_{\text{dis}}}\left(P_{\text{dis},t}^{\text{ESS}}\times \Delta t\right)$ (27)

式中，${\eta }_{ch}$ 和${\eta }_{dc}$ 为储能系统的充放电效率；$S_t^{\text{ESS}}$ 表示t时刻储能系统中电能量。

为了防止储能系统过度充放电和延长使用寿命，需要满足以下约束条件：

${\text{ S}}_{\min }^{\text{ESS}}\le S_{t+1}^{\text{ESS}}\le S_{\max }^{\text{ESS}}$ (28)

式中，${\text{ S}}_{\min }^{\text{ESS}}、S_{\max }^{\text{ESS}}$ 分别为储能系统在一个调度周期内允许的最小剩余容量和最大剩余容量。

储能系统在运行中受到输出功率的限制：

$\begin{array}{l}0\le P_{\text{ch},t}^{\text{ESS}}\le P_{\text{ch},\max }^{\text{ESS}}\\ 0\le P_{\text{dis},t}^{\text{ESS}}\le P_{\text{dis},\max }^{\text{ESS}}\end{array}$ (29)

式中，$P_{\text{ch},t}^{\text{ESS}}、P_{\text{ch},\max }^{\text{ESS}}$ 分别为储能系统的充电功率和最大可充电功率；其中，$P_{\text{dis},t}^{\text{ESS}}、P_{\text{dis},\max }^{\text{ESS}}$ 为储能系统的放电功率和最大可放电功率。

为保证每个调度周期最终的剩余容量固定，使每个周期的调度方案从同一状态启动，需要满足以下约束条件：

$S_T^{\text{ESS}}=S_0^{\text{ESS}}$ (30)

式中，$S_0^{\text{ESS}}$ 为储能系统中的初始储能量；$S_T^{\text{ESS}}$ 为调度周期结束时的储能量。

3.5. 吸收式制冷机模型

$P_{c,t}^{\text{EC}}={\mu }_c^{\text{EC}}{P}_{e,t}^{\text{EC}}$ (31)

$0\le P_{c,t}^{\text{EC}}\le P_{c,\max }^{\text{EC}}$ (32)

式中，$C_{\text{EC}}$ 表示电制冷机的运行成本；$P_{c,t}^{\text{EC}}、P_{e,t}^{\text{EC}}$ 表示电制冷机产生的冷能和吸收的电能；${\mu }_c^{\text{EC}}$ 表示电制冷机制冷效率；$P_{i,c,\max }^{\text{EC}}$ 表示电制冷机单位时间输出的功率上限。

3.6. 配电网交互模型

虚拟电厂与配电网交互成本可表示为：

$C_{\text{grid}}={\displaystyle {\sum }_{t=1}^T{\lambda }_{\text{buy,t}}{P}_{\text{buy},t}}+{\displaystyle {\sum }_{t=1}^T{\lambda }_{\text{sell,t}}{P}_{\text{sell},t}}$ (33)

虚拟电厂与配电网之间的交互功率满足线路传输约束：

$\begin{array}{l}0\le P_{\text{buy},i,t}\le P_{\text{buy},i,t}^{\max }\\ 0\le P_{\text{sell},i,t}\le P_{\text{sell},i,t}^{\max }\end{array}$ (34)

式中，$P_{grid,\max }^{buy}、P_{grid,\max }^{sell}$ 分别为配电网与虚拟电厂交互功率的上限和下限。

3.7. DCVPP能量平衡模型

为保持数据中心虚拟电厂的电力和冷能供需平衡，需满足的约束条件如下：

$P_{e,t}^{\text{GT}}+P_t^{PV}+P_{\text{buy},t}+P_{\text{dis},t}^{ESS}=P_{\text{sell},i,t}+P_{\text{ch},t}^{ESS}+P_{\text{IT,}t}+P_{of,t}$ (35)

$P_{c,t}^{\text{AC}}+P_{c,t}^{\text{CW}}={\kappa }_{\text{CW}}{P}_{\text{IT},t}$ (36)

4. 两阶段分布鲁棒优化模型

4.1. 第一阶段目标函数

第一阶段考虑燃气轮机、储能设备的装机容量，各个设备的投资成本计算方法如下：

$F_{\text{GT}}=c_{\text{GT}}{S}_{\text{GT}}$ (37)

$F_{\text{ESS}}=c_{\text{ESS}}{S}_{\text{ESS}}$ (38)

式中，$F_{\text{GT}}$ 与$F_{\text{ESS}}$ 分别为燃气轮机与储能设备的安装总成本；$c_{\text{MTG}}$ 与$c_{\text{ESS}}$ 分别为燃气轮机与储能装置的单位成本；$S_{\text{GT}}$ 与$S_{\text{ESS}}$ 分别为燃气轮机与储能设备的安装容量。不同单元对应的装机容量约束如下：

$S_{\text{GT,min}}\le S_{\text{GT}}\le S_{\text{GT,max}}$ (39)

$S_{\text{ESS,min}}\le S_{\text{ESS}}\le S_{\text{ESS,max}}$ (40)

式中，$S_{\text{GT,max}}$ 和$S_{\text{GT,min}}$ 分别为微型燃气轮机装机容量的上限与下限；$S_{\text{ESS,max}}$ 和$S_{\text{ESS,min}}$ 分别为储能设备容量的上限与下限。

第一阶段目标函数如下：

$f_1=\frac{r{\left(1+r\right)}^s}{{\left(1+r\right)}^s-1}\left(F_{\text{GT}}+F_{\text{ESS}}\right)$ (41)

式中，r为折现率；s为折现年数。

4.2. 第二阶段目标函数

第二阶段的目标函数包括微型燃气轮机运行成本、配电网购售电成本、储能设备运行成本、弃光成本、CO₂排放处罚成本和数据中心运行成本。该目标函数描述如下：

$f_2\left(x\right)=F_{DC}+C_{\text{GT}}+C_{ESS}+C_{grid}+C_{Loss}+C_{C{O}_2}$ (42)

$C_{Loss}={\displaystyle {\sum }_{t=1}^T\left[c_{Loss}\left({\widehat{P}}_t^{PV}-P_t^{PV}\right)\right]}$ (43)

$C_{C{O}_2}=c_{c{o}_2}{\displaystyle {\sum }_{t=1}^T\left(k_{MTG}{P}_{EL,t}^{MTG}+k_{grid}{P}_{grid,t}^{buy}\right)}$ (44)

式中，$C_{Loss}$ 为光伏机组在t时刻的利用率成本；${\widehat{P}}_t^{PV}$ 为光伏预测出力；$C_{C{O}_2}$ 为二氧化碳排放成本；$c_{c{o}_2}$ 为CO₂排放单位罚款价格；$k_{MTG}$ 和$k_{grid}$ 分别是微型燃气轮机和电网产生的每单位电力的二氧化碳排放量。

4.3. 两阶段分布于鲁棒优化框架

第一阶段是投资阶段，确定合理的燃气轮机、储能设备、光伏机组和吸收式制冷机的装机容量，其决策变量为：$U=\left[\begin{array}{c}{S}_{\text{G}},S_{\text{ESS}}\end{array},S_{\text{PV}},S_{\text{AC}}\right]$ ；第二阶段是运行阶段，在已知第一阶段投资方案下，发现光伏不确定性后制定相应的调度计划，包括燃气轮机、吸收式制冷机输出、从配电网购电和售电、光伏实际出力、中央空调出力，其决策变量为：$V=\left[P_{e,t}^{\text{GT}},P_{c,t}^{\text{AC}},P_{\text{buy},t},P_{\text{sell},t},P_t^{PV},P_{\text{ch},t}^{\text{ESS}},P_{\text{dis},t}^{\text{ESS}},P_{c,t}^{\text{CW}}\right]$ 。

上述两阶段分布鲁棒优化模型表示如下：

$\underset{u\in U}{\text{min}}Au+\underset{p_s\in {\Omega }^p}{\text{max}}{\displaystyle \sum _{s=1}^{N_s}{p}_s}\underset{v_s\in V\left(u,{\xi }_s\right)}{min}B{v}_s$ (45)

$s.t.\{\begin{array}{c}Du\le d\\ E{v}_s=e\\ F{v}_s\le f\\ Gu+H{v}_s\le g\end{array}$ (46)

式中，${\xi }_s$ 为场景聚类筛选后第s个离散场景下的光伏出力值；$p_s$ 为聚类后情景s出现的概率；${\Omega }^p$ 为情景概率分布的集合区间，即综合范数的置信区间；n为聚类后的场景总数；$A~H$ 和$c~g$ 是常系数矩阵。

本文采用基于概率密度的分布鲁棒优化模型，并通过多离散场景的概率分布构建模糊集。$\xi$ 表示光伏出力向量，可在已获得的M个实际样本中通过场景聚类得到有限的N_s个离散场景值，如用$\left({\xi }_1,{\xi }_{\text{2}},\cdots ,{\xi }_s\right)$ 来表征光伏出力可能值，利用$\left({\xi }_1,{\xi }_{\text{2}},\cdots ,{\xi }_s\right)$ 来对日前预测的辐照度场景进行拟合，并由此获得各离散场景的初始概率分布。

由1-范数和$\infty$ -范数组成的复合范数概率置信区间表示如下：

$\mathrm{Pr}\left\{{\displaystyle \sum _{s=1}^{N_s}\left|p_s-p_{s0}\right|}\le {\theta }_1\right\}\ge 1-2N_s{e}^{-\frac{2M{\theta }_1}{N_s}}$ (47)

$\mathrm{Pr}\left\{\underset{{}_{1\le s\le N_s}}{\max }\left|p_s-p_{s0}\right|\le {\theta }_{\infty }\right\}\ge 1-2N_s{e}^{-2M{\theta }_{\infty }}$ (48)

式中：$\text{Pr}(·)$ 为求概率函数；$p_{s0}$ 为概率分布的预测值；${\theta }_1$ 和${\theta }_{\infty }$ 为1-范数和$\infty$ -范数对应约束下的概率允许偏差极限。

将(49)和(50)不等式右侧分别设为${\alpha }_1$ 和${\alpha }_{\infty }$ ，则${\alpha }_1$ 和${\alpha }_{\infty }$ 表示概率分布分别基于1-范数和$\infty$ -范数所满足的置信度，则${\theta }_1$ 和${\theta }_{\infty }$ 可分别表示为：

${\theta }_1=\frac{N_s}{2M}\ln \frac{2N_s}{1-{\alpha }_1}$ (49)

${\theta }_{\infty }=\frac{1}{2M}\ln \frac{2N_s}{1-{\alpha }_{\infty }}$ (50)

由上述公式可推导出概率分布的置信集合为[6]：

${\Omega }^p=\left\{\left\{p_s\right\}|\begin{array}{c}{p}_s\ge 0,s=1,\cdots ,N_s,{\displaystyle \sum _{s=1}^{N_s}{p}_s}=1\\ {\displaystyle \sum _{s=1}^{N_s}\left|p_s-p_{s0}\right|}\le {\theta }_1,\underset{1\le s\le N_s}{\max }\left|p_s-p_{s0}\right|\le {\theta }_{\infty }\end{array}\right\}$ (51)

4.4. 两阶段分布于鲁棒优化模型求解

本文采用列和约束生成算法(C&CG)进行求解式(45)。通过将其分解为主问题(master problem, MP)和子问题(sub problem, SP)进行交替求解，求解过程中通过不断向主问题中引入和子问题相关的变量和约束，来获得更加紧凑的目标函数的下界，满足目标迭代精度[17]。

1) 主问题

主问题是在已知初始场景概率分布$p_{s0}$ 的前提下，求解满足系统经济性的最优解，并给出式(52)的下界值，主问题可表示为：

$\underset{u\in U,{\nu }_s^w\in V\left(u,{\xi }_s\right),\eta }{\min }Au+\eta$ (52)

$s.t.\{\begin{array}{c}\eta \ge \underset{p_s\in {\Omega }^p}{\max }{\displaystyle \sum _{s=1}^{N_s}{p}_s^w}\underset{{\nu }_s^w\in V\left(u,{\xi }_s\right)}{\min }B{v}_s^{w}w=1,2,\cdots ,W\\ Gu+H{v}^w\le g\end{array}$ (53)

式中，$w$ 为迭代次数；$W$ 为最大迭代次数；$\eta$ 为引入的辅助变量。

2) 子问题

子问题在主问题给定第一阶段变量$u_w^{*}$ 下，寻找实时运行下的最恶劣概率分布，返回给主问题，并为式(54)提供上界值，子问题可表示为：

$F\left(u^{\ast }\right)=\underset{p_s\in {\Omega }^p}{\max }{\displaystyle \sum _{s=1}^{N_s}{p}_s^w}\underset{{\nu }_s\in V\left(u^{\ast },{\xi }_s\right)}{\min }B{v}_s^w,w=1,2,\cdots ,W$ (54)

由于子问题中离散场景概率值和第二阶段变量之间相互独立，故可将子问题分为两步求解，先求解子问题中内层最小值问题，然后求解子问题外层问题。具体求解流程如下：

Step 1：设定下界为$LB=0$ ，上界为$UB=\infty$ ，置迭代次数为1，收敛精度为$\alpha$ ，应用初始概率分布$p_{s0}$ ；

Step 2：求解主问题，得到第一阶段最优决策变量，更新下界值$\max \left\{LB,A{u}_w^{\ast }+{\eta }^{\ast }\right\}$ ；

Step 3：保持第一阶段变量$u_w^{\ast }$ 不变，求解子问题，获得最恶劣的场景概率分布$p_s^{*}$ 和目标函数值$F\left(u_w^{\ast }\right)$ ，并更新上界值$\min \left\{UB,A{u}_w^{\ast }+F\left(u_w^{\ast }\right)\right\}$ ；

Step 4：如果$UB-LB\le \alpha$ ，停止迭代，返回最优解；反之，更新最坏概率值，更新主问题中的第二阶段决策变量及相关约束，重复以上步骤。

5. 算例分析

本文以图2所示的含数据中心的虚拟电厂为算例，验证所提出的分布鲁棒优化模型和求解算法的有效性，并使用MATLAB R2023和GUROBI对优化程序进行求解，表1是本文设备运行参数。辐照度–温度预测场景及其光伏预测出力如图3所示。

Table 1. Equipment operating parameters

表1. 设备运行参数

参数	数值	参数	数值	参数	数值
$P_{IT}^{\text{max}}/\text{kW}$	1200	${\eta }_{\text{WH}}$	0.35	r	5
${\kappa }_{\text{ec}}$	0.3	${\kappa }_{ESS}$	0.02	s	20
${\kappa }_{\text{CW}}$	3.3	$P_{grid,\max }^{buy}$	1000	$k_{MTG}$	0.49
a	0.004	${\eta }_{ch}$	0.95	$k_{grid}$	0.8
b	1.2	${\eta }_{dc}$	0.95	$R_{EL}^{\text{MTG}}/\text{kW}$	120
c	39	$c_{\text{MTG}}$	1200	$c_{\text{ESS}}$	800

(a) 辐照度–温度 (b) 光伏出力

Figure 3. Irradiance-temperature and photovoltaic output prediction diagram

图3. 辐照度–温度和光伏出力预测图

5.1. 扩散模型性能评估

1) DDPM模型训练评估

本文算例中的辐照度和温度的数据来源于澳大利亚的DKASC测试平台。为了研究辐照度和温度的联合场景，选取了实验室公开数据集中0.8 MW的光伏发电单元在2016~2022年的辐照度和温度数据，采样时间间隔5 min。数据集80%用于训练，20%用于测试。

为验证所提DDPM算法在场景生成上的准确性和多样性，本文将DDPM方法和GAN方法进行对比，DDPM和GAN在训练过程中的损失函数变化曲线如图4(a)和图4(b)所示。此外，本文将DDPM和GAN训练速度进行了对比，选取Epoch = 500进行对比，DDPM训练时间5小时，而GAN训练时间8小时。结果表明，DDPM训练速度比GAN有较大提升。

2) DDPM和GAN的性能指标对比

① 最大化均值差异(maximum mean discrepancy, MMD)距离度量[18]：

最大平均差异(MMD)是希尔伯特空间中两个分布之间差异的度量。它可以用来度量生成的数据集与原始数据集之间的距离。MMD表示某固定核函数$k$ 的$P_d$ 与$P_g$ 之差，定义如下：

$MM{D}^2\left(P_d,P_g\right)=E_{x_d,{x^{\prime }}_d~P_d,x_g,{x^{\prime }}_g~P_g}\left[k\left(x_d,{x^{\prime }}_d\right)-2k\left(x_d,x_g\right)+k\left(x_g,{x^{\prime }}_g\right)\right]$ (55)

式中，$k(\cdot )$ 表示核函数；$x_d、x_g$ 表示原始数据集和生成数据集的分布；${x^{\prime }}_d、{x^{\prime }}_g$ 表示$x_d、x_g$ 映射到希尔伯特空间中的分布。MMD越小，两个分布越接近。

DDPM和GAN在训练过程中的MMD的变化过程如图5(a)和图5(b)所示。

② Fréchet inception distance (FID)距离度量：

FID是通过计算两个分布之间的Fréchet距离来衡量生成模型和真实数据分布之间的差异。它考虑了两个分布之间的均值和协方差矩阵，可以更好地描述两个分布之间的差异，定义如下：

$FlD\left(P_x,P_g\right)=\Vert {\mu }_x-{\mu }_g\Vert +Tr\left[{\sigma }_x^2+{\sigma }_g^2-2{\left({\sigma }_x^2{\sigma }_g^2\right)}^{\frac{1}{2}}\right]$ (56)

式中：${\mu }_x$ 和${\mu }_g$ 分别为真实数据和生成样本分布的平均值；${\sigma }_x^2$ 和${\sigma }_g^2$ 分别为实际数据和生成样本分布的协方差。FID越小，两个分布越接近。

DDPM和GAN在训练过程中的FID的变化过程如图6(a)和图6(b)所示。

(a) U-Net损失函数图 (b) GAN损失函数图

Figure 4. Comparison diagram of DDPM and GAN loss functions

图4. DDPM和GAN损失函数对比图

通过MMD训练过程的变化对比图可知，GAN的MMD最终收敛于0.5，并且几乎没有波动性，说明GAN生成的数据准确性较好但是在生成数据的多样性方面较差。DDPM最终大部分数据在于0~0.5之间，在准确性方面更好，并且曲线波动性较强，说明DDPM在生成数据的多样性方面较生成数据的多样性方面更好。同理，通过MMD训练过程的变化对比图可知，GAN最终收敛于10并且几乎没有波动。DDPM最终数据大部分存在于5~10之间，并且波动性更大。

综上所述，DDPM在生成数据的准确性和多样性方面比GAN有更大的优势。

(a) U-Net (b) GAN

Figure 5. Comparison of MMD changes during training

图5. 训练过程MMD变化对比

(a) U-Net (b) GAN

Figure 6. Comparison of FID changes during training

图6. 训练过程FID变化对比

3) 生成场景和聚类结果

图7显示通过DDPM生成的场景和通过GMM聚类后获得的典型场景。

Figure 7. Generate scene and clustering results

图7. 生成场景和聚类结果

5.2. 对比结果分析

5.2.1. 与确定性方法对比分析

将所提出的DRO方法与传统确定性模型的调度结果进行对比。分布鲁棒优化M = 2000，${\alpha }_1={\alpha }_{\infty }=0.99$ ，对比结果如表2所示。

从表2中可知，采用DRO方法相比确定性方法具有更好的收益，并且提高了清洁能源的消纳率。这是因为考虑光伏不确定性的波动后，系统会更加保守，当出现不确定性波动较大的情况时，不至于出现较多的弃电情况。另外，由于DRO方法中也考虑了场景概率分布的不确定性，即双重不确定性，因而获得了最恶劣的概率分布期望，表现出更优的经济性能。这说明采用DRO方法进行决策在清洁能源消纳和经济性方面更具有优势。

5.2.2. 与其他不确定性方法对比分析

将所提出的DRO方法与传统的随机优化和鲁棒优化进行对比。随机优化采用典型场景和相应的概率分布；鲁棒优化采用盒型不确定性集，波动区间设置为20%；分布鲁棒优化M = 2000，${\alpha }_1={\alpha }_{\infty }=0.99$ ，对比结果如表2所示。

从表2中可知，DRO优化的含数据中心的虚拟电厂的运营成本低于鲁棒优化，而高于随机优化。这是因为随机优化的概率分布没有考虑最恶劣的情况，因此保守性较差；鲁棒优化过于关注恶劣情况，导致保守性太强，在光伏的消纳和经济性上较差。DRE方法在经济性和保守性上达到较好的平衡，最大程度降低运营成本的同时，也提高了清洁能源的消纳率，在应对不确定性规划方面更具优势。

Table 2. Algorithm comparison results

表2. 算法对比结果

优化结果		优化方法
		确定性算法	随机优化	鲁棒优化	本文所提方法
成本和 利润/(10^3/¥)	燃气轮机成本	8.1266	7.3020	8.4596	7.9887
	买电利润	3.6243	3.5882	3.7566	3.6972
	卖电利润	1.1101	1.2129	1.1506	1.0998
	储能设备成本	0.1605	0.1563	0.1663	0.1370
	吸收式制冷机成本	1.1349	1.1757	1.1398	1.1673
	需求响应成本	0.8862	0.7436	0.9186	0.8149
	碳排放成本	1.4742	1.5091	1.5280	1.4976
	弃光成本	0.7905	0.5915	0.8193	0.5018
	总运营成本	15.0727	13.8418	15.6228	14.6982
弃电率/%	弃光率	20	10.86	13.66	8.16

5.2.3. 模糊集参数对比分析

两阶段分布鲁棒优化模型中共有四个参数：场景生成数M、场景聚类数K、1-范数约束以及$\infty$ -范数约束。下面对各参数对调度结果的影响进行分析，设定1-范数和$\infty$ -范数为0.95，研究场景生成数M以及场景聚类数K对调度结果的影响。对比结果如表3所示。

Table 3. Comparative analysis of parameters K and M

表3. 参数K、M对比分析

M	调度成本/(10^3/¥)
	K = 5	K = 10	K = 15	K = 20
2000	14.5431	14.6982	14.7215	14.7934
5000	14.4965	14.5624	14.6424	14.7631
8000	14.4001	14.4863	14.5862	14.6934

从表3可得，随着场景生成数M的逐渐增大，DCVPP调度成本逐渐减小，随着场景缩减数K的逐渐增大，调度成本逐渐增大。当场景生成数M逐渐增大时，随着场景数的增多，会使得缩减场景的初始概率越来越接近于真实场景分布，减小在优化过程中场景的偏移程度，从而降低了不确定性对DCVPP优化调度的影响，使得成本降低，提高系统调度经济性。当缩减场景数增多时，首先增加了系统的计算负担，其次随着缩减场景数目的增多，会使得缩减场景对最恶劣光伏出力场景的刻画逐渐变得清晰，最恶劣光伏出力场景可成为典型的缩减场景，系统保守程度会随着缩减场景数目的变多而加深，从而因为提高系统的鲁棒性而降低系统调度经济性水平。

6. 结论

针对包含数据中心的虚拟电厂的调度问题，本文提出了一种基于扩散模型的两阶段分布鲁棒优化模型，通过算例验证，得出如下结论：

1) 本文为了更准确地描述光伏的不确定性，提出了一种基于DDPM的场景生成方法，其不仅避免了GAN的梯度消失和模态崩溃问题，使训练更加稳定，而且在生成场景的准确性和多样性方面更有优势。

2) 本文从数据中心运营商的角度出发，考虑虚拟电厂的优化调度问题，建立考虑时间限制的批处理负载转移模型，实现了调峰填谷以及充分响应了分时电价，降低了运营成本。

3) 本文提出的DRO方法与随机规划和鲁棒规划相比，DRO方法兼具随机规划的经济性和鲁棒规划的保守性，在经济性和保守性之间达到了良好的均衡。

基金项目

国家自然科学基金(62473151, 61973116)、国家重点研发计划(2019YFB1505400)、中央高校基本科研业务费专项资金(2023JC001)资助项目。

NOTES

^*通讯作者。