天然气的冷热电联供系统（combined cooling, heating and power，CCHP）是连接电网与气网的耦合系统，也是区域综合能源系统（regional integrated energy system，RIES）中最具发展前景的一种运营模式^[1]. 目前，RIES的研究多以优化不同效益目标，得到系统各设备的运行策略为主. Wei等^[2]构建电转气的峰值负荷转移模型，从理论上证明了电气耦合系统具有较好的削峰填谷的效果；Guandalini等^[3]对电气耦合系统进行了评价，结果证明，该系统可以提高可再生能源的可调度性；张儒峰等^[4]提出合理利用弃风的电-气综合能源系统，实现互联系统之间的双向耦合；Qu等^[5]利用电转气实现电力系统与天然气系统的双向能量流动，是促进风电消纳平滑功率需求的有效途径. 上述文献均未考虑发、用电波动性对系统带来的收益损失风险. 张虹等^[6-8]引入CVaR计算一定收益下系统要承担的收益风险，但均将CVaR转化为离散情况下最差CVaR进行求解，该方法结果的主观性强. 上述研究均仅考虑单维不确定变量，不适用于同时考虑多种不确定关系耦合下的建模与分析.

综合需求响应利用冷热负荷的惯性特征，是平衡新电改下各市场主体利益诉求的绝佳手段^[9]，但RIES中考虑需求响应的调度方法较少涉及. 张虹等^[6]让需求侧互动资源主动提供用电意愿，根据系统调度灵活选择用电行为；王文超等^[9]将电价型需求响应应用于系统优化运行中；徐业琰等^[10]通过电价型、激励型和博弈方法协同作用，实现对用户侧的联合调度. 来自于荷、源双侧（如风电、光伏、负荷等）的多重不确定性是RIES运行时面临的主要挑战. 在描述发电与用电不确定性上，场景法^[11]、点估计法^[12]、随机机会约束规划^[13]、模型预测^[14]都有较好的应用，但随机法与点估计法均需要实际中的大量样本数据，且场景法结果受场景个数的制约，机会约束规划难以保证求解效率与精度.

鉴于以上分析，本文建立基于Copula-RCVaR的区域综合能源系统两阶段鲁棒博弈优化调度模型. Copula-RCVaR模型能够对多个不确定变量耦合、不同决策需求下系统的收益损失风险进行分析与评估. 考虑综合需求响应，利用CCHP机组和虚拟能量厂（virtual energy plant，VEP）平抑RIES内发电、用电波动性. 采用鲁棒理论，对系统内不确定变量建立不确定性合集，剖析不确定变量与系统经济性、保守性的动态相依关系，探索在不同决策需求下最经济可靠的调度方案. 针对模型特点，利用基于滤子技术的多目标鲸鱼算法进行求解. 以修改的IEEE33节点配电网与CCHP系统耦合形成RIES为例，验证模型能够在保证安全稳定的前提下，平衡各层主体利益，实现电力经济的可持续性发展.

如图1所示为RIES结构及运行方式示意图. 在经典CCHP系统组成的RIES中，加入能量集线器与冷/热/电储能装置. 燃汽轮机是系统中的主要源动设备，发电量与RIES在能量交易中心向上级电网的购电量（包括在日前市场与实时市场的购电）共同承担系统负荷用电，通过余热转换装置与锅炉向系统内用户提供热负荷需求，系统的冷负荷需求由电制冷机和吸收式制冷机提供，电制冷机由电能驱动，吸收式制冷机由热能驱动.

系统中，内燃气轮机和锅炉运行所需的天然气由RIES在能量交易中心向上级气网购得（仅在实时市场）. 为了减少天然气的消耗，在系统中加入可再生能源电站（图1中的风电场），承担系统内部分电负荷与热负荷需求，可再生能源发电不接受调控且不计发电成本. RIES在日前市场向上级电网购买电量，能量盈余或亏空通过实时市场与上级电、气网的能量交换，调控CCHP机组组合出力、VEP、各能量转换装置得到平衡. 如图1所示，VEP包括各储能系统与各种类型的可控负荷. 其中可控负荷根据特性分为以下4类^[15]. 1）常规负荷（CL），具有较大的随机性与波动性，且不可调控. 2）迎峰负荷（LSI），切负荷量较低，一般为该类型总量的15%，补偿价格指数较高. 3）避峰负荷（LSII），该类型负荷用电灵活性较大，切负荷量较高，为总量的30%，且补偿价格指数较低. 以上3种类型仅有电负荷CL-e/LSI-e/LSII-e. 4）可转移负荷（TL），在不影响使用舒适度的前提下转移，补偿价格系数较低，但转移前、后的负荷总量不变，分为TL-h/TL-c，表示热/冷负荷.

该模型将RIES与VEP作为电力系统中不同的市场主体，针对运营体系及特点，采用双层多目标鲁棒优化对混合系统进行建模. 其中RIES位于上层，VEP位于下层. 优化时，先由RIES向VEP发送调度计划，VEP在满足自身运行约束的前提下调控管辖内的可控资源（各储能系统、可控负荷）对该计划实行初步响应；将自身优化的结果反馈至上层，RIES根据反馈结果进一步调整计划. 过程中，上、下两层信息互相更新与传递，在尽可能满足各系统电力需求的前提下，经济性、社会性最好. CCHP机组与VEP的参与可平抑发电与用电的波动性，将盈余电量在实时市场较稳定地外送，使RIES获利，该运营模式在一定程度上可以提高可再生能源的并网能力. 由于冷/热网中的冷/热惯性，使得冷/热负荷中的不确定性能够被各自传输管道中的管存能力缓解^[16]，模型只考虑发、用电不确定性.

CCHP装置互相耦合，与上级电、气网共同实现对RIES能源的供应，各装置按如下方式建模.

1）燃汽轮机. 出力与耗气量为二次函数.

(1) $G_t^{{\rm{GT}}} = {y_{{\rm{G}}{{\rm{T}}_t}}} \left( {{a_1}{{\left( {P_{{\rm{G}}{{\rm{T}}_t}}^{\rm{e}}} \right)}^2} + {b_1}\left( {P_{{\rm{G}}{{\rm{T}}_t}}^{\rm{e}}} \right) + {c_1}{y_{{\rm{G}}{{\rm{T}}_t}}}} \right).$

式中：a₁、b₁、c₁为燃汽轮机的耗气常数； ${y_{_{{{\rm{GT}}_{_t}}}}} $、 $P_{{\rm{G}}{{\rm{T}}_t}}^{\rm{e}} $、 $G_t^{{\rm{GT}}} $为t时刻燃气轮机的运行状态变量、出力和耗气量，其中 ${{{y}}_{{\rm{G}}{{\rm{T}}_{{t}}}}} $=1为运行. 燃气轮机应满足相关的运行约束如下.

(2) ${y_{{\rm{G}}{{\rm{T}}_t}}}P_{{\rm{G}}{{\rm{T}}_{\min }}}^{\rm{e}} \leqslant P_{{\rm{G}}{{\rm{T}}_t}}^{\rm{e}} \leqslant {y_{{\rm{G}}{{\rm{T}}_t}}}P_{{\rm{G}}{{\rm{T}}_{\max }}}^{\rm{e}},\;\forall t.$

