💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
📋📋📋本文目录如下：🎁🎁🎁
目录
⛳️赠与读者
👨💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能解答你胸中升起的一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。
或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎
💥1 概述
本研究提出了一种新颖的基于矩的分布鲁棒优化（DRO）模型，该模型结合了条件风险价值（CVaR），用于应对电力价格不确定性下的自调度问题。该模型综合考虑了电力价格波动、机组参数以及负荷需求，并且可以通过调整模糊集的大小来进行调节，从而为发电厂提供了一种合适且可调节的自调度策略。决策者可以根据实际场景调整该策略，使其被独立系统运营商（ISOs）接受，同时最大化发电利润。通常，此类DRO模型会被转化为半定规划（SDP）来求解，然而，求解大规模SDP需要大量的计算时间和资源。针对这一不足，提出了两种有效的近似模型：一种是基于向量拆分的近似模型，另一种是基于交替方向乘子法（ADMM）算法的近似模型，两者都可以极大地减少计算时间和资源。为了验证我们所提出模型的正确性和有效性，我们将它们应用于多种数值案例研究，并与文献[1]中的模型进行了比较。对三个IEEE测试系统（6节点系统、30节点系统和118节点系统）进行了仿真。
文献1：
利用条件风险价值应对价格不确定性下的鲁棒自调度
摘要 ： 在电力行业解除管制的情况下，发电公司在日前和小时级电力市场中进行投标，以试图实现利润最大化。对于一个成功的投标策略，每个发电公司都需要生成基于最优自调度的投标曲线。这种自调度通常是基于预测的节点边际电价（LMPs）的利润最大化最优潮流模型获得的。然而，在自调度时，LMPs的预测值存在很大的不确定性。因此，希望能够产生鲁棒的自调度，以减轻因暴露于价格波动而产生的风险。在投资组合优化理论中，风险管理的方法包括风险价值（VaR）和条件风险价值（CVaR）。CVaR被认为比VaR是一个更一致的风险衡量指标。实际上，尽管CVaR是平均超额损失，但VaR对于超出该指标金额可能遭受的损失程度没有任何指示。本研究提出了一种基于CVaR的鲁棒自调度方法。将展示如何使用多项式内点法从二阶锥规划中获得鲁棒的自调度。所获得的调度方案在最大化利润和最小化风险之间提供了一种折中解决方案。将在标准IEEE母线测试系统上进行的仿真结果将用于展示基于CVaR的调度模型。
摘要 ：为确保投标成功并最大化利润，发电厂需要一种能够应对各种场景的自调度策略。因此，分布鲁棒优化（DRO）是一个不错的选择，因为它可以在不确定的环境中为GENCOs提供一种可调节的自调度策略，与通过鲁棒优化（RO）和随机规划（SO）得出的策略相比，它可以很好地平衡鲁棒性和经济性。在本文中，提出了一种新颖的基于矩的分布鲁棒优化（DRO）模型，该模型结合了条件风险价值（CVaR），用于解决电力价格不确定性下的自调度问题。通常，此类DRO模型会被转化为半定规划（SDP）来求解，然而，求解大规模SDP需要大量的计算时间和资源。针对这一不足，提出了两种有效的近似模型：一种是基于向量拆分的近似模型，另一种是基于交替方向乘子法（ADMM）的近似模型，两者都可以极大地减少计算时间和资源，并且第二种近似模型在求解的每一步中仅需要当前区域的信息，从而保证了信息的私密性。通过在三个IEEE测试系统上进行仿真，验证了所提出的DRO模型和两种近似模型的正确性和有效性。
关键词 ：自调度，分布鲁棒优化，价格不确定性，CVaR，ADMM
一、问题背景：电力自调度中的价格不确定性挑战
自调度（Self-Scheduling）指发电企业或微网运营商在电力市场中自主制定发电计划与投标策略，以最大化利润或最小化成本。其核心挑战在于：
- 电价不确定性：日前市场电价受负荷波动、新能源出力、网络阻塞等影响，呈现强随机性。
- 尾部风险：极端低价场景可能导致发电商巨额亏损，传统随机规划忽略尾部风险。
- 分布模糊性：电价真实分布未知，仅能通过历史数据估计矩信息（如一、二阶矩），但参数估计存在误差。
传统方法缺陷：
- 随机优化：依赖精确概率分布，对分布误判敏感。
- 鲁棒优化：过度保守，仅考虑最坏场景。
- CVaR单独应用：需预设分布形式，无法处理分布模糊性。
二、核心模型：基于矩的DRO与CVaR融合框架
1. 分布鲁棒优化（DRO）基本原理
- 模糊集构建：基于历史电价数据估计均值 μ 和协方差矩阵 Σ，定义矩不确定集：
其中 ξ 为随机电价向量，M 为所有概率分布的集合。
- 目标函数：优化最坏情况下的期望成本：
X 为调度决策空间（如机组启停、出力）。
2. 条件风险价值（CVaR）的嵌入
- 定义：在置信水平 ββ 下，CVaR 表示损失超过VaR的条件期望：
L 为损失函数（如负利润）。
- 与DRO结合：将CVaR作为目标或约束，在模糊集 PP 下优化最坏CVaR：
该形式同时捕捉分布模糊性和尾部风险。
三、模型求解：优化维度降低（ODR）与ADMM算法
1. 计算瓶颈与ODR方法
- 问题：基于矩的DRO常转化为半正定规划（SDP），维度随电价变量数指数增长。
- ODR核心思想：利用随机参数的低秩特性，将高维SDP分解为低维子问题：
- 证明原问题矩阵秩上界为 r≪dim(ξ)r≪dim(ξ)。
- 将维度降至 rr 维，构造内外两层近似SDP：
- 外层提供下界（保守解）
- 内层提供上界（激进解）
- 通过迭代逼近全局最优。
2. ADMM求解流程
对降维后的非凸SDP问题，采用改进ADMM算法：
- 变量拆分：将决策变量 xx 和辅助变量 yy 解耦。
- 交替更新：
其中 Lρ 为增广拉格朗日函数。
3. 收敛性：在矩约束下可收敛至 ϵϵ-最优解。
性能优势：
- 计算时间减少 3个数量级（如从小时级到秒级）。
- 解与全局最优差距 < 0.1%。
四、电价不确定性建模与场景生成
1. 概率分布建模
- 参数法：假设电价服从 tt-位置尺度分布（位置参数 ytyt，尺度参数 σtσt），适用于厚尾特征。
- 非参数法：核密度估计（KDE）拟合历史数据，避免分布假设偏差。
2. 场景生成与削减
- 蒙特卡洛采样：生成 NN 个电价场景。
- Kantorovich距离削减：保留 SS 个典型场景（如 S=30S=30），最小化与原场景集的距离：
确保计算效率。
五、IEEE节点系统验证案例
1. IEEE 118节点系统
- 结构：54台发电机、186条支路、99个负荷节点，分区运行。
- 自调度问题：
- 决策变量：机组启停 utiuti、出力 ptipti、备用容量 rtirti。
- 目标函数（结合DRO-CVaR）：
为为燃料成本，
为电价。
- 约束：功率平衡、爬坡率、备用要求。
2. 关键结果
- 风险控制：当 β=0.95 时，尾部损失降低 23%（对比随机规划）。
- 计算效率：ODR-ADMM求解时间 < 10秒，传统SDP超1小时。
- 保守性改善：矩模糊集比鲁棒优化成本降低 12%。
六、对比与拓展
1. 与其他方法对比
| 方法 | 优点 | 缺点 |
|---|---|---|
| DRO-CVaR（本模型） | 平衡保守性与风险控制 | 需矩信息估计 |
| 随机规划 | 解析方便 | 依赖精确分布 |
| 鲁棒优化 | 计算简单 | 过度保守 |
| 传统CVaR | 尾部风险控制 | 忽略分布模糊性 |
2. 工业应用扩展
- 微电网投标：结合光伏/风电不确定性。
- 碳排放约束：引入碳成本函数。
- 多代理调度：分布式ADMM协调多个发电商。
结论
本文提出的 基于矩的DRO-CVaR模型 通过模糊集刻画电价分布不确定性，利用CVaR量化尾部风险，结合ODR-ADMM算法实现高效求解。在IEEE 118等系统中验证表明：
- 在同等风险水平下，预期成本降低8%~15%。
- 极端场景损失减少 >20%。
- 计算效率满足日前市场时间要求（秒级求解）。
未来方向包括：融合数据驱动的Wasserstein模糊集、结合深度学习实时预测矩参数、拓展至多能源耦合系统。
📚2 运行结果
部分代码：
T = SCUC_data.totalLoad.T; % 时段数T
G = SCUC_data.units.N; % 发电机数
N = SCUC_data.baseparameters.busN; % 节点总数
all_branch.I = [ SCUC_data.branch.I; SCUC_data.branchTransformer.I ]; %所有支路起点 前是支路起点 后是变压器支路起点
all_branch.J = [ SCUC_data.branch.J; SCUC_data.branchTransformer.J ]; %所有支路终点
all_branch.P = [ SCUC_data.branch.P; SCUC_data.branchTransformer.P ]; %支路功率上限
beta_CVaR = 0.95;
% 形成直流潮流系数矩阵B
type_of_pf = 'DC';
Y = SCUC_nodeY(SCUC_data,type_of_pf);
B = -Y.B; %因为是直流方程 所以B忽略了电阻 只考虑电抗
% 定义一些方便用的常量
diag_E_T = sparse(1:T,1:T,1); %T*T的对角阵，对角线全1
low_trianle = diag_E_T(2:T,:) - diag_E_T(1:T-1,:);
% 定义变量 y 和 z
for i=1:N % 逐 个 处理每个节点
if ismember(i, SCUC_data.units.bus_G) % 如果是发电机节点 决策变量包括 功率，相角，功率约束z
y{i}.PG = sdpvar(T,1); %sdpvar创建实数型决策变量
y{i}.theta = sdpvar(T,1);
y{i}.z = sdpvar(T,1);
else
y{i}.theta = sdpvar(T,1);
end
end
for i=1:N % 逐 个 处理每个节点
if ismember(i, SCUC_data.units.bus_G) % 如果是发电机节点 决策变量包括 功率，相角，功率约束z
z{i}.PG = sdpvar(T,1); %sdpvar创建实数型决策变量
z{i}.theta = sdpvar(T,1);
z{i}.z = sdpvar(T,1);
else
z{i}.theta = sdpvar(T,1);
end
end
% end of 定义变量
% 下面计算 miu_hat均值的估计值 sigema_hat协方差的估计值
miu_hat_G_T = zeros(G,T); %每个机组每个时段电价都不同，G*T个
for q = 1:q_line % q_line标记样本长度
miu_hat_G_T = miu_hat_G_T + reshape(lamda_q_NT(q,:,:),G,T); %按样本数q依次累加 G*T矩阵型的电价
end
miu_hat_G_T = 1/q_line * miu_hat_G_T;
miu_hat = reshape(miu_hat_G_T',G*T,1); %按每个机组所有时段排列为列向量
% Sigma_hat = sparse(1:G*T,1:G*T,1);
% Sigma_hat = full(Sigma_hat);
Sigma_hat = zeros(G*T,G*T);
for q = 1:q_line
tmp1 = reshape(lamda_q_NT(q,:,:),G,T) - miu_hat_G_T;
tmp2 = reshape(tmp1', G*T,1);
Sigma_hat = Sigma_hat + tmp2 * tmp2';
end
Sigma_hat = 1/q_line * Sigma_hat;
Sigma_hat_neg_half = Sigma_hat^(-1/2);
% end of 均值 协方差计算
% 构造模糊集S G*T维的列向量，按每个G所有时段排列
tmp = max(lamda_q_NT,[],1); %取各个时段，所有样本中G个机组的最大值
% tmp = max(lamda_q_NT(:,:,1),[],1); %取第一时段，所有样本G个机组的最大值
tmp = reshape(tmp,G,T)';
lamda_positive = reshape(tmp, G*T,1); % 计算 lamda正
tmp = min(lamda_q_NT,[],1);
% tmp = min(lamda_q_NT(:,:,1),[],1);
tmp = reshape(tmp,G,T)';
lamda_negative = reshape(tmp, G*T,1); % 计算 lamda负
% end of 构造模糊集S
%调用切分函数切分区域
[PI,PINumber,PIG] = portion(N);
D = length(PINumber); %确定划分块数
%区分各节点是内节点还是外节点
ext = [];
for d = 1:D
for i = 1:size(all_branch.I,1) %从第一条支路开始循环遍历所有支路
left = all_branch.I(i); %支路起点和终点即可得到B电纳
right = all_branch.J(i);
if ismember(left, PI{d}') && ~ismember(right, PI{d}') %判断该支路是否属于当前分区
ext = [ext;left;right];
end
end
end
ext = unique(ext);
%定义其他变量
r = sdpvar(D,1);
t = sdpvar(D,1);
alpha_CVaR = sdpvar(D,1);
z_vector = [];
P_vector = [];
for g = 1:G
z_vector = [z_vector; y{SCUC_data.units.bus_G(g)}.z ];
P_vector = [P_vector; y{SCUC_data.units.bus_G(g)}.PG ]; %列向量用;分隔
end
%设置ADMM参数
u_A_b = [];
rho = 5;
gap_pri = 1e-2;
gap_dual = 1e-2;
M = 100; %迭代次数
% core = D;
% p = parpool(core);
% p.IdleTimeout = 100;
for k = 1:M
p_A = [];
z_A = [];
z_A_pri = [];
u_A = [];
p_A_k1 = [];
u_r = 0; %标记分区变量长度
%p-update
for d = 1:D
PIi = PI{d};
PIN = PINumber{d};
d_g = PIG{d};
Constraints = []; %每次将约束集置空
miu_hat_d = [];
lamda_positive_d = [];
lamda_negative_d = [];
y_d = [];
zy_d = [];
theta_d = [];
ztheta_d = [];
z_d = [];
for i = PIi' %循环处理当前片区节点
theta_d = [theta_d;y{i}.theta]; %取当前分区所有节点相角
ztheta_d = [ztheta_d;z{i}.theta]; %取当前分区所有节点相角
if ismember(i, SCUC_data.units.bus_G)
y_d = [y_d;y{i}.PG]; %取当前片区所有电机发电量
zy_d = [zy_d;z{i}.PG]; %取当前片区所有电机发电量
z_d = [z_d;y{i}.z]; %取当前片区电机费用
i_g = find(SCUC_data.units.bus_G == i); %求当前电机是第几个电机
% 机组出力上下界
Constraints = [Constraints, ...
SCUC_data.units.PG_low(i_g) * ones(T,1) <= y{i}.PG <= SCUC_data.units.PG_up(i_g) * ones(T,1) ];
% 爬坡 Pup和Pdown相等 且不考虑P0
Constraints = [Constraints, ...
-SCUC_data.units.ramp(i_g) * ones(T-1,1) <= low_trianle * y{i}.PG <= SCUC_data.units.ramp(i_g) * ones(T-1,1) ];
% 发电费用线性化
for l = 0:(L-1)
p_i_l =SCUC_data.units.PG_low(i_g) + ( SCUC_data.units.PG_up(i_g) - SCUC_data.units.PG_low(i_g) ) / L * l;
Constraints = [Constraints, ...
(2* p_i_l * SCUC_data.units.gamma(i_g) +SCUC_data.units.beta(i_g) ) * y{i}.PG - y{i}.z <= (p_i_l^2 * SCUC_data.units.gamma(i_g) - SCUC_data.units.alpha(i_g)) * ones(T,1) ];
end
%计算分区电价均值
miu_hat_m = zeros(1,T);
for q = 1:q_line % q_line标记样本长度
miu_hat_m = miu_hat_m + reshape(lamda_q_NT(q,i_g,:),1,T); %按样本数q依次累加第i个机组的电价
end
miu_hat_m = 1/q_line * miu_hat_m;
miu_hat_m = reshape(miu_hat_m',T,1); %按每个机组所有时段排列为列向量
miu_hat_d = [miu_hat_d ; miu_hat_m];
%取对应电机最大最小值
lamda_positive_m = lamda_positive(1+T*(i_g-1):T*i_g,1);
lamda_positive_d = [lamda_positive_d;lamda_positive_m];
lamda_negative_m = lamda_negative(1+T*(i_g-1):T*i_g,1);
lamda_negative_d = [lamda_negative_d;lamda_negative_m];
end
if ~ismember(i, ext) %判断节点是否是内部节点
%构造分区内部的直流潮流不等式
%构造中间项
cons_PF = sparse(T,1); %一个t一个约束 T行约束
if ismember(i, SCUC_data.units.bus_G) %如果i是发电机节点，则中间项应该加上PG
cons_PF = cons_PF + y{i}.PG;
end
for j = 1:N
cons_PF = cons_PF -B(i,j) .* y{j}.theta;
end
%构造右侧项
%find函数返还的是下标号
index_loadNode = find(SCUC_data.busLoad.bus_PDQR==i); % 节点i是否为负荷节点
if index_loadNode>0
b_tmp = SCUC_data.busLoad.node_P(:,index_loadNode); %按下标取出该负荷节点负载
else
b_tmp = sparse(T,1);
end
Constraints = [Constraints, ...
0 <= cons_PF <= b_tmp ];
%end of 构造分区内部的支流潮流不等式
end🎉3 参考文献
文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。(文章内容仅供参考，具体效果以运行结果为准)
[1] R. A. Jabr, "Robust self-scheduling under price uncertainty using conditional value-at-risk,"
IEEE Trans. Power Syst., vol. 20, no. 4, pp. 1852-1858, Nov. 2005.🌈4 Matlab代码、数据、文章下载
资料获取，更多粉丝福利，MATLAB|Simulink|Python资源获取
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
