主要内容
程序复现文献《高比例可再生能源电力系统的调峰成本量化与分摊模型》，从净负荷波动的角度出发，建立了调峰成本的量化与分摊模型，构造了无调峰需求的替代场景，将负荷和可再生能源出力曲线分别转换为无波动的均值线。其次，建立了含深度调峰和抽水蓄能的调度优化模型，用于计算不同场景下的调峰成本，并将有无调峰需求两种场景下的系统调峰成本之差作为单一主体导致的边际调峰成本，然后采用 Shapley值计算不同主体导致的调峰成本。最后，根据成本的引发程度分摊调峰成本。程序采用matlab+cplex编写，程序已经深度优化，求解速度很快，注释清晰，代码免费分享给大家学习参考！
部分代码
for t=1:24C=[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),== 0;== 0;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:24C=[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), %正旋转备用约束:,:,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:24for i = 1:18C = [C,:,t))<= u300Git(i,t),:,t))<= u600Git(i,t),=P200_it(i,t),P200_it(i,t)<=u200Git(i,t)*200,=P300_Diot(i,1,t),P300_Diot(i,1,t)<=uP300_Diot(i,1,t)*300,=P300_Diot(i,2,t),P300_Diot(i,2,t)<=uP300_Diot(i,2,t)*300*0.5,=P300_Diot(i,3,t),P300_Diot(i,3,t)<=uP300_Diot(i,3,t)*300*0.45,=P300_Diot(i,4,t),P300_Diot(i,4,t)<=uP300_Diot(i,4,t)*300*0.4,=P600_Diot(i,1,t),P600_Diot(i,1,t)<=uP600_Diot(i,1,t)*600,=P600_Diot(i,2,t),P600_Diot(i,2,t)<=uP600_Diot(i,2,t)*600*0.5,=P600_Diot(i,3,t),P600_Diot(i,3,t)<=uP600_Diot(i,3,t)*600*0.45,=P600_Diot(i,4,t),P600_Diot(i,4,t)<=uP600_Diot(i,4,t)*600*0.4,];endfor i = 19:20C = [C, u200Git(i,t)*200*0.5<=P200_it(i,t),P200_it(i,t)<=u200Git(i,t)*200,];endend%公式（11-12）不用再写了，已经隐含在上式（10）里了。公式（13-14）%约束的含义是，持续运行需要受到上下爬坡约束，但启停不受爬坡约束。火电机组爬坡率为2%/min,小时调度体现不出来for i= 1:18for t = 2:24C = [C,:,t))-sum(P300_Diot(i,:,t-1)) <= 0.3*300+(1-u300Git(i,t-1))*300,:,t-1))-sum(P300_Diot(i,:,t)) <= 0.3*300+(1-u300Git(i,t))*300,:,t))-sum(P600_Diot(i,:,t-1)) <= 0.3*600+(1-u600Git(i,t-1))*600,:,t-1))-sum(P600_Diot(i,:,t)) <= 0.3*600+(1-u600Git(i,t))*600,<= 0.3*200+(1-u200Git(i,t-1))*200,<= 0.3*200+(1-u200Git(i,t))*200,];endfor i = 19:20C = [C, P200_it(i,t)-P200_it(i,t-1) <= 0.3*200+(1-u200Git(i,t-1))*200,<= 0.3*200+(1-u200Git(i,t))*200,];endend公式（15-16） 火电机组的启停状态标识位置C = [C,:,2:24)-z200_it(:,2:24)==u200Git(:,2:24)-u200Git(:,1:23), %公式(15):,2:24)-z300_it(:,2:24)==u300Git(:,2:24)-u300Git(:,1:23),:,2:24)-z600_it(:,2:24)==u600Git(:,2:24)-u600Git(:,1:23),:,1)==0,y300_it(:,1)==0,y600_it(:,1)==0,:,1)==0,z300_it(:,1)==0,z600_it(:,1)==0,=1,%公式(16)=1,=1,公式（17-18） 最短运行时间，最短停机时间限值，取20h和10h%需要注意的是，调度的维度只有24，因此，编写代码的时候，需要避免调用越界24Ton = 20;%最短开机时间Toff = 10;%最短关机时间for t=1:(25-Ton) %公式(17)for i=1:18C = [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)];endfor i=19:20C = [C,sum(u200Git(i,t:(t+Ton-1))) >=Ton*y200_it(i,t)];endendfor t=1:(25-Toff) %公式(18)for i=1:18C = [C,sum(1-u200Git(i,t:(t+Toff-1))) >= Toff*z200_it(i,t),:(t+Toff-1))) >= Toff*z300_it(i,t),:(t+Toff-1))) >= Toff*z600_it(i,t)];endfor i=19:20C = [C,sum(1-u200Git(i,t:(t+Toff-1))) >=Toff*z200_it(i,t)];endend公式（19） 可再生能源出力约束C = [C, 0<=P_PV,P_PV<=E_PV,=P_PVcur,P_PVcur<=E_PV,==E_PV,=P_Wind,P_Wind<=E_Wind,=P_Windcur,P_Windcur<=E_Wind,==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,=P_PuG,P_PuG<=U_PuG*300, %公式（20）=P_PuS,P_PuS<=U_PuS*300, %公式（21）=1, %公式（22）for t=1:24C = [C,:,t)<=u_PuG(1,t), %公式（23）:,t)<=u_PuS(1,t), %公式（24）end
结果一览
“阅读原文”获取程序源码