(3) $ \left. \begin{aligned} & P_{{\rm{G}}{{\rm{T}}_t}}^{\rm{e}} - P_{{\rm{G}}{{\rm{T}}_{(t - 1)}}}^{\rm{e}} \leqslant {y_{{\rm{G}}{{\rm{T}}_t}}}R_{{\rm{G}}{{\rm{T}}_{\rm{U}}}}^{\rm{e}} + \left( {1 - {y_{{\rm{G}}{{\rm{T}}_{(t - 1)}}}}} \right)P_{{\rm{G}}{{\rm{T}}_{\max }}}^{\rm{e}} + \\& \qquad P_{{\rm{G}}{{\rm{T}}_{\min }}}^{\rm{e}}\left( {{y_{{\rm{G}}{{\rm{T}}_t}}} - {y_{{\rm{G}}{{\rm{T}}_{(t - 1)}}}}} \right),\;t \geqslant 2;\\& P_{{\rm{G}}{{\rm{T}}_{(t - 1)}}}^{\rm{e}} - P_{{\rm{G}}{{\rm{T}}_t}}^{\rm{e}} \leqslant {y_{{\rm{G}}{{\rm{T}}_t}}}R_{{\rm{G}}{{\rm{T}}_{\rm{D}}}}^{\rm{e}} + \left( {1 - {y_{{\rm{G}}{{\rm{T}}_{(t - 1)}}}}} \right)P_{{\rm{G}}{{\rm{T}}_{\max }}}^{\rm{e}} + \\& \qquad P_{{\rm{G}}{{\rm{T}}_{\min }}}^{\rm{e}}\left( {{y_{{\rm{G}}{{\rm{T}}_{(t - 1)}}}} - {y_{{\rm{G}}{{\rm{T}}_t}}}} \right),\;t \geqslant 2. \end{aligned} \right\} $

(4) $\left. {\begin{gathered} - {y_{{\rm{G}}{{\rm{T}}_{\left( {t - 1} \right)}}}} + {y_{{\rm{G}}{{\rm{T}}_t}}} - {y_{{\rm{G}}{{\rm{T}}_k}}} \leqslant 0,\;\forall t;\;t \leqslant k \leqslant T_{{\rm{GT}}}^{{\rm{on}}} + t - 1 . \\ {y_{{\rm{G}}{{\rm{T}}_{\left( {t - 1} \right)}}}} - {y_{{\rm{G}}{{\rm{T}}_t}}} + {y_{{\rm{G}}{{\rm{T}}_k}}} \leqslant 1,\;\forall t;\;t \leqslant k \leqslant T_{{\rm{GT}}}^{{\rm{off}}} + t - 1 . \\ \end{gathered} } \right\}$

式（2）为发电功率约束，式（3）为爬坡约束，式（4）为最小启停时间约束. 式中： $P_{{\rm{G}}{{\rm{T}}_{\min }}}^{\rm{e}} $、 $P_{{\rm{G}}{{\rm{T}}_{\max }}}^{\rm{e}} $为出力上、下界限； $R_{{\rm{G}}{{\rm{T}}_{\rm{U}}}}^{\rm{e}} $、 $R_{{\rm{G}}{{\rm{T}}_{\rm{D}}}}^{\rm{e}} $分别为向上和向下最大爬坡功率， $T_{{\rm{GT}}}^{{\rm{on}}} $、 $T_{{\rm{GT}}}^{{\rm{off}}} $分别为最小开机和停机时间.

2） 余热回收装置. 该装置输出热量与燃汽轮机的出力有关：

(5) $P_{{\rm{WH}}{{\rm{R}}_t}}^{\rm{h}} = {y_{{\rm{G}}{{\rm{T}}_t}}} \left( {{a_2}{{\left( {P_{{\rm{G}}{{\rm{T}}_t}}^{\rm{e}}} \right)}^2} + {b_2}\left( {P_{{\rm{G}}{{\rm{T}}_t}}^{\rm{e}}} \right) + {c_2}{y_{{\rm{G}}{{\rm{T}}_t}}}} \right).$

式中：a₂、b₂、c₂为耗量系数， $P_{{\rm{WH}}{{\rm{R}}_t}}^{\rm{h}} $为转换的可用热量.

3） 电制冷机、吸收式制冷机和锅炉.

(6) $\left. {\begin{split} & P_{{\rm{AB}}{{\rm{S}}_t}}^{\rm{c}} = {\eta _{{\rm{ABS}}}}P_{{\rm{AB}}{{\rm{S}}_t}}^{\rm{h}},\\ & P_{{\rm{AS}}{{\rm{R}}_t}}^{\rm{c}} = {\eta _{{\rm{ASR}}}}P_{{\rm{AS}}{{\rm{R}}_t}}^{\rm{e}},\\ & P_{{{\rm{B}}_t}}^{\rm{h}} = {\eta _{\rm{B}}}G_t^{\rm{B}}. \end{split}} \right\}$

式中：η_ABS、η_ASR分别为吸收式制冷机与电制冷机的效率，η_B为锅炉热效率， $P_{{\rm{AS}}{{\rm{R}}_t}}^{\rm{c}} $、 $P_{{\rm{AS}}{{\rm{R}}_t}}^{\rm{e}} $分别为电制冷机的制冷量与耗电量， $P_{{\rm{AB}}{{\rm{S}}_t}}^{\rm{c}} $、 $P_{{\rm{AB}}{{\rm{S}}_t}}^{\rm{h}} $分别为吸收式制冷机的制冷量与吸热量， $P_{{{\rm B}_t}}^{\rm{h}}$、 $G_{{t}}^{\rm{B}} $分别为锅炉的产热量与耗气量. 此外，应满足相关运行约束如下：

(7) $\left. { \begin{split} & 0\leqslant {P}_{{\rm{WHR}}_{t}}^{\rm{h}}\leqslant {{P}_{\rm{WHR}_\mathrm{max}}^{\rm{h}}},\\& 0\leqslant {P}_{{\rm{ABS}}_{t}}^{\rm{h}}\leqslant {P}_{{\rm{ABS}}_{\mathrm{max}}}^{\rm{h}},\\& 0\leqslant {P}_{{\rm{ASR}}_{t}}^{\rm{e}}\leqslant {P}_{{\rm{ASR}}_{\mathrm{max}}}^{\rm{e}},\\& 0\leqslant {P}_{{\rm{B}}_{t}}^{\mathrm{h}}\leqslant {P}_{{\rm{B}}_{\mathrm{max}}}^{\mathrm{h}},\;\forall t.\end{split}} \right\}$

　　1）目标函数1. RIES运营利润最高为

(8) $ \begin{split}&\!\!\!\!\!\!\max \;{f_{1.1}} = \sum\limits_{t = 1}^T {\left[ {C_t^{{\rm{S - e}}}\left( {P_{{{\rm{L}}_t}}^{\rm{e}} - \sum\limits_{i = 1}^{N_{\rm{e}}^{{\rm{VEP}}}} {P_{i,t}^{{\rm{VEP - e}}}} } \right) + } \right.} \\ &C_t^{{\rm{S - h}}}\left( {P_{{{\rm{L}}_t}}^{\rm{h}} - \sum\limits_{i = 1}^{N_{\rm{h}}^{{\rm{VEP}}}} {P_{i,t}^{{\rm{VEP - h}}}} } \right) + C_t^{{\rm{S - c}}}\left( {P_{{{\rm{L}}_t}}^{\rm{c}} - \sum\limits_{i = 1}^{N_{\rm{c}}^{{\rm{VEP}}}} {P_{i,t}^{{\rm{VEP - c}}}} } \right) - \\ &{C_{\rm{b}}^{\rm{e}}{\eta _t}P_{{{\rm{L}}_t}}^{{\rm{e - s}}} - {C^{\rm{g}}}\left( {G_t^{\rm{B}} + G_t^{{\rm{GT}}}} \right) + C_t^{{\rm{re}}}P_t^{{\rm{grid}}}} \Bigg]{\rm{ - }}\\ &{C_{{\rm{VEP}}}} - {C_{{\rm{conv}}}} - {C_{{\rm{ek}}}}. \end{split} $

