【鲁棒优化、机会约束】具有分布鲁棒联合机会约束的能源和储备调度研究（Matlab代码实现）

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码、数据、文章讲解

💥1 概述

文献来源：

本文具有分布鲁棒联合机会约束的能源和储备调度研究，用于可再生能源渗透率高的联合电力和天然气系统的能源和储备调度。数据驱动的分布稳健的机会约束确保没有减载和高概率的可再生溢出。我们使用条件风险值近似和线性决策规则有效地解决了这个问题。样本外实验表明，该模型在没有机会约束的情况下主导了相应的随机程序，明确模拟了减载和可再生溢出的影响。

原文摘要：

Abstract

We develop a two-stage stochastic program for energy and reserve dispatch of a joint power and gas system with a high penetration of renewables. Data-driven distributionally robust chance constraints ensure that there is no load shedding and renewable spillage with high probability. We solve this problem efficiently using conditional value-at-risk approximations and linear decision rules. Out-of-sample experiments show that this model dominates the corresponding stochastic program without chance constraints that models the effects of load shedding and renewable spillage explicitly.

风能、太阳能和潮汐能等可再生能源部署的增加极大地改变了发电结构。虽然有利于可持续性，但可再生能源由于其间歇性和有限的可预测性而损害了输电系统的稳定性[27]。因此，灵活的燃气发电厂（GFPP）被广泛建造，以取代退役的火力发电厂或核电站，以平衡不稳定的可再生能源发电。因此，可以预见电力和天然气系统的更紧密耦合[28]。这种观点促使我们同时研究这两个系统。

在存在不确定的可再生能源的情况下，如果实际可再生能源产量偏离其预测，则可能无法实施提前调度。在这种情况下，常规发电厂需要实时调整其生产水平。如果在某些极端情况下，这些工厂不够灵活，无法恢复输电系统的完整性，那么可再生能源溢出或减载可能是必要的。然而，这种严厉的措施会产生高昂的经济成本。

能源和储备调度问题可以通过鲁棒优化[3]，随机规划[26]和机会约束规划[24]的方法来解决。稳健的优化模型在规定的不确定性集[4]，[35]内实现可再生预测误差的最坏情况下，最大限度地减少了日前调度和纠正追索措施的成本。然而，对最坏情况的关注以及在线性决策规则中解决较大问题实例的必要性可能会导致过于保守的解决方案。另一方面，随机规划模型在规定的离散分布[21]，[23]下最小化预期成本，但其解决方案可能会显示出较差的样本外性能，除非离散化点的数量随问题维度呈指数增长。最后，机会约束的编程模型不考虑可再生溢出和减载，但确保系统在没有这些严厉的追索措施的情况下保持稳定，概率很高[5]，[18]，[25]。当难以甚至不可能以有意义的方式将追索权和不可行性货币化时，机会受限的编程公式特别有吸引力。不利的一面是，这些配方通常不能保证产生具有可行追索权的可实施的第一阶段解决方案。

在本文中，我们将能量和储备调度问题表述为一个数据驱动的分布鲁棒机会约束程序，而不是一个位于 Wasserstein 模糊集上的程序，即概率分布空间中的一个球，相对于 1 型 Wasserstein 度量，以给定训练数据集上的经验分布为中心 [20].在最小化最坏情况下的预期成本的同时，我们的模型防止了违反储备裕度，并通过分配稳健的联合机会约束来执行电力和天然气网络的传输容量限制。这意味着对模糊性集中的所有分布强制执行基础经典机会约束。分布鲁棒优化结合了随机规划的特殊性（通过概率分布对不确定性进行建模）和鲁棒优化的保守性（通过对冲分布不确定性）。所提出的方法承认可用的统计数据通常可以用许多不同的分布来解释，并减轻了基于单个分布的经典随机程序的过拟合效应特征，该分布总是被估计误差破坏。因此，如果相对于不确定问题参数的维数而言，训练样本很少，则基于 Wasserstein 模糊性集的分布稳健方法特别有用，这在大型电力系统中通常就是这种情况。

