👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆下载资源链接👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆
《《《《《《《《更多资源还请持续关注本专栏》》》》》》》
论文与完整源程序_电网论文源程序的博客-CSDN博客
https://blog.csdn.net/liang674027206/category_12531414.html
请注意使用该代码前需要对主程序以及调用函数的路径改成自己的！！！！！
高比例可再生能源接入对电力系统调峰能力提出了更高要求，完善调峰辅助服务成本分摊 机制有利于促进可再生能源消纳。文中从净负荷波动的角度出发，建立了调峰成本的量化与分摊模型，构造了无调峰需求的替代场景，将负荷和可再生能源出力曲线分别转换为无波动的均值线。其次，建立了含深度调峰和抽水蓄能的调度优化模型，用于计算不同场景下的调峰成本，并将有无调峰需求两种场景下的系统调峰成本之差作为单一主体导致的边际调峰成本，然后采用Shapley值计算不同主体导致的调峰成本。最后，根据成本的引发程度分摊调峰成本。算例表明，提出的调峰成本量化模型能够反映不同主体导致的调峰成本或贡献的调峰价值，成本分摊机制能够传递调峰成本信息。
部分代码展示：
% function []= Wind_PV_b2( )
clc
clear
close all
addpath('E:\code\78-高比例可再生能源电力系统的调峰成本量化与分摊模型\第二步多主体集合调用程序包\函数包');
load('kONOFF');
addpath('E:\code\78-高比例可再生能源电力系统的调峰成本量化与分摊模型\第一步整理源荷数据');
load('E_PV');load('E_Wind');load('Eload');load('D_Eload_FW');load('D_Eload_JM');load('D_Eload_ZZ');
% D_Eload_JM=zeros(1,24);
% D_Eload_FW=zeros(1,24);
% D_Eload_ZZ=zeros(1,24);
load('L_jun');
load('kWind');
load('kPV');
E_PV=kPV*E_PV;
E_Wind=kWind*E_Wind;
%% 公式（6） 目标函数：调峰成本
Cost = sdpvar(1,1); %系统总调峰成本
%公示6在编写的时候，需要将i属于Gc这种，将其中的属于号，用向量或者矩阵运算给等效
%此外，这里目标函数仅仅只是调峰成本，是不是应该将目标函数定为煤炭成本+调峰成本更合适呢（又或者是文中无弃风弃光，无失负荷，所以总发电量必然不变，等效为一阶线性模型总煤耗量不变，所以就没必要算呢）
%为了方便功率累加写电平衡等式，这里将功率和调峰价格的基础部分也占用一行
lamda300_Do = repmat([600;782;833;936]',18,1,24); %为 机 组 i 在 第 o 段 深 度 调 峰 的 报 价,分3段，即35%。多加了第一行基础调峰的范围。
lamda600_Do = repmat([550;731;780;883]',18,1,24); %0.3元一度电的成本价格，售价基准价格按照0.6元/度计算，加上深度调峰的惩罚价格
%注意，这里的深度调峰功率档位是互斥事件，即最多只有一个是非零量，这样更好编程。
%这里在犹豫要不要用三维的矩阵去写【机组序号，时刻，各调峰区间出力】的代码，不知可不可行
P300_Diot = sdpvar(18,4,24,'full');
P600_Diot = sdpvar(18,4,24,'full');
y300_it = binvar(18,24,'full'); % 机组i在时段t是否开启动作的0-1变量
z300_it = binvar(18,24,'full'); % 机组i在时段t是否停机动作的0-1变量
y600_it = binvar(18,24,'full');
z600_it = binvar(18,24,'full');
%200MW级别的火电机组不参与深度调峰（见表A2）
P200_it = sdpvar(20,24,'full');
y200_it = binvar(20,24,'full');
z200_it = binvar(20,24,'full');
%火电机组的话停机重启一次成本比较高，66万千瓦机组一次大概一百万人民币左右。https://www.zhihu.com/question/330394129
C_SU300 = 400000*kONOFF; %算例源荷用的是MW单位，调峰成本用万元做单位
C_SD300 = 100000*kONOFF;
C_SU600 = 800000*kONOFF;
C_SD600 = 200000*kONOFF;
C_SU200 = 300000*kONOFF;
C_SD200 = 100000*kONOFF;
%公式6只计算抽水蓄能电站的抽水成本，没计算其发电成本（我猜是原作认为抽水相当于加负荷，可以缓解深度调峰压力）
%算例分析中指出，抽 水 蓄 能 电 站 的 综 合 效 率 为 0.75，报 价 为100元/ ( MW·h )。
C_PuS = 100*1;
P_PuS = sdpvar(4,24,'full');
%风光削减惩罚
V_LE = 600*1; %单位能源削减成本
P_PV = sdpvar(1,24,'full');
P_Wind = sdpvar(1,24,'full');
P_PVcur = sdpvar(1,24,'full');
P_Windcur = sdpvar(1,24,'full');
C=[];
%火电机组深度调峰费用
%火电机组启停费用
%抽水蓄能调峰费用
%削减风光的补偿费用
P_guding = sdpvar(1,1); %固定出力
load('kGD.mat','kGD');
% %% 修改②删去弃光惩罚
C=[C,
%Cost==sum(sum(V_LE*(P_Windcur+P_PVcur)))+sum(sum(sum(lamda300_Do.*P300_Diot+lamda600_Do.*P600_Diot)))+sum(sum(y300_it*C_SU300+z300_it*C_SD300+y600_it*C_SU600+z600_it*C_SD600))+sum(sum(y200_it*C_SU200+z200_it*C_SD200))+sum(sum(C_PuS*P_PuS)),
Cost==sum(sum(sum(lamda300_Do.*P300_Diot+lamda600_Do.*P600_Diot)))+700*sum(sum(P200_it))+sum(sum(y300_it*C_SU300+z300_it*C_SD300+y600_it*C_SU600+z600_it*C_SD600))+sum(sum(y200_it*C_SU200+z200_it*C_SD200))+sum(sum(C_PuS*P_PuS))+200*sum(P_PV+P_Wind+P_guding)+sum(sum(V_LE*(P_PVcur+P_Windcur))),
P_guding == 5000*kGD,
];
%反思一下，原文里的火电机组深度调峰的成本应该是阶梯成本，但也没完全说清楚，这里为了简化运算，就直接用常值成本来计算了，因此，50-100%之间的单位成本最好也给算上
%% 公式（7） 电平衡约束
%需要注意的是，AA(i,j,t)，无法通过AA(1,1,1:24)的方式，调取出一个向量
%只能通过AA(1,1,1)的方式，调取出一个数值，范例如下：
% A = ones(18,4,24);
% AA = sum(sum(A(:,:,1)));
% AAA=AA(1,1,1:24);
P_PuG = sdpvar(4,24,'full'); %抽水蓄能发电
%因此，必须得用一个1-24的循环来写各时刻的电平衡等式约束了。
%注意sum(sum(A))后的维度是1*1*24,注意sum(sum(A(:,:,1)))后的维度是1*1,
% 论文里用的是固定功率，但主体个数太少的话，就会出现负荷小于固定出力，因此这里，将固定出力设置为单价较低的处理范围
for t=1:24
C=[C, sum(sum(P300_Diot(:,:,t)))+sum(sum(P600_Diot(:,:,t)))+sum(P200_it(:,t))+P_PV(t)+P_Wind(t)+sum(P_PuG(:,t)-P_PuS(:,t))+P_guding == L_jun(t)+D_Eload_FW(t)+D_Eload_JM(t)+D_Eload_ZZ(t),
];
end
%% 公式（8-9）系统正负旋转备用约束
u200Git = binvar(20,24,'full'); %uit 表示机组运行状态的0-1变量,1开机，0关机
u300Git = binvar(18,24,'full');
u600Git = binvar(18,24,'full');
R_Ureq = (Eload*0.05+E_Wind*0.05+E_PV*0.05); %正旋转备用
R_Dreq = (Eload*0.05+E_Wind*0.05+E_PV*0.05); %负旋转备用
for t=1:24
C=[C, 200*sum(u200Git(:,t))+300*sum(u300Git(:,t))+600*sum(u600Git(:,t))-sum(sum(P300_Diot(:,:,t)))-sum(sum(P600_Diot(:,:,t)))-sum(P200_it(:,t)) >= R_Ureq(t), %正旋转备用约束
sum(sum(P300_Diot(:,:,t)))+sum(sum(P600_Diot(:,:,t)))+sum(P200_it(:,t))-200*0.5*sum(u200Git(:,t))-300*0.35*sum(u300Git(:,t))-600*0.35*sum(u600Git(:,t)) >= R_Dreq(t), %负旋转备用约束
];
end
%% 公式（10-12）火电机组功率大小约束，这里换一种方法去写
uP300_Diot = binvar(18,4,24,'full');
uP600_Diot = binvar(18,4,24,'full');
for t= 1:24
for i = 1:18
C = [C,
sum(uP300_Diot(i,:,t))<= u300Git(i,t),
sum(uP600_Diot(i,:,t))<= u600Git(i,t),
u200Git(i,t)*200*0.5<=P200_it(i,t),P200_it(i,t)<=u200Git(i,t)*200,
uP300_Diot(i,1,t)*300*0.5<=P300_Diot(i,1,t),P300_Diot(i,1,t)<=uP300_Diot(i,1,t)*300,
uP300_Diot(i,2,t)*300*0.45<=P300_Diot(i,2,t),P300_Diot(i,2,t)<=uP300_Diot(i,2,t)*300*0.5,
uP300_Diot(i,3,t)*300*0.4<=P300_Diot(i,3,t),P300_Diot(i,3,t)<=uP300_Diot(i,3,t)*300*0.45,
uP300_Diot(i,4,t)*300*0.35<=P300_Diot(i,4,t),P300_Diot(i,4,t)<=uP300_Diot(i,4,t)*300*0.4,
uP600_Diot(i,1,t)*600*0.5<=P600_Diot(i,1,t),P600_Diot(i,1,t)<=uP600_Diot(i,1,t)*600,
uP600_Diot(i,2,t)*600*0.45<=P600_Diot(i,2,t),P600_Diot(i,2,t)<=uP600_Diot(i,2,t)*600*0.5,
uP600_Diot(i,3,t)*600*0.4<=P600_Diot(i,3,t),P600_Diot(i,3,t)<=uP600_Diot(i,3,t)*600*0.45,
uP600_Diot(i,4,t)*600*0.35<=P600_Diot(i,4,t),P600_Diot(i,4,t)<=uP600_Diot(i,4,t)*600*0.4,
];
end
for i = 19:20
C = [C, u200Git(i,t)*200*0.5<=P200_it(i,t),P200_it(i,t)<=u200Git(i,t)*200,
];
end
end
%公式（11-12）不用再写了，已经隐含在上式（10）里了。
%% 公式（13-14）
%约束的含义是，持续运行需要受到上下爬坡约束，但启停不受爬坡约束。火电机组爬坡率为2%/min,小时调度体现不出来
for i= 1:18
for t = 2:24
C = [C,
sum(P300_Diot(i,:,t))-sum(P300_Diot(i,:,t-1)) <= 0.3*300+(1-u300Git(i,t-1))*300,
sum(P300_Diot(i,:,t-1))-sum(P300_Diot(i,:,t)) <= 0.3*300+(1-u300Git(i,t))*300,
sum(P600_Diot(i,:,t))-sum(P600_Diot(i,:,t-1)) <= 0.3*600+(1-u600Git(i,t-1))*600,
sum(P600_Diot(i,:,t-1))-sum(P600_Diot(i,:,t)) <= 0.3*600+(1-u600Git(i,t))*600,
P200_it(i,t)-P200_it(i,t-1) <= 0.3*200+(1-u200Git(i,t-1))*200,
P200_it(i,t-1)-P200_it(i,t) <= 0.3*200+(1-u200Git(i,t))*200,
];
end
for i = 19:20
C = [C, P200_it(i,t)-P200_it(i,t-1) <= 0.3*200+(1-u200Git(i,t-1))*200,
P200_it(i,t-1)-P200_it(i,t) <= 0.3*200+(1-u200Git(i,t))*200,
];
end
end
%% 公式（15-16） 火电机组的启停状态标识位置
C = [C,
y200_it(:,2:24)-z200_it(:,2:24)==u200Git(:,2:24)-u200Git(:,1:23), %公式(15)
y300_it(:,2:24)-z300_it(:,2:24)==u300Git(:,2:24)-u300Git(:,1:23),
y600_it(:,2:24)-z600_it(:,2:24)==u600Git(:,2:24)-u600Git(:,1:23),
y200_it(:,1)==0,y300_it(:,1)==0,y600_it(:,1)==0,
z200_it(:,1)==0,z300_it(:,1)==0,z600_it(:,1)==0,
y200_it+z200_it<=1,%公式(16)
y300_it+z300_it<=1,
y600_it+z600_it<=1,
];
%% 公式（17-18） 最短运行时间，最短停机时间限值，取20h和10h
%需要注意的是，调度的维度只有24，因此，编写代码的时候，需要避免调用越界24
Ton = 20;%最短开机时间
Toff = 10;%最短关机时间
for t=1:(25-Ton) %公式(17)
for i=1:18
C = [C,sum(u200Git(i,t:(t+Ton-1))) >= Ton*y200_it(i,t)];
C = [C,sum(u300Git(i,t:(t+Ton-1))) >= Ton*y300_it(i,t)];
C = [C,sum(u600Git(i,t:(t+Ton-1))) >= Ton*y600_it(i,t)];
end
for i=19:20
C = [C,sum(u200Git(i,t:(t+Ton-1))) >=Ton*y200_it(i,t)];
end
end
for t=1:(25-Toff) %公式(18)
for i=1:18
C = [C,sum(1-u200Git(i,t:(t+Toff-1))) >= Toff*z200_it(i,t),
sum(1-u300Git(i,t:(t+Toff-1))) >= Toff*z300_it(i,t),
sum(1-u600Git(i,t:(t+Toff-1))) >= Toff*z600_it(i,t)];
end
for i=19:20
C = [C,sum(1-u200Git(i,t:(t+Toff-1))) >=Toff*z200_it(i,t)];
end
end
%% 公式（19） 可再生能源出力约束
C = [C, 0<=P_PV,P_PV<=E_PV,
0<=P_PVcur,P_PVcur<=E_PV,
P_PV+P_PVcur==E_PV,
0<=P_Wind,P_Wind<=E_Wind,
0<=P_Windcur,P_Windcur<=E_Wind,
P_Wind+P_Windcur==E_Wind,
];
%% 公式（20-24） 抽水蓄能电站功率约束
U_PuG = binvar(4,24,'full');
U_PuS = binvar(4,24,'full');
u_PuG = binvar(1,24,'full');
u_PuS = binvar(1,24,'full');
C = [C,
U_PuG*300*0.3<=P_PuG,P_PuG<=U_PuG*300, %公式（20）
U_PuS*300*0.3<=P_PuS,P_PuS<=U_PuS*300, %公式（21）
u_PuG+u_PuS<=1, %公式（22）
];
for t=1:24
C = [C,
U_PuG(:,t)<=u_PuG(1,t), %公式（23）
U_PuS(:,t)<=u_PuS(1,t), %公式（24）
];
end
%% 公式（25-27） 抽水蓄能电站容量约束
%抽水蓄能电站的综合效率为0.75，取抽水效率为0.9 ,发电效率为0.833 ，nG=1.2
%抽水蓄能可以直接用电量表示容量，也可以用蓄水量代表容量
%https://baijiahao.baidu.com/s?id=1709750463608980862&wfr=spider&for=pc
%仔细看公式（25）可知，文中的抽水蓄能是四个抽水机，四个发电机，共用一个蓄水池
nG=1.2;
nS=0.9;
QU = sdpvar(1,24,'full');
C = [C,
QU(1,2:24)==QU(1,1:23)+nS*sum(P_PuS(:,2:24))-nG*sum(P_PuG(:,2:24)), %公式（25）
300*4*1<=QU,QU<=300*4*5, %公式（26）
nG*sum(sum(P_PuG))==nS*sum(sum(P_PuS)), %始末容量相等约束%公式（27）
% 用这两个公式表示始末容量相等，也行的
% QU(1,1)== 3000*2+nS*sum(P_PuS(:,1))-nG*sum(P_PuG(:,1)),
% QU(1,24) == 3000*2,
];
%% 调用cplex求解器求解
ops=sdpsettings('verbose', 2, 'solver', 'cplex');
% ops = sdpsettings('solver', 'cplex+', 'verbose', 2, 'debug', 1, 'cplex.NonConvex', -1);
load('mipgap');
ops.cplex.MIPGap =mipgap;
ops.cplex.TimeLimit = 200;
ops.cplex.MIPFocus=3;
ops.cplex.Presolve=2;
% ops.cplex.DualReductions = 0;
%进行求解计算
result=optimize(C,Cost,ops);
if result.problem==0
fprintf('求解成功 ！！！\n');
else
error('求解出错！请查找错误来源');
end
% [model,recoverymodel] = export(C,Cost,sdpsettings('solver','cplex'));
% iis = cplex_iis(model);
% cplex_write(model, 'TestModel.lp');
% sum(iis.Arows)
% find(iis.Arows==1)-1
%% 绘图
Cost = value(Cost);
%% 电平衡图
P200_it=value(P200_it);
P300_Diot=value(P300_Diot);
P600_Diot=value(P600_Diot);
P_PV=value(P_PV);
P_Wind=value(P_Wind);
P_PuS=value(P_PuS);
P_PuG=value(P_PuG);
f_P_PuS=sum(P_PuS);
f_P_PuG=sum(P_PuG);
f_P200_it = sum(P200_it);
for t=1:24
f_P300_Diot(1,t) = sum(sum(P300_Diot(:,:,t)));
f_P600_Diot(1,t) = sum(sum(P600_Diot(:,:,t)));
end
figure
Plot_E_yuan=zeros(7,24);
for t=1:24
Plot_E_yuan(1,t)=P_guding;
Plot_E_yuan(2,t)=f_P200_it(t);
Plot_E_yuan(3,t)=f_P300_Diot(t);
Plot_E_yuan(4,t)=f_P600_Diot(t);
Plot_E_yuan(5,t)=P_PV(t);
Plot_E_yuan(6,t)=P_Wind(t);
Plot_E_yuan(7,t)=f_P_PuG(t);
end
Plot_E_fuhe =zeros(2,24);
for t=1:24
Plot_E_fuhe(1,t) = -f_P_PuS(t);
Plot_E_fuhe(2,t) = -(L_jun(t)+D_Eload_JM(t)+D_Eload_FW(t)+D_Eload_ZZ(t));
end
bar(Plot_E_yuan','stacked');
hold on
bar(Plot_E_fuhe','stacked');
xlabel('时间/h');
ylabel('功率/MW');
title('电功率平衡图');
legend('+廉价固定出力','+200MW火电机组','+300MW火电机组','+600MW火电机组','+光电','+风电','+抽水蓄能发电', ...
'-抽水蓄能耗电','-电负荷耗电');
xlim([0 32 ]);
set(gca,'FontSize',12 );
%% 处于深度调峰位置的火电机组出力标识位绘图
uP300_Diot=value(uP300_Diot);
uP600_Diot=value(uP600_Diot);
Nu300Git=zeros(18,24);
Nu600Git=zeros(18,24);
for i = 1:18
for t=1:24
Nu300Git(i,t)=sum(uP300_Diot(i,:,t)*[1;2;3;4]);
Nu600Git(i,t)=sum(uP600_Diot(i,:,t)*[1;2;3;4]);
end
end
figure
[XX,YY] = meshgrid(1:24,1:18 );
mesh(XX,YY,Nu600Git);
xlabel('时刻（h）');
ylabel('机组序号');
zlabel('（0是关机，1是基础调峰，2、3、4分别是深度调峰1、2、3挡）');
title('600MW机组的深度调峰档位图');
figure
[XX,YY] = meshgrid(1:24,1:18 );
mesh(XX,YY,Nu300Git);
xlabel('时刻（h）');
ylabel('机组序号');
zlabel('（0是关机，1是基础调峰，2、3、4分别是深度调峰1、2、3挡）');
title('300MW机组的深度调峰档位图');
%% 四个主体时的总成本是
fprintf('五个主体时的总成本是%d\n',Cost);
%% 四个主体时的深度调峰成本是
P300_Diot=value(P300_Diot);
P600_Diot=value(P600_Diot);
y300_it=value(y300_it);
z300_it=value(z300_it);
y600_it=value(y600_it);
z600_it=value(z600_it);
y200_it=value(y200_it);
z200_it=value(z200_it);
P_PVcur=value(P_PVcur);
P_Windcur=value(P_Windcur);
P_guding=value(P_guding);
deep_lamda300_Do = repmat([600;782;833;936]'-600,18,1,24); %为 机 组 i 在 第 o 段 深 度 调 峰 的 报 价,分3段，即35%。多加了第一行基础调峰的范围。
deep_lamda600_Do = repmat([550;731;780;883]'-550,18,1,24); %0.3元一度电的成本价格，售价基准价格按照0.6元/度计算，加上深度调峰的惩罚价格
deepCost=sum(sum(sum(deep_lamda300_Do.*P300_Diot+deep_lamda600_Do.*P600_Diot)))+sum(sum(y300_it*C_SU300+z300_it*C_SD300+y600_it*C_SU600+z600_it*C_SD600))+sum(sum(y200_it*C_SU200+z200_it*C_SD200))+sum(sum(C_PuS*P_PuS))+sum(sum(V_LE*(P_PVcur+P_Windcur)));
fprintf('五个主体时的深度调峰成本成本是%d\n',deepCost);
%% 需要分摊的深度调峰成本总量是：
FentanSum_deepCost = deepCost;
save('E:\code\78-高比例可再生能源电力系统的调峰成本量化与分摊模型\第三步求取五个波动主体时的深度调峰成本\FentanSum_deepCost.mat','FentanSum_deepCost');
fprintf('需要分摊的深度调峰成本总量是：%d\n',FentanSum_deepCost);
% end
效果展示：
75号资源-复现源程序：论文可在知网下载《高比例可再生能源电力系统的调峰成本量化与分摊模型》本人博客有解读资源-CSDN文库https://download.csdn.net/download/LIANG674027206/89139140👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆下载资源链接👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆
《《《《《《《《更多资源还请持续关注本专栏》》》》》》》
论文与完整源程序_电网论文源程序的博客-CSDN博客https://blog.csdn.net/liang674027206/category_12531414.html
810
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