式中：η_t为日前市场电量购买比例；C_t^S-e、C_t^S-h、C_t^S-c分别为电、热、冷能出售价格； $C_{\rm{b}}^{\rm{e}} $、 $C_{{t}}^{\rm{re}} $分别为日前、实时市场向上级电网的购电价格；C^g为向上级气网的购气价格； $P_{{{\rm{L}}_t}}^{\rm{e}} $、 $P_{{{\rm{L}}_t}}^{\rm{h}} $、 $P_{{{\rm{L}}_t}}^{\rm{c}} $分别为t时刻电、热、冷负荷； $P_{{{\rm{L}}_t}}^{{\rm{e - s}}} $为电负荷的预测值， $ P_{i,t}^{{\rm{VEP - e}}}$、 $P_{i,t}^{{\rm{VEP - h}}} $、 $P_{i,t}^{{\rm{VEP - c}}} $为RIES对第i个VEP下达的调度计划，上标e、h、c表示VEP类型， $P_{i,t}^{{\rm{VEP - e}}} $>0表示向系统注入能量； $P_t^{{\rm{grid}}} $为系统与上级电网之间的电量交换，P_t^grid>0表示RIES向上级电网售电，反之为购电；C_VEP为VEP运行成本，由下层模型计算得出返回至上层；C_conv、C_ek为转换装置运行成本及旋转备用成本，

(9) ${C_{{\rm{conv}}}} = \sum\limits_{t = 1}^T {\left( {{S_{{\rm{WHR}}}}P_{{\rm{WH}}{{\rm{R}}_t}}^{\rm{h}} + {S_{{\rm{ABS}}}}P_{{\rm{AB}}{{\rm{S}}_t}}^{\rm{c}} + {S_{{\rm{ASR}}}}P_{{\rm{AS}}{{\rm{R}}_t}}^{\rm{c}}} \right)} ,$

(10) ${C_{{\rm{ek}}}} = \sum\limits_{t = 1}^T {\left\{ {\sum\limits_{i = 1}^{{N_{\rm{W}}}} {{K_{{{\rm{W}}_i}}}{\gamma _{{{\rm{W}}_i}}}\Delta P_{{{\rm{W}}_{i,t}}}^{\rm{u}} + \sum\limits_{i = 1}^{N_{\rm{L}}^{\rm{e}}} {{K_{{{\rm{L}}_i}}}{\gamma _{{{\rm{L}}_i}}}\Delta P_{{{\rm{L}}_{i,t}}}^{{\rm{e - u}}}} } } \right\}} .$

式中：S_WHR、S_ABS、S_ASR分别为余热回收装置、吸收式制冷机、电制冷机的成本系数，N_W、 $N_{\rm{L}}^{\rm{e}} $分别为风电站及常规电负荷总数， ${K_{{{\rm{W}}_i}}} $/ ${K_{{{\rm{L}}_i}}} $、 ${\gamma _{{{\rm{W}}_{{i}}}}}$/ ${\gamma _{{{\rm{L}}_{{i}}}}}$、Δ $P_{{{\rm{W}}_{i,t}}}^{\rm{u}} $/Δ $P_{{{\rm{L}}_{i,t}}}^{\rm{e-u}} $分别为风电/常规电负荷的旋转备用惩罚系数、功率偏差系数及各功率偏差上限.

2）目标函数2. RIES的收益损失风险最小值min f_1.2，详细见2章的风险能量管理模型，此处不赘述.

3）目标函数3. RIES与VEP调度偏差最小.

(11) $\begin{split} \min\; {f_{1.{\rm{3}}}} =& \sum\limits_{t = 1}^T \left( \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{e}}} {\left| {P_{i,t}^{{\rm{VEP - e}}} - P_{i,t}^{{\rm{VEP - }}{{\rm{e}}^*}}} \right|} + \right. \\ & \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{h}}} {\left| {P_{i,t}^{{\rm{VEP - h}}} \!-\! P_{i,t}^{{\rm{VEP - }}{{\rm{h}}^*}}} \right| }\!+\! \left. \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{c}}} {\left| {P_{i,t}^{{\rm{VEP - c}}}\! -\! P_{i,t}^{{\rm{VEP - }}{{\rm{c}}^*}}} \right|} \right) . \end{split}$

式中： $P_{i,t}^{{\rm{VEP}} - {\rm{e^*}}}$、 $P_{i,t}^{{\rm{VEP}} - {\rm{h^*}}}$、 $P_{i,t}^{{\rm{VEP}} - {\rm{c^*}}}$分别为通过下层模型优化返回至上层，第i个虚拟电、热、冷厂在t时刻的调度功率； $N_{\rm{v}}^{\rm{e}} $、 $N_{\rm{v}}^{\rm{h}} $、 $N_{\rm{v}}^{\rm{c}} $分别为各虚拟电、热、冷厂的个数.

上层模型除式（2）~（4）、（7）外，还需满足如下约束.电集线器平衡约束为

(12) $\begin{split} & {\eta _t}P_{{{\rm{L}}_t}}^{{\rm{e - s}}} + P_t^{{\rm{VEP - e}}} + \sum\limits_{i = 1}^{{N_{\rm{W}}}} {P_{{{\rm{W}}_{i,t}}}^{\rm{e}}} + {y_{{\rm{G}}{{\rm{T}}_t}}}P_{{\rm{GT}}{}_t}^{\rm{e}} - \\& \sum\limits_{i = 1}^{N_{\rm{L}}^{\rm{e}}} {P_{{{\rm{L}}_{i,t}}}^{\rm{e}}} - p_t^{{\rm{grid}} } - P_{{\rm{AS}}{{\rm{R}}_t}}^{\rm{e}} = 0,\;\forall t. \end{split}$

热集线器转换约束为

(13) $P_t^{{\rm{VEP - h}}} + p_{{{\rm{B}} _t}}^{\rm{h}} + P_{{\rm{WH}}{{\rm{R}}_t}}^{\rm{h}} \geqslant P_{{L_t}}^{\rm{h}}{\rm{ + }}P_{{\rm{AB}}{{\rm{S}}_t}}^{\rm{h}},\;\forall t.$

冷集线器转换约束为

(14) $P_{{\rm{AS}}{{\rm{R}}_t}}^{\rm{c}} + P_{{\rm{AB}}{{\rm{S}}_t}}^{\rm{c}} + P_t^{{\rm{VEP - c}}} \geqslant P_{{{\rm{L}}_t}}^{\rm{c}},\;\forall t.$

为了防止RIES与上级电网之间的联络线功率毛刺过多，使其能够运行平稳，将 $P_t^{{\rm{grid}}} $离散成10的整数倍，设置上、下功率界限为 $P_{\rm{max}}^{{\rm{grid}}} $、 $P_{\rm{min}}^{{\rm{grid}}} $，最大爬升功率为120 kW，最小保持功率时间为2 h.

VEP将RIES下达的调用计划分解至各个可控单元上，使得两层之间的调度计划偏差最小，VEP达到最大的经济效益与社会效益.