现有文献中大多数分布稳健的能量模型都基于矩模糊集[9]，其中包含共享相同平均向量和协方差矩阵或满足一组广义矩约束[1]，[2]，[29]，[31]，[32]，[33]的所有分布。分布稳健的个体机会约束在[2]，[33]中研究，而双侧机会约束在[31]中解决。真正的联合机会约束强制要求以高概率同时满足多个安全条件。它们比多个个人机会约束更具表现力，更不保守。例如，100条输电线路中任何一条超载的风险最多为1%（联合机会约束）的操作制度比每条输电线路最多1%的过载风险（100个单独的机会约束）的制度更可取。事实上，在后一种情况下，任何线路过载的风险可能高达100%（例如，如果每天只有一条传输线以随机顺序发生故障）。不幸的是，联合机会约束在计算上不如单个机会约束容易处理。通过使用邦弗朗尼不等式将联合机会约束分解为几个单独的机会约束来获得可处理的保守近似。例如，在我们的示例中，通过将每条传输线的单个过载概率限制在 1.0% 以下，可以将任何传输线的联合过载概率推低到 01% 以下。不幸的是，这种邦弗朗尼近似可能过于保守。通过将非凸关节机会约束替换为凸条件风险值（CVaR）约束，可以获得替代的可处理保守近似[22]。对于具有平均协方差模糊设置的分布稳健联合机会约束，最佳CVaR近似是精确的，但识别这种最佳近似仍然很困难[36]。精确可处理的重新表述仅适用于涉及一阶离散测度的限制性矩模糊集[14]。

分布稳健的联合机会约束程序在Wasserstein模糊集合上是NP难的，但如果决策变量和不确定参数是可分的，它们允许精确的混合整数圆锥重新表述[6]，[15]，[30]。[11]中提出了一个可处理的鲁棒近似。

本文的贡献总结如下：

(1)
我们为GFPP制定了一个两阶段的能源和储备调度随机计划，其中包含燃料约束。具有Wasserstein模糊集的分布稳健联合机会约束确保不需要高概率的负荷卸减和可再生泄漏。为了使这个问题易于处理，我们用线性决策规则来解决它，并为机会约束开发了一系列CVaR近似。我们还设计了一种顺序凸优化算法来调整CVaR近似。
(2)
我们开发了一个基于实时最佳潮流模型的模拟环境，具有减载和可再生溢出功能，使我们能够评估不同日前调度策略的样本外性能。
(3)
通过数值模拟表明，所提出的CVaR近似帕累托在样本外成本和机会约束的经验违反概率方面在邦弗朗尼近似中占主导地位。我们还表明，考虑到样本外成本的均值和可变性，所提出的模型比没有机会约束的两阶段随机程序更可取，后者明确地模拟减载和可再生溢出。

我们的研究结果表明，为了最小化样本外成本，求解分布鲁棒机会约束程序可能比解决（看似更现实的）两阶段分布鲁棒优化问题更好。这可能令人惊讶，因为机会约束程序忽略了所有低概率情景，在这些情景中，电力系统的稳定性只能通过昂贵的追索行动（如减载和可再生泄漏）来维持，而两阶段问题考虑了所有可能的情况，并如实模拟了必要的追索行动及其成本。我们认为，从两阶段问题中获得的操作策略的样本外性能较差，源于有效解决两阶段问题所需的最新决策规则近似导致的最优性显着损失。因此，忽略触发昂贵追索权行动的低概率情景是有好处的，这些行动在两阶段模型中总是会产生较大的近似误差。我们论文的这一关键见解可以作为机会约束能量模型流行的可能解释，这些模型通常在没有严格理由的情况下临时使用。

📚2 运行结果

部分代码：

% Inverse of logit-normal transformation (Eq. (2) in ref. [31])

R = chol(sigma_m);
y = repmat(mu,Nscen,1) + randn(Nscen,size(WindDATA,1))*R;
Wind = (1+exp(-y)).^(-1);

% Checking correlation, mean and true mean of data

corrcoef(Wind);
mean(Wind);
true_mean_Wind = (1+exp(-mu)).^(-1);

% Reshaping the data structure

nWind = Wind';
nWind = reshape(nWind,size(WindDATA,1), N_max+OOS_max, IR_max);

% Initializing the matrices to gather final results

Joint_CVaR_Obj_IR = zeros(IR_sim, OOS_sim, length(rho_vectorC));
CVaR_Obj_IR = zeros(IR_sim, OOS_sim, length(rho_vectorJC));
ICC_TC = NaN(IR_sim,length(rho_vectorC));
JCC_TC = NaN(IR_sim,length(rho_vectorJC));

% Loop for each individual run for 100 coupled datasets

for j = 1:IR_sim
display('out of sample iteration:');
j

