基于灰狼算法的综合能源系统多时间尺度优化调度（MATLAB程序）

function Obj = fun_cost(x,flag) if nargin<2 flag=0; end x=x'; Obj = zeros(2,1); PMT=x(1,49:72);%燃气轮 P_dianjie=x(1,25:48);%电解槽 P_tanbu=x(1,1:24);%碳捕 pv=[0 0 0 0 0 0.6 2.4 10.5 30 69 69.9 95.4 129.9 111 120.9 99 71.4 39.9 12.9 0.9 0 0 0 0]'; %%夏季典型日 % PW=[60 50 48 30 32 28 29 45 46 55 45 67 65 45 60 55 54 52 45 50 55 65 55 50];%wind % PV=[0 0 0 0 0 5 10.2 25 30 35 45 48 60 61 45 43 30 15 0 0 0 0 0 0 ];%PV % % P_H=[80 78 95 92 90 105 65 75 85 90 100 95 80 75 80 90 80 85 80 72 70 75 60 65 ];%热 % P_C=[25 20 35 30 35 34 40 45 55 60 65 70 80 90 110 120 100 90 125 110 100 90 25 40];%冷 %%电价 Grid_P_S=[0.21 0.21 0.21 0.21 0.21 0.39 0.39 0.39 0.65 0.65 0.65 0.65 0.65 0.39 0.39 0.39 0.65 0.65 0.65 0.39 0.39 0.39 0.39 0.21 ]; Grid_P_B=[0.31 0.31 0.31 0.31 0.31 0.47 0.47 0.47 0.71 0.71 0.71 0.71 0.71 0.47 0.47 0.47 0.71 0.71 0.71 0.47 0.47 0.47 0.47 0.31 ]; % % % %%冬季典型日 P_H=[125 120 123 125 126 128 110 100 95 75 72 68 70 75 76 78 79 80 90 95 105 108 110 115];%热 PW=[50 30 40 20 10 15 50 22 50 30 35 30 35 25 30 35 34 32 25 20 25 35 25 30];%wind PV=[0 0 0 0 0 0 10.2 25 30 35 45 46 47 48 40 33 20 10 0 0 0 0 0 0 ];%PV P_C=[25 20 35 30 35 34 40 45 55 60 65 70 80 90 110 120 100 90 125 110 100 90 25 40]./2*0;%冷%冷 %%%阶段一：优先安排清洁能源机组 PDD=[700,750,850,950,1000,1100,1150,1200,1300,1400,1450,1500,1400,1300,1200,1050,1000,1100,1200,1400,1300,1100,900,800]*0.3;%负荷 [Pcd1,P_heat]=P_Elect_config(x(1,73:120));%%SOC ->MW KCF P_cH2=zeros(1,24); %%安排光伏风电 PL=PDD-PV-PW; %%安排MT、储能 PLD=PL-PMT+P_tanbu+P_dianjie-Pcd1; if flag~=1 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%气体平衡 %电解槽产生的Ｈ２ P_dianjie_H=P_dianjie.*0.75.*3.41./0.342;%(电_气综合能源系统环境经济调度研究_何良策　　式中2-6) %燃气轮机消耗的气体　　（２－１７） % a=[1000]; % b=[16.19]; % c=[0.00048]; MT_gas=1000+16.16.*PMT+0.00048.*PMT.*PMT; %%甲烷化的Ｈ２ Ch4_gas=(P_dianjie_H-P_cH2)/2;%%燃气轮机接在天然气管网，甲烷化的气体先于ＭＴ机组消耗的对冲，然后输入天然气管 %甲烷化消耗的CO2量化计算 CO2_gas=(P_dianjie_H-P_cH2)/2; %预热回收 PXSL_cool=P_H*0.2;%%工业预热回收供热 %电锅炉 Guol_dian=(P_H-P_heat-PXSL_cool)*0.8/2.42;%%%%热电值是2.42 %%电空调的出力 % P_ecooler=P_C-PXSL_cool; %%阶段二，利用火电进行实时调度，考虑风电光伏的误差及备用容量+ %%不考虑火电 Pgrid=PLD+Guol_dian; Fobject=sum(Pgrid.*Grid_P_B); UC_G=[]; % [Fobject,UC_G]=syst_optimal(PLD+P_ecooler,25); % [Fobject,UC_G]=syst_optimal(PLD+P_ecooler+0.03*PV,25); %%消耗的气量 %%消耗的气量 Gas_tol=sum(MT_gas-Ch4_gas.*1.05)*0;%不考虑燃气 %%计算成本 %运行成本 cost1=Fobject+Gas_tol.*2.47; %%%%%%%%%%%%%%%%%%%%%%%%%%%CO2平衡 Pc_co2=P_tanbu.*3.25;%碳捕系统单位电量捕获的cO2系数为3.25 %%环境成本 % %火电机组排放及治理 % co2=sum(sum(UC_G,2)).*0.21-sum(Pc_co2)-sum(CO2_gas)-sum(PMT).*0.2; % v_co2=max([0,co2]);%%看是否全部处理，如果全部处理，则可取0 % cost2=(v_co2.*8.5);%kg/元 % %%Ｓｏ２　０．.８０　氮　０．８５ % %燃气轮机排放及治理 % % cost2=cost2+sum(PMT).*0.2*62.96;%kg/元 % cost3=sum(abs(PLD-sum(UC_G,2)'))-sum(P_tanbu)-sum(P_dianjie);%%弃DG量 % cost3_1=max([0,cost3]); %%弃DG量 %%阶段三：日内滚动调度 %%设计场景：利用碳捕集系统去抑制风光波动，即当风光波动时，调整碳捕系统的出力。 %%滚动调度以每5分钟滚动一次 %01-模拟风光误差 cost1= cost1+sum(P_tanbu)*0.015; P_tc_time=zeros(1,96); for tt=1:24 %光伏波动 PV=pv(1,tt).*ones(1,12) *0.03; %%碳捕出力扩展成15min PP=repmat(P_tanbu(1,tt),1,12); PP=PV+PP; P_tc_time=[P_tc_time,PP]; end %%重新计算碳捕的出力及费用 %%%%%%%%%%%%%%%%%%%%%%%%%%%CO2平衡 Pc_co2=P_tc_time.*3.25;%碳捕系统单位电量捕获的cO2系数为3.25 %%环境成本 cost_evr=(sum(sum(UC_G,2)).*0.21+sum(PMT).*0.2)*45-sum(pv)*0.97*45-sum(PMT).*0.39*45;%排放处罚 %机组排放及治理,二氧化碳存储后交易 co2=(sum(Pc_co2)./12-sum(CO2_gas))*25; % v_co2=max([0,co2]);%%看是否全部处理，如果全部处理，则可取0 %%----------------------计算碳排放权交易量------------------------ %%清洁能源+储能发电量+碳补系统处理的量+电解装置消耗的量 %%碳排放价格按照广州碳排放交易中心确定（45/吨）（http://www.cnemission.cn/article/jydt/） %%碳排放指标参考（排放指标配额技术指南发改委）https://www.sohu.com/a/442954915_651733 % 燃气轮机：０．３９ｔ／ＭＷ;燃煤机组０．９７ｔ／ＭＷ cost2=co2; %%Ｓｏ２　０．.８０　氮　０．８５ %燃气轮机排放及治理 % cost2=cost2+sum(PMT).*0.2*62.96;%kg/元 cost3=sum(P_tanbu)-sum(P_dianjie)*0;%%弃DG量 cost3_1=max([0,cost3]); %% 运维成本 % cost4=sum(pv+PW)*0.045; % cost4=cost4+sum(abs(Pcd1))*0.05; % cost4=cost4+sum(P_tc_time*0.025); % cost4=cost4+sum(P_heat)*0.015+sum(Guol_dian*0.012)+sum(PXSL_cool)*0.002; Obj(1) = cost1+cost3_1*150+cost_evr; Obj(2) =cost2; end if flag==1 cost4=0; %预热回收 PXSL_cool=P_H*0.2;%%工业预热回收供热 %电锅炉 Guol_dian=(P_H-P_heat-PXSL_cool)*0.8/2.42;%%%%热电值是2.42 %% %工业预热回收 PXSL_cool=P_H*0.2; Pgrid=PLD+Guol_dian; Fobject=sum(Pgrid.*Grid_P_B); UC_G=[]; %%计算成本 %运行成本 cost1=Fobject; %%%%%%%%%%%%%%%%%%%%%%%%%%%CO2平衡 Pc_co2=P_tanbu.*3.25;%碳捕系统单位电量捕获的cO2系数为3.25 %%弃DG量 %%阶段三：日内滚动调度 %%设计场景：利用碳捕集系统去抑制风光波动，即当风光波动时，调整碳捕系统的出力。 %%滚动调度以每5分钟滚动一次 %01-模拟风光误差 P_tc_time=[];PVV=[];PV1=[];PW11=[]; for tt=1:24 %风电光伏波动 PW1=PW(1,tt).*ones(1,12) *0.01; PV=pv(1,tt).*ones(1,12) *0.03; %%碳捕出力扩展成5min PP=repmat(P_tanbu(1,tt),1,12); PP=PV+PP+PW1; P_tc_time=[P_tc_time,PP]; PVV=[PVV,repmat(pv(1,tt),1,12)*1.03;]; PV1=[PV1,repmat(pv(1,tt),1,12)]; PW11=[PW11,repmat(PW(1,tt),1,12)]; end %%重新计算碳捕的出力及费用 %%%%%%%%%%%%%%%%%%%%%%%%%%%CO2平衡 Pc_co21=P_tc_time.*3.25;%碳捕系统单位电量捕获的cO2系数为3.25 %%环境成本 %机组排放及治理 co2=sum(sum(Pc_co21)./12); % v_co2=max([0,co2]);%%看是否全部处理，如果全部处理，则可取0 %%----------------------计算碳排放权交易量------------------------ %%清洁能源+储能发电量+碳补系统处理的量+电解装置消耗的量 %%碳排放价格按照广州碳排放交易中心确定（45/吨）（http://www.cnemission.cn/article/jydt/） %%碳排放指标参考（排放指标配额技术指南发改委）https://www.sohu.com/a/442954915_651733 % 燃气轮机：０．３９ｔ／ＭＷ;燃煤机组０．９７ｔ／ＭＷ cost2=(-co2)*45./1000+sum(pv)*0.97*45./1000+sum(PMT).*0.39*45./1000; %%Ｓｏ２　０．.８０　氮　０．８５ %燃气轮机排放及治理 % cost2=cost2+sum(PMT).*0.2*62.96;%kg/元 cost3=-sum(P_tanbu);%%弃DG量 cost3_1=max([0,cost3]); %% 运维成本 cost4=cost4+sum(pv+PW)*0.045; cost4=cost4+sum(abs(Pcd1))*0.05; cost4=cost4+sum(P_tc_time*0.025); cost4=cost4+sum(P_heat)*0.15+sum(Guol_dian*0.12)+sum(PXSL_cool)*0.2; %%输出图像 figure;%风光 plot(1:24,pv,'b-o','LineWidth',2);hold on tt=24/288:24/288:24; plot(tt,PVV,'k--o');hold off ylabel('功率／ＭＷ'); xlabel('时间／Ｔ'); legend('日前光伏','实时光伏') figure;%风光 plot(1:24,PW,'b-o','LineWidth',2);hold on tt=24/288:24/288:24; plot(tt,PW11,'k--o');hold off ylabel('功率／ＭＷ'); xlabel('时间／Ｔ'); legend('日前风电','实时风电') figure; bar([pv;PW;Pgrid;Pcd1;-P_tanbu]','stack');hold on %bar(Pcd1)% %bar(-P_tanbu) plot(1:24,pv+PW-P_tanbu+Pgrid,'k-->','LineWidth',2); ylabel('功率／ＭＷ'); xlabel('时间／Ｔ'); leg