1）目标函数1. 调度计划偏差最小.

(15) $\left. {\begin{split} & \min \;{f_{2.1}} = \sum\limits_{t = 1}^T \left( \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{e}}} {\left| {P_{i,t}^{{\rm{VEP - }}{{\rm{e}}^*}} - P_{i,t}^{{\rm{VEP - e}}}} \right|} + \right. \\ & \left. \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{h}}} {\left| {P_{i,t}^{{\rm{VEP - }}{{\rm{h}}^*}} - P_{i,t}^{{\rm{VEP - h}}}} \right| + } \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{c}}} {\left| {P_{i,t}^{{\rm{VEP - }}{{\rm{c}}^*}} - P_{i,t}^{{\rm{VEP - c}}}} \right|} \right) ; \\ & P_{i,t}^{{\rm{VEP - }}{{\rm{e}}^{\rm{*}}}} = P_{{\rm{LS}}{{\rm{{\rm I}}}_{i,t}}}^{\rm{e}} + P_{{\rm{LSI}}{{\rm{I}}_{i,t}}}^{\rm{e}} + P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{e}}, \\ & P_{i,t}^{{\rm{VEP - }}{{\rm{h}}^{\rm{*}}}} = P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{h}} + P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{h}} ,\\ & P{_{i,t}^{{\rm{VEP - c^{*}}}}{}} = P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{c}} + P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{c}} . \end{split} } \right\}$

式中： $P_{i,t}^{{\rm{VEP}} - {\rm{e}}} $、 $P_{i,t}^{{\rm{VEP}} - {\rm{h}}} $、 $P_{i,t}^{{\rm{VEP}} - {\rm{c}}} $由上层模型优化所得并传递至下层.

2）目标函数2. 经济效益最好，调度成本最小.

(16) $\left. {\begin{split} &\!\!\!\!\!\!\!\!\!\! \min\; {f_{2.{\rm{2}}}} = \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{e}}} {C_{{\rm{VE}}{{\rm{P}}_i}}^{\rm{e}}} {\rm{ + }}\sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{h}}} {C_{{\rm{VE}}{{\rm{P}}_i}}^{\rm{h}} + } \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{c}}} {C{{_{{\rm{VEP}}}^{\rm{c}}}_{_i}}} ; \\ & C_{{\rm{VE}}{{\rm{P}}_i}}^{\rm{e}} \!=\! \sum\limits_{t = 1}^T {\left( {\xi _{{\rm{LS{\rm I}}}}^{\rm{e}} P{{_{{\rm{LS{\rm I}}}}^{\rm{e}}}_{_{i,t}}} + \xi _{{\rm{LSII}}}^{\rm{e}} P_{{\rm{LSI}}{{\rm{I}}_{i,t}}}^{\rm{e}} + \xi _{{\rm{ESS}}}^{\rm{e}} \left| {P_{{\rm{ES}}{{\rm{S}}_{_{i,t}}}}^{\rm{e}}} \right|} \right)} , \\ &C_{{\rm{VE}}{{\rm{P}}_i}}^{\rm{h}} = \sum\limits_{t = 1}^T {\left( {\xi _{{\rm{LT}}}^{\rm{h}} \left| {P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{h}}} \right| + \xi _{{\rm{ESS}}}^{\rm{h}} \left| {P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{h}}} \right|} \right)}, \\ &C_{{\rm{VE}}{{\rm{P}}_i}}^{\rm{c}} = \sum\limits_{t = 1}^T {\left( {\xi _{{\rm{LT}}}^{\rm{c}} \left| {P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{c}}} \right| + \xi _{{\rm{ESS}}}^{\rm{c}} \left| {P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{c}}} \right|} \right)}. \\ \end{split}} \right\} $

3）目标函数3. 社会效益最高，受文献[15]中以用电舒适度表征虚拟电厂社会效益方式的启发，以用能舒适度来表征VEP的社会效益，即负荷切出率和转移率较低，社会效益较好.

(17) $ {\begin{split}{\rm{min}}\;{f_{2.3}} =& 3^{-1}\times\left\{ { \left[ { 1 - \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{e}}} \left( { \frac{{\lambda _{{\rm{LSI}}}^{\rm{e}}}}{{\lambda _{{\rm{LSI}}}^{\rm{e}} + \lambda _{{\rm{LSII}}}^{\rm{e}}}}\sum\limits_{t = 1}^T {\left( {\frac{{P{{_{{\rm{LSI}}}^{\rm{e}}}_{_{i,t}}}}}{{{\rm{Lim}}_{i.{\rm{max}}}^{{\rm{LSI}}}}}} \right)}+ } \right.} \right. } \right.\\ &\left. {\left. {\frac{{\lambda _{{\rm{LSII}}}^{\rm{e}}}}{{\lambda _{{\rm{LSI}}}^{\rm{e}} + \lambda _{{\rm{LSII}}}^{\rm{e}}}} \sum\limits_{t = 1}^T {\left( {\frac{{P_{{\rm{LSI}}{{\rm{I}}_{i,t}}}^{\rm{e}}}}{{{\rm{Lim}}_{i.{\rm{max}}}^{{\rm{LSII }}}}}} \right)} } \right)} \right] \times \\ &\left. {{\left[ {1 - \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{h}}} {\sum\limits_{t = 1}^T {\left( {\frac{{\left| {P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{h}}} \right|}}{{P_{{\rm{L}}{{\rm{T}}_i}}^{{\rm{max - h}}}}}} \right)} } } \right] + {1 - \sum\limits_{i = 1}^{N_{\rm{v}}^{\rm{c}}} {\sum\limits_{t = 1}^T {\left( {\frac{{\left| {P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{c}}} \right|}}{{P_{{\rm{L}}{{\rm{T}}_i}}^{{\rm{max - c}}}}}} \right)} } } }} \right\}. \end{split}} $

式中： $\lambda _{{\rm{LSI}}}^{\rm{e}} $、 $\lambda _{{\rm{LSII}}}^{\rm{e}} $、 $\lambda _{{\rm{LT}}}^{\rm{h}} $、 $\lambda _{{\rm{LT}}}^{\rm{c}} $分别为各类型负荷占该类总负荷比， ${\rm{Lim}}_{i,\max }^{{\rm{LSI}}} $、 ${\rm{Lim}}_{i,\max }^{{\rm{LSII}}} $、 ${P}_{{{\rm LT}}_i }^{{\rm{max-h}}}$、 ${P}_{{{\rm LT}}_i }^{{\rm{max-c}}}$分别为各类型可控负荷的总量. 从式（17）可以看出，min f_2,3的取值为[0, 1.0]，当各可控负荷在调度周期内完全不调用时， $P_{{\rm{LS}}{{\rm{{\rm I}}}_{i,t}}}^{\rm{e}}$、 $P_{{\rm{LSI}}{{\rm{I}}_{i,t}}}^{\rm{e}}$、 $P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{h}}$、 $P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{c}}$均为0，此时用电舒适度最高，min f_2,3=1；当各可控负荷调度总量达到上限时，用电舒适度最低，min f_2,3=0.

下层优化模型须满足各储能系统的相关约束. 其中储电约束如下.

(18) ${{{\rm{SOC}}}_{\min }}E_{{\rm{ES}}{{\rm{S}}_R}}^{\rm{e}} \leqslant E_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{e}} \leqslant {{{\rm{SOC}}}_{\max }}E_{{\rm{ES}}{{\rm{S}}_R}}^{\rm{e}},\;\forall t.$