% For each coupled dataset, we pick N and N' samples
WPf_max = nWind(:,1:N_max,j)';
WPr_max = nWind(:,N_max+1:N_max+OOS_max,j)';
WPf = WPf_max(1:N,:);
WPr = WPr_max(1:OOS_sim,:);

% Build the corresponding data related to wind power production
all = [1:N];
system_info.Wscen = WPf(all,:)';
system_info.mu = mean(system_info.Wscen,2);
system_info.xi = system_info.Wscen - repmat(system_info.mu, 1, size(system_info.Wscen,2));

% Calculation of A,B,C,b matrices for joint chance constraints
CC_jcc = CC_matrices(system_info, DRO_param);
jcc = CC_jcc.jcc;

% Loop for each value of \rho in P vector
for i = 1:length(rho_vectorC)

% optimize for each value of rho for Bonferroni approximation

DRO_param.rho = rho_vectorC(i);

DRO_ICC_CVaR = DRO_CVaR_ICC(system_info, DRO_param, jcc);
ICC_p_DA{j, i} = DRO_ICC_CVaR.p;
ICC_ru{j, i} = DRO_ICC_CVaR.ru;
ICC_rd{j, i} = DRO_ICC_CVaR.rd;
ICC_obj{j, i} = DRO_ICC_CVaR.Obj;
ICC_flag{j, i} = DRO_ICC_CVaR.Flag;
CVaR_Y{j,i} = DRO_ICC_CVaR.Y * system_info.xi;
CVaR_Qy{j,i} = DRO_ICC_CVaR.q;
CVaR_QY{j,i} = DRO_ICC_CVaR.qY * system_info.xi;
CVaR_Fy{j,i} = DRO_ICC_CVaR.fy;
CVaR_FY{j,i} = DRO_ICC_CVaR.fY * system_info.xi;

% optimize for each value of rho for Zymler approximation

DRO_param.rho = rho_vectorJC(i);

DRO_JCC_CVaR = DRO_JCVaR_All(system_info, DRO_param, jcc);
JCC_p_DA{j, i} = DRO_JCC_CVaR.p;
JCC_ru{j, i} = DRO_JCC_CVaR.ru;
JCC_rd{j, i} = DRO_JCC_CVaR.rd;
JCC_obj{j, i} = DRO_JCC_CVaR.Obj;
JCC_flag{j, i} = DRO_JCC_CVaR.Flag;
Joint_CVaR_Y{j,i} = DRO_JCC_CVaR.Y * system_info.xi;
Joint_CVaR_Qy{j,i} = DRO_JCC_CVaR.q;
Joint_CVaR_QY{j,i} = DRO_JCC_CVaR.qY * system_info.xi;
Joint_CVaR_Fy{j,i} = DRO_JCC_CVaR.fy;
Joint_CVaR_FY{j,i} = DRO_JCC_CVaR.fY * system_info.xi;

% Loop for each out-of-sample realization
for k = 1:OOS_sim
system_info.Wreal = WPr(k,:)';
system_info.DWreal = system_info.Wreal - system_info.mu;

% Solve real-time optimal power flow for the solution of Bonferroni
% approximation
RT_solution_CVaR = RT_solve_R(system_info,DRO_ICC_CVaR.p,DRO_ICC_CVaR.ru,DRO_ICC_CVaR.rd);
CVaR_Obj_IR(j,k,i) = RT_solution_CVaR.Obj_RT;
CVaR_lshed{j,k,i} = RT_solution_CVaR.lshed_RT;
CVaR_flow{j,k,i} = DRO_ICC_CVaR.fy + DRO_ICC_CVaR.fY * system_info.DWreal;
CVaR_p{j,k,i} = DRO_ICC_CVaR.Y * system_info.DWreal;
CVaR_q{j,k,i} = DRO_ICC_CVaR.q + DRO_ICC_CVaR.qY * system_info.DWreal;
CVaR_flag(j,k,i) = RT_solution_CVaR.Flag;

% Solve real-time optimal power flow for the solution of Zymler
% approximation
RT_solution_Joint_CVaR = RT_solve_R(system_info,DRO_JCC_CVaR.p,DRO_JCC_CVaR.ru,DRO_JCC_CVaR.rd);
Joint_CVaR_Obj_IR(j,k,i) = RT_solution_Joint_CVaR.Obj_RT;
Joint_CVaR_lshed{j,k,i} = RT_solution_Joint_CVaR.lshed_RT;