(19) $ - P_{{\rm{ES}}{{\rm{S}}_t}}^{{\rm{e - ch}}} \leqslant P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{e}} \leqslant P_{{\rm{ES}}{{\rm{S}}_t}}^{{\rm{e - diss}}},\;\forall t.$

(20) $\sum\limits_{t = 1}^T {P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{e}}} = 0.$

(21) $\left. {\begin{gathered} E_{{{\rm{ESS}}_{i,t}}}^{\rm{e}} = \left( {1 - \rho _{{\rm{ESS}}}^{\rm{e}}} \right)E_{{{\rm{ESS}}_{i,\left( {t - 1} \right)}}}^{\rm{e}} - \Delta E_{{{\rm{ESS}}_{i,t}}}^{\rm{e}} ; \\ \Delta E_{{{\rm{ESS}}_{i,t}}}^{\rm{e}} = \left\{ \begin{gathered} P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{e}} \eta _{\rm{c}}^{\rm{e}},\;P_{{\rm{ES}}{{\rm{S}}_{_{i,t}}}}^{\rm{e}} > 0 ; \\ P_{{\rm{ES}}{{\rm{S}}_{_{i,t}}}}^{\rm{e}}/\eta _{\rm{d}}^{\rm{e}},\;P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{e}} \leqslant 0 . \\ \end{gathered} \right. \\ \end{gathered} } \right\}$

式中：SOC_min、SOC_max分别表示最小、最大充电状态， $P_{{\rm{ES}}{{\rm{S}}_t}}^{{\rm{e - ch}}} $、 $P_{{\rm{ES}}{{\rm{S}}_t}}^{{\rm{e - diss}}} $分别为最大充、放电功率， $\rho _{{\rm{ESS}}}^{\rm{e}} $、 $\eta _{\rm{c}}^{\rm{e}} $、 $\eta _{\rm{d}}^{\rm{e}} $分别为自放电率、充电率、放电率.

储冷储热系统运行方式相同，储热为例，约束如下：

(22) $\sum\limits_{t = 1}^T {P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{h}}} = 0.$

(23) $ - P_{{\rm{ES}}{{\rm{S}}_t}}^{{\rm{h - ch}}} \leqslant P_{{\rm{ES}}{{\rm{S}}_{i,t}}}^{\rm{h}} \leqslant P_{{\rm{ES}}{{\rm{S}}_t}}^{{\rm{h - diss}}},\;\forall t.$

式中： $P_{{\rm{ES}}{{\rm{S}}_t}}^{{\rm{h}} - {\rm{ch}}} $、 $P_{{\rm{ES}}{{\rm{S}}_t}}^{{\rm{h}} - {\rm{diss}}} $为最大充放电量.

LSI、LSII运行方式类似，以LSI为例，运行约束如下：

(24) $P_{{\rm{LS}}{{\rm{I}}_{\min ,t}}}^{\rm{e}} \leqslant P_{{\rm{LS}}{{\rm{{\rm I}}}_{i,t}}}^{\rm{e}} \leqslant P_{{\rm{LS}}{{\rm{I}}_{\max ,t}}}^{\rm{e}},\;\forall t.$

可转移的冷热负荷运行类似，以热负荷为例：

(25) $\sum\limits_{t = 1}^T {P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{h}}} = 0.$

(26) $P_{{{{\rm{LT}} }_{\min ,t}}}^{\rm{h}} \leqslant P_{{\rm{L}}{{\rm{T}}_{i,t}}}^{\rm{h}} \leqslant P_{{{{\rm{LT}} }_{\max ,t}}}^{\rm{h}},\;\forall t.$

对1.3.1节上层模型的目标函数2进行建模. 鉴于发电、用电的不确定性，RIES收益具有风险特征.

CVaR度量损失的平均情况可以描述尾部风险^[6]，CVaR为

(27) $\left. {\begin{split} & {{V_{{\rm{CVaR}},\;\beta }}\left( x \right) = E\left[ {{f^{{\rm{c {\text{-}} l}}}}\left( {x,y} \right)|{f^{{\rm{c {\text{-}} l}}}}\left( {x,y} \right) \geqslant {V_{{\rm{VaR}},\beta }}} \right]} =\\& {\qquad\qquad\quad \frac{1}{{1 - \beta }}{\smallint _{{f^{{\rm{c {\text{-}} l}}}}\left( {x,y} \right) \geqslant {V_{{\rm{VaR}},\beta }}}}{f^{{\rm{c {\text{-}} l}}}}\left( {x,y} \right) \times {f^{{\rm{PDF}}}}\left( y \right){\rm{d}}y},\\& {{V_{{\rm{VaR}},\beta }}\left( x \right) = {\rm{min}}\;\left\{ {{C_{\text{α}} } \in {\bf{R}}:P\left\{ {{f^{{\rm{c {\text{-}} l}}}}\left( {x,y} \right) \leqslant {C_\alpha }} \right\} \geqslant \beta } \right\}}=\\& {\qquad\qquad\quad {\rm{min}}\;\left\{ {{C_{\text{α}} } \in {\bf{R}}:\phi \left( {x,{C_\alpha }} \right) \geqslant \beta } \right\}}. \end{split}} \right\}$

式中：E(.)为期望函数；x∈Ω^D为决策变量；y∈Ω^R为随机变量，概率密度函数为f^PDF(y)，f^c-l(x,y)为RIES的收益损失函数，且E(|f^c-l(x,y)|)<+∞；C_α为损失值的阈值；VaR为在给定置信度β下，RIES可能遭受的最大损失值. 引入辅助函数计算CVaR，表示如下：

(28) $\left. {\begin{gathered} {V_{{\rm{CVaR}},\beta }}\left( x \right) = \mathop {\min }\limits_{{C_{\text{α}} } \in {\bf{R}}} {F_\beta }\left( {x,{C_{\text{α}} }} \right), \\ {F_\beta }\left( {x,{C_{\text{α}} }} \right)\! =\! {C_{\text{α}} } \!+\! \frac{1}{{1 \!+\! \beta }}{\int_{y \in {\varOmega ^{\rm{R}}}} {\left[ {{f^{{\rm{c {\text{-}} l}}}}\left( {x,y} \right) \!-\! {C_{\text{α}}}} \right]} ^ + }{f^{{\rm{PDF}}}}\left( y \right){\rm{d}}y. \\ \end{gathered} } \right\}$

式中：[t]⁺=max {t, 0}.

在RIES中，考虑发电与用电的双重不确定性，根据Copula函数的性质^[16]，根据单个随机变量的概率密度函数，可得多个随机变量耦合关系下的联合概率密度函数. 建立2种随机变量情况下的Copula-CVaR模型^[17]：

(29) $ {\rm{Copula - }}{V_{{\rm{CVaR}},\beta }} = \mathop {\min }\limits_{{C_{\text{α}} } \in {\bf{R}}} \left\{ {{C_{\text{α}} } + \frac{1}{{1 + \beta }}\int_{{y^{\rm{I}}} \in {\varOmega ^{\rm{RI}}}} {\left[ {\int_{{y^{{\rm{II}}}} \in {\varOmega ^{\rm{RII}}}} {\left( \begin{gathered} {\left[ {{f^{{\rm{c - l}}}}\left( {x,y} \right) - {C_\alpha }} \right]^ + } \times {f_1}\left( {{y^{\rm{I}}}} \right){f_2}\left( {{y^{{\rm{II}}}}} \right) \times C\left( {{F_1}\left( {{y^{\rm{I}}}} \right),{F_2}\left( {{y^{{\rm{II}}}}} \right)} \right) \end{gathered} \right){\rm{d}}{y^{\rm{I}}}} } \right]} {\rm{d}}{y^{{\rm{II}}}}} \right\} .$

式中：F₁(y^I)、F₂(y^II)、f₁(y^I)、f₂(y^II)分别为随机变量y^I、y^II的累计概率密度函数与概率密度函数；Ω^RI、Ω^RII由鲁棒优化理论进行构建，分别为描述风电、常规电负荷随机性的不确定合集.

以风电出力为例，利用鲁棒理论，对各时段的输出功率构建加法不确定合集：

(30) $\left. { \begin{split} & {\varOmega ^{{\rm{RW}}}} = \left\{ {P_{{{\rm{W}}_{i,t}}}^{\rm{e}} = P_{{{\rm{W}}_{i,t}}}^{{\rm{e - s}}} + {\gamma _{{{\rm{W}}_{i,t}}}} \Delta P_{{{\rm{W}}_{i,\;t}}}^{{\rm{e - u}}},\;\forall i,\;t;} \right.\\& \;\;\;\;\;\;\;\;\;\;\;\;\; \left. {{{\left\| {{{{\gamma}} _{{{\rm{W}}_{i,\;t}}}}} \right\|}_\infty } \leqslant 1,{{\left\| {{{{\gamma}} _{{{\rm{W}}_{i,t}}}}} \right\|}_1} \leqslant {{{\varGamma}} _{{\rm{W}},\;t}}} \right\};\\& {\left\| {{{{\gamma}} _{{{\rm{W}}_{i,t}}}}} \right\|_\infty } = \mathop {\max }\limits_{1 \leqslant i \leqslant {N_{\rm{W}}}} \left| {{\gamma _{{{\rm{W}}_{i,\;t}}}}} \right|,\;{\left\| {{{{\gamma}} _{{{\rm{W}}_{i,t}}}}} \right\|_1} = \sum\limits_{i = 1}^{{N_{\rm{W}}}} {\left| {{\gamma _{{{\rm{W}}_{i,\;t}}}}} \right|} . \end{split}} \right\} $

式中： $P_{{{\rm{W}}_{i,t}}}^{{\rm{e - s}}} $、Δ $P_{{{\rm{W}}_{i,t}}}^{{\rm{e }}} $分别为风电场i在t时段的预测出力与出力偏差；Δ $P_{{{\rm{W}}_{i,t}}}^{{\rm{e - u}}} $为出力偏差的上限； ${\gamma _{{{\rm{W}}_{i,t}}}} $为出力偏差系数；||·||_∞为无穷范数；||·||₁≤ ${\varGamma}_{{\rm{W}},t} $表示1范数约束对应不确定变量的空间集群效应，既在某个调度时段各风电场的出力偏差不可能同时达到最大，由此引入空间约束参数 ${\varGamma}_{{\rm{W}},t}$来调整不确定合集的边界. 若 ${\delta _{{{\rm{W}}_{i,t}}}} $=| ${\gamma _{{{\rm{W}}_{i,t}}}} $|，Δ $P_{{{\rm{W}}_{i,t}}}^{\rm{e}} $独立且服从正态分布，记期望和方差为0和σ_W^*，利用Lindeberg-Levy中心极限定理，可得

(31) $ {\varGamma }_{\rm{W}}\!=\!{N}_{\rm{W}} \sqrt{\frac{2{\sigma }_{\rm{W}}{}^{{*}^{2}}}{\text{π}}}+{\varPhi }^{-1}\left({\alpha }_{\rm{W}}\right)\sqrt{{N}_{\rm{W}} \left({\sigma }_{\rm{W}}{}^{*}-\frac{2{\sigma }_{\rm{W}}{}^{{*}^{2}}}{\text{π} }\right)};\;\forall t.$

式中：Φ⁻¹（·）为正态分布密度函数的反函数，α_W为风电置信概率.

通过构造拉格朗日函数与线性对偶理论可知，考虑在t时段的最极端情况，风电场出力达到不确定合集下限. 此时仅有一个风电场出力的偏差系数不足1，设该风电场为j，如下：

(32) $\begin{split} \sum\limits_{i = 1}^{{N_{\rm{W}}}} {P_{{{\rm{W}}_{i,t}}}^{\rm{e}}} = & \sum\limits_{i = 1}^{{N_{\rm{W}}} - 1} {\left( {P_{{{\rm{W}}_{i,t}}}^{{\rm{e - s}}} - \Delta P_{{{\rm{W}}_{i,t}}}^{{\rm{e - u}}}} \right) + } P_{{{\rm{W}}_{i,t}}}^{{\rm{e - s}}} - \\& \left( {{\varGamma _{{\rm{W}},t}} -\left\lfloor { {{\varGamma _{{\rm{W}},t}}} } \right\rfloor} \right)\Delta P_{{{\rm{W}}_{j,t}}}^{{\rm{e - u}}}. \end{split}$

同理建立用电不确定性合集，可得空间约束参数 ${\varGamma _{_{{\rm{L}},t}}^{\rm{e}}}$与极端功率情况，如下：

(33) $\begin{split} \sum\limits_{i = 1}^{N_{\rm{L}}^{\rm{e}}} {P_{{{\rm{L}}_{i,t}}}^{\rm{e}}} = & \sum\limits_{i = 1}^{N_{\rm{L}}^{\rm{e}} - 1} {\left( {P_{{{\rm{L}}_{i,t}}}^{\rm{e}} + \Delta P_{{{\rm{L}}_{i,t}}}^{{\rm{e - u}}}} \right)} + P_{{{\rm{L}}_{j,t}}}^{{\rm{e - s}}} - \\ & \left( {\varGamma _{_{{\rm{L}},t}}^{\rm{e}} - \left\lfloor { {\varGamma _{_{{\rm{L}},t}}^{\rm{e}}} } \right\rfloor } \right)\Delta P_{{{\rm{L}}_{j,t}}}^{{\rm{e - u}}}. \end{split}$

为了定量分析RIES的收益风险，建立基于Copula-RCVaR的多能流收益风险模型. 模型中，x为上层目标的决策变量，随机变量y^I为风电出力偏差Δ $P_{{{\rm{W}}_{i,t}}}^{\rm{e}} $，y^II为常规电负荷偏差Δ $P_{{{\rm{L}}_{i,t}}}^{\rm{e}} $，定义系统运行时的损失函数为利润函数的负数，f^c-l(x,y)=−f_1.1，上层模型的目标函数2（min f_1.2）为式（29）的形式.

鲸鱼算法（WOA）具有参数设置少、寻优性能强等特点，在求解精度和收敛速度上均优于粒子群算法PSO^[18]，已成功应用于大规模优化问题上. 标准WOA存在不能有效平衡全局与局部搜索能力，导致在迭代后期算法的多样性丧失，收敛能力不足的问题，如在文献[18]测试问题F2和F21上，算法在迭代最终收敛. 提出相关的改进策略，改进的鲸鱼算法（improved WOA，IWOA）伪代码如下.

算法：IWOA

输入：Np（种群规模）；D（维度）；G（最大迭代次数）；A_constant; X（初始种群）

输出：x^*（最优个体）

1.F←计算X适应度；x^*←从X中选择最优个体；

2.while （迭代停止条件不满足） do

3. 通过式（34）、（35）更新a, A, C, l；

4. if |A|≥A_constant

5. 　　在X中随机选择不同5个个体（x_r1, x_r2, x_r3, x_r4, x_r5）；

6.　　通过式（36）更新X；

7. else 在X中随机选择不同2个个体（x_r1, x_r2）；

8. 通过式（37）更新X；end if

9. 越界处理；计算F；更新x^*；end while

10. return x^*

(34) $\left. {\begin{split} & a = 2\left[ {1 - {{\left( {{{{G_{{\rm{iter}}}}}}\left/{{{{G_{\max }}}} }\right.} \right)}^{ {{l}/{{1.5}}} }}} \right], \\& A = 2 a r - a. \end{split}} \right\} $

(35) $\left. {\begin{split} & C = 2 r ,\\ & l = 1 - \left( {\frac{{G_{\rm{i ter}}{\rm{ + }}{G_{\max }}}}{{{G_{\max }}}} - 1} \right){{r}}. \end{split}} \right\} $

式中：G_iter、G_max分别为当前迭代次数与最大迭代次数；r为（0，1.0）的随机数；系数A、C均由收敛因子a计算，随着迭代次数由2减小到0；l为螺旋系数. 设置探索固定值A_constant，当A≥ A_constant时执行全局搜索，反之为局部. 借助差分进化算法中个体的合作与竞争指导优化搜索，分别进行螺旋运动和直线运动，更新方式如下：

(36) $\left. {\begin{split} & {{x}}(t + 1) = \left\{ \begin{aligned} & {{{x}}_{{\rm{r}}1}}\left( t \right) - A {{{D}}_1} ,\;\;\;\;A \geqslant {\rm{A\_constan}} {\rm{t}};\\& {{{D}}_1} {{\rm{exp}}\;(bl)} \cos \left( {2{\text{π}} l} \right) + {{{x}}_{{\rm{r}}1}}\left( t \right) ,\;A\! <\! {\rm{A \_ constant}}; \end{aligned} \right. \\& {{{D}}_1} = {C\left( {{{{x}}_{{\rm{r}}2}}\left( t \right) + {{{x}}_{{\rm{r}}3}}\left( t \right) - {{{x}}_{{\rm{r}}4}}\left( t \right) - {{{x}}_{{\rm{r}}5}}\left( t \right)} \right)} . \end{split} } \right\}$

(37) $\left. {\begin{split} & {{x}}(t + 1) = \left\{ \begin{aligned} & {{{x}}^{\rm{*}}}\left( t \right) - A {{{D}}_2},\;\;\;\;A \geqslant {\rm{A\_constan}} {\rm{t}} ; \\ & {{{D}}_2} {{\rm{exp}}\;(bl)} \cos\; \left( {2{\text{π}} l} \right) \!+\! {{{x}}^{\rm{*}}}\left( t \right)\! ,\;A\! < \!{\rm{A \_ constan}} {\rm{t}}; \end{aligned} \right. \\ & {{{D}}_2} = {C \left[ { {{{{x}}^{\rm{*}}}\left( t \right){\rm{ - }}{{x}}\left( t \right)} + {{{{x}}_{{\rm{r}}1}}\left( t \right) - } {{{x}}_{{\rm{r}}2}}\left( t \right)} \right]} . \end{split} } \right\}$

多目标鲸鱼算法借鉴NSGAII中的精英保留策略，利用外部存档保存进化过程中已经发现的非占优解. 当外部存档超出设定的最大容量时，采用拥挤熵的方式对Pareto解集进行裁剪^[19]. 该方法考虑相邻解的分布情况，能够合理反映非支配解之间的拥挤程度. 从问题的实际出发，需要得到一个满足各个目标的解，使用模糊数学的方式提取最优折中解，选择线性函数作为隶属度函数.

针对RIES两阶段风险能量管理模型中复杂的等式与不等式约束，采用滤子技术对约束条件进行处理. 构造由目标函数与约束违反度组成的数对（F, G）来表示滤子^[20]，

(38) $G = \sum\limits_{i \in m} {\max \left( {0,{g_i}({{Y}})} \right) + \sum\limits_{j \in n} {\left| {{h_n}({{Y}})} \right|} } .$

式中：g_i(Y)、h_n(Y)分别为不等式与等式约束，m、n为对应的个数. 借助Pareto理论，在最小值问题上有如下定义.

定义1 　若F(Y_i)≤F(Y_j)，G(Y_i)≤G(Y_j)，则称滤子(F(Y_i))，G(Y_i))支配(F(Y_j)，G(Y_j)).

定义2 　滤子集内的滤子互不支配.

将上层模型中各集线器能量约束（式（12）~（14））与目标函数构造滤子对；其他约束均可以作为边界条件，直接利用元启发式算法处理. 下层模型储电侧电荷约束（18）与目标函数构造滤子对；可以转移冷热负荷、储能系统的可持续运行约束（20）、（22）、（25），采用动态可松弛约束处理方式^[21]. 以储电为例，计算约束违反程度记为ε_ESS-e，根据边界条件计算松弛度，根据松弛度确定调整量；其他约束可以作为边界条件.

模型整体求解包括约束处理流程，如图2所示.

以修改的IEEE33节点配电网与CCHP系统耦合形成RIES，CCHP内设备及参数见表1、2. 表中，P_max、P_min分别为功率的上、下界，η为能效，C_c为成本价格，GT为燃气轮机，WHR为余热回收装置，ABS为吸收式制冷机，ASR为电制冷机，BO为锅炉。RIES包含3个虚拟能量厂，分别实现RIES内电、热、冷负荷的需求响应，虚拟能量厂的相关参数如表3所示. 表中，VEP-e、VEP-h、VEP-c分别表示电、热、冷的虚拟能量厂，P_t为占比，SOC为容量，P_ESS为最大充放电功率，SOC_pu为归一化后的容量。配电网中1为根节点，与上级电网相连，节点15接入总容量为55 MW的风电厂群. 区域内电、冷、热负荷及风电出力预测见图3. 图中，P_L为预测电荷。电负荷的85%购自日前市场，设购买价格为0.4 kW·h/美元，电能在实时市场的交易价格与用电量有关，如图4所示，天然气购买价格为0.22 Kcf·h/美元. RIES的电、热、冷售价见图4. 图中，C_c为价格。风电、常规电负荷的惩罚系数为0.65、0.60 kW/美元.

鲁棒优化是在不确定变量的极端情况下系统进行的优化调度. 根据2.3节的分析，可以推出系统在所考虑的极端情况之外运行的概率：

(27) $\begin{split} & P_{{\rm{RIES}}}^{{\rm{out}}} \leqslant \\& {\exp\; \left( { - \frac{{{\varGamma _{{\rm{W}},t}}^2}}{{2{N_{\rm{W}}}}}} \right) + \exp\; \left( { - \frac{{\varGamma {{_{_{{\rm{L}},t}}^{\rm{e}}}^2}}}{{2N_{\rm{L}}^{\rm{e}}}}} \right) + \exp\; \left( { - \frac{{{\varGamma _{{\rm{W}},t}}^2}}{{2{N_{\rm{W}}}}} - \frac{{\varGamma {{_{_{{\rm{L}},t}}^{\rm{e}}}^2}}}{{2N_{\rm{L}}^{\rm{e}}}}} \right)} .\end{split} $

为了分析风电与负荷的置信概率与总数和系统运行在所考虑极端情况外的概率P_OE的关系，分别针对单个不确定变量与多不确定变量互相耦合的情况进行研究. 图5中，α_W、α_L分别为风电置信概率和常规电负荷置信概率，N_W、N_L分别为风电场数量和常规电负荷总数。如图5（a）所示为单个不确定变量（以风电为例），如图5（b）、（c）所示分别为2个不确定变量耦合. 图5（a）中，常规电负荷总数、置信概率固定分别为20、0.6；图5（b）中，风电常规电负荷总数均为15；图5（c）中，风电常规电负荷置信概率均为15.

从图5（a）可以看出，不确定变量的置信概率增大，超出极端情况的概率降低；不确定变量总个数减小，该概率升高. 多个不确定变量耦合下超出极端情况概率的等高线间距增加且不等，说明该情况下不确定变量对系统的影响更加复杂. 较图5（c）、（b）中小概率等高线包含区域较小，置信概率对系统超出极端情况的概率影响较明显.

基于建立的Copula-RCVaR模型，对以下4个算例进行分析. 算例1：RIES含虚拟冷/热/电厂；算例2：RIES仅含虚拟热厂与虚拟冷厂；算例3：RIES仅含虚拟电厂；算例4：RIES完全不含虚拟能量厂，可调度仅为燃气轮机与锅炉. 设发电与用电偏差服从正态分布（预测精度为68.27%），考虑空间集群效应，总数量均为20，置信概率均为0.6. 运行结果见表4，虚拟能量厂优化方案见图6. 表4中，B_RIES为RIES利润，D_VEP-e、D_VEP-h、D_VEP-c分别为VEP-e、VEP-h、VEP-c调度偏差功率，C_VEP-e、C_VEP-h、C_VEP-c分别为VEP-e、VEP-h、VEP-c调度成本，S_VEP-e、S_VEP-h、S_VEP-c分别为归一化后的VEP-e、VEP-h、VEP-c社会成本。图6中，P为功率，LSI-e、LSII-e分别表示迎峰电负荷和避峰电负荷，ESS-e表示储电系统，LT-h表示可转移热负荷，ESS-h表示储热系统，LT-c表示可转移冷负荷，ESS-c表示储冷系统。

从表4可以看出，随着不同类型VEP的加入，对更多种类的可调度资源与储能装置集中管理，系统内包含的可调度资源种类增加，调度变得更加灵活，偏差随之减小，情况1（虚拟冷、热、电厂全参与的情况）下的偏差较情况3（仅有虚拟电厂参与的情况）下减小6%. 虚拟冷/热厂中包含的可控负荷主要为TL-h/TL-c，基于该类型负荷转移前后负荷总量不变的强约束条件，使得可转移负荷数量增加，RIES与VEP之间的偏差大大降低. 随着可调控资源数量的增加，分摊了VEP在调控时的经济与社会成本，各类可控资源充分全面参与调度，VEP的经济运行成本在VEP全参与下较仅有虚拟电厂时减少36.3%，较虚拟冷、热厂参与时分别减少17.1%、6%；社会成本相应提高，用户的用电舒适度增高；RIES利润逐渐增加，VEP全参与下的经济成本较不含VEP降低1.9%.

从图6可以看出，在用电高峰时段（11时—13时、19时—22时），VEP向RIES注入能量，保证供需与电量平衡；RIES将盈余电量以较高的实时电价，在能量交易中心通过实时市场较平稳外送至上级电网，在保证大电网稳定运行的前提下解决负荷集中地区的高峰用电需求. 在低耗电时期，VEP向RIES吸收能量，以满足自身区域内可控资源的运行需求. 对比图6中各算例VEP的调度方案可知，VEP全参与下的计划较其他2种方式更平稳，图6（a）的累积调度相对集中在[−600, 600] kW，与表4的结果吻合.

分析各不确定合集的置信概率对RIES利润、收益损失风险的影响. 在相同风险阈值下，CVaR置信度为0.95；不确定变量的预测精度均为68.27%；风电场、常规电负荷总数分别为20、32，不同置信概率α下的RIES利润、收益损失风险见表5. 表中，收益损失风险为归一化后的数值. 可以看出，随着置信概率的不断减小，不确定合集区间逐渐收缩，系统所需旋转备用成本不断减小，RIES利润随之升高；系统的收益损失风险为仅考虑系统不确定合集内的不确定性计算而来，由于不确定合集收缩，RIES收益损失风险逐渐减小，意味着系统运行时面临的风险逐渐减小. 盲目减小置信概率，会使得系统运行在极端情况外的概率大大增加，当置信概率降至20%时，该概率为100%. 当置信概率为30%~45%时，RIES利润增加最快，收益损失风险下降最快；当置信概率为45%~60%时，极端情况外运行概率处于可接受的低概率段.

如图7所示为当α=55.5%时，上层与下层的Pareto有效前沿. 图中，CVaR_RIES为RIES收益风险，D_B为RIES利润的相反数。可以看出，利用改进的多目标鲸鱼算法得到的Pareto解集较均匀地分布在Pareto前沿上，具有较好的分布性. Pareto解集中的每一点对应在该利润与收益损失风险下的RIES及各VEP的优化运行策略. 系统调度员可以根据实际中的不同情况，平衡RIES风险与利润、各VEP的偏差与成本进行决策，寻找合适的最优折中解.

为了说明Pareto最优解集为有效解，当α=55.5%时双层模型中各目标函数的收敛情况如图8所示. 图中，D_VEP为调度偏差功率，N_it为迭代次数。由于上层目标的收益损失风险与下层目标的社会效益度量尺度较小，将社会效益乘以10 000后，分别与调度计划偏差、各VEP经济偏差相加. 从图8可以看出，Pareto解集是在双层模型各目标函数均收敛下得到的有效解.

分析不确定变量的空间集群效应对RIES利润、收益损失风险的影响. 将接入的风电场与常规电负荷在总功率不变的情况下细分为5种. 1）情况1：风电场/常规电负荷总数均为2；2）情况2：风电场/常规电负荷总数均为10；3）情况3：风电场、常规电负荷总数为20、10；4）情况4：风电场、常规电负荷总数均为20；5）情况5：风电场、常规电负荷总数均为30.

分析在上述5种不同情况时，不确定集合的置信概率与系统在极端情况外运行的概率关系，结果见图9. 各情形在不同不确定集合置信概率下的RIES利润与收益损失风险值见图10.

随着不确定变量总数的增加，分布空间更广阔，空间约束参数相应增大，各不确定变量的波动性被更加细致地描述，表现为系统旋转备用成本、VEP调度成本降低，RIES利润增加. 由图9可知，系统在情况5下，超出极端情况的概率随着置信概率的变化更多处在数值较小区域。在情况1下，该概率更多地处在较大区域，且不管置信概率如何增加，也不会出现在低风险区域. 从图10可以看出，随着不确定变量总数的增加，在相同置信概率下的系统收益损失风险从情况1到情况5呈递增趋势.

（1） 该系统能够在保证系统安全稳定运行的前提下，平衡各层级的主要收入，有效地实现电力系统的可持续发展.

（2） 引入空间约束参数，弥补了传统鲁棒优化过于保守的不足，降低不确定性决策的盲从性.

（3） 分析不确定集合置信度概率、空间集群效应对系统利润、收益损失风险和极端场景外运行概率的影响. 在考虑系统保守性和提高经济性的同时，对合理选择敏感性参数有理论指导意义.

（4） 利用提出算法有效求解了包含复杂约束的混合整数非凸非线性规划问题.