77号资源-源程序：论文可在知网下载《基于合作博弈的综合能源系统利益分配优化调度》本人博客有解读

clc clear tic global F1 F2 F3 F12 F13 F23 F %% 基础数据 T=24; K_gas=37.72; %天然气低位燃烧热值（MJ/m^3） nita_ch4=0.4; %甲烷反应器转换效率 nita_H2=0.87; %电解槽转换效率 nita_CB=0.7; %电解槽转换效率 alpha_co2=1.76; %生成单位天然气所需CO2的量（kg） CHP_h=0.8; %CHP热效率 CHP_e=0.35; beta1=0.2; %弃风惩罚费用 beta2=0.2; %弃光惩罚费用 alpha_net=0.0046*10^(-3)/7.2; %过网费成本系数alpha beta_net=1.548*10^(-3)/7.2; %过网费成本系数beta pw=0.039; %风光运维成本 c_chp=0.021; %CHP运维成本 c_p2g=20*10^(-3)/7.2; %P2G运维费用 load_e=xlsread('data.xlsx','C2:C25')'; load_h=xlsread('data.xlsx','D2:D25')'; load_g=xlsread('data.xlsx','E2:E25')'; load_H=xlsread('data.xlsx','G2:G25')'; C_e=xlsread('data.xlsx','A2:A25')'; C_g=xlsread('data.xlsx','B2:B25')'; S_w=xlsread('data.xlsx','F2:F25')'; S_p=xlsread('data.xlsx','H2:H25')'; C=[]; %% 决策变量定义 P_wt=sdpvar(1,T); %风机消纳功率 P_pv=sdpvar(1,T); %光伏消纳功率 V_gas=sdpvar(1,T); %P2G产生天然气体积 Q_co2=sdpvar(1,T); %制气厂消耗二氧化碳的量 P_H2=sdpvar(1,T); %电解池产生的氢气量 P_P2G=sdpvar(1,T); %P2G消耗电功率量 P_h2cha=sdpvar(1,T); %储氢罐充气功率 P_h2dis=sdpvar(1,T); %储氢罐放气功率 H_h2cha=binvar(1,T); %储氢罐充气变量 H_h2dis=binvar(1,T); %储氢罐放气变量 P_H2total=sdpvar(1,T); E_h=sdpvar(1,T); %储氢罐放气功率 P_CHPe=sdpvar(1,T); %CHP发电功率 P_CHPh=sdpvar(1,T); %CHP发热功率 P_CHPg=sdpvar(1,T); %CHP耗气功率 Q_co2out=sdpvar(1,T); %CHP燃烧产生的co2 P_be=sdpvar(1,T); %购电功率 P_bg=sdpvar(1,T); %购气功率 V_CB=sdpvar(1,T); %燃气锅炉进气量 Q_CB=sdpvar(1,T); %燃气锅炉发热量 P_MRh=sdpvar(1,T); %甲烷反应器输入氢功率 P_MRg=sdpvar(1,T); %甲烷反应器输出天然气功率 %% 风光模型 C=[C,0<=P_wt<=S_w]; C=[C,0<=P_pv<=S_p]; %% P2G设备运行模型 C=[C,Q_co2==alpha_co2*nita_ch4*load_H]; C=[C,P_H2==nita_H2*P_P2G]; %电解槽产生氢气功率 C=[C,P_H2total==P_H2+P_h2dis-P_h2cha]; C=[C,0<=P_H2<=500]; %% 甲烷反应器运行模型 C=[C,P_MRg==nita_ch4*P_MRh]; %电解槽产生氢气功率 C=[C,0<=P_MRg<=250]; %% 储氢罐模型 C = [C, H_h2dis + H_h2cha <= 1]; C = [C, 0 <= P_h2cha <=H_h2cha*50]; C = [C, 0 <= P_h2dis<= H_h2dis*50]; for t = 2:24 C = [C, E_h(t) == E_h(t-1) + 0.9*P_h2cha - P_h2dis/0.9]; end C = [C, 0.1 <= E_h(:)/200 <= 0.9]; C = [C, E_h(1)==E_h(24)]; %% CHP运行约束 C=[C,P_CHPe==P_CHPg*CHP_e]; C=[C,P_CHPh==P_CHPg*CHP_h]; C=[C,Q_co2out==P_CHPg*alpha_co2]; C=[C,0<=P_CHPg<=600]; for t=1:1:23 C=[C,-50<=P_CHPg(t+1)-P_CHPg(t)<=50]; end %% 燃气锅炉约束 C=[C,Q_CB==K_gas*nita_CB*V_CB]; C=[C,0<=Q_CB<=800]; for t=1:1:23 C=[C,-100<=Q_CB(t+1)-Q_CB(t)<=100]; end %% 功率平衡约束 C=[C,0<=P_be<=1000]; C=[C,0<=P_bg<=1000]; %电功率平衡 C=[C,load_e+P_P2G==P_wt+P_pv+P_CHPe+P_be]; %热功率平衡 C=[C,load_h==P_CHPh+Q_CB]; %气功率平衡 C=[C,load_g+P_CHPg+V_CB==P_MRg+P_bg]; %氢功率平衡 C=[C,load_H+P_MRh==P_H2+P_h2dis-P_h2cha]; %% 目标函数 Fe=sum(load_e.*C_e); Fh=sum(0.65*load_h); Fg=sum(load_g.*C_g); CDE=sum(P_CHPe*c_chp+c_p2g*P_P2G+pw*(P_wt+P_pv)); CQ=sum(beta1*(S_w-P_wt)+beta2*(S_p-P_pv)); Cbe=sum(C_e.*P_be); Cbg=sum(C_g.*P_bg); Cnet=sum(alpha_net*(P_P2G).^2+beta_net*(P_P2G)+alpha_net*(load_e).^2+beta_net*(load_e)); Cco2=0.5*sum((Q_co2out-Q_co2+0.92*P_be)); Ftotal=Fe+Fh+Fg-(CDE+CQ+Cbe+Cbg+Cnet); F=0.5*Ftotal-0.5*Cco2; %% 5.设求解器 ops = sdpsettings('verbose', 1, 'solver', 'cplex'); sol = optimize(C,-F); %% if sol.problem == 0 disp('succcessful solved'); else disp('error'); yalmiperror(sol.problem) end toc % [model,recoveryalmip,diagnostic,internalmodel]=export(C,-F,ops);%转为cplex模型 % milpt=Cplex('milp for htc'); % milpt.Model.sense='minimize'; % milpt.Model.obj=model.f; % milpt.Model.lb=model.lb; % milpt.Model.ub=model.ub; % milpt.Model.A=[model.Aineq;model.Aeq]; % milpt.Model.lhs=[-inf*ones(size(model.bineq,1),1);model.beq]; % milpt.Model.rhs=[model.bineq;model.beq]; % milpt.Model.ctype=model.ctype; % milpt.writeModel('ab.lp');%输出cplex模型（注意大小写） % milpt.solve();%模型求解 P_wt=value(P_wt); %风机消纳功率 P_pv=value(P_pv); %光伏消纳功率 V_gas=value(V_gas); %P2G产生天然气体积 Q_co2=value(Q_co2); %制气厂消耗二氧化碳的量 P_H2=value(P_H2); %电解池产生的氢气量 P_P2G=value(P_P2G); %P2G消耗电功率量 P_h2cha=value(P_h2cha); %储氢罐充气功率 P_h2dis=value(P_h2dis); %储氢罐放气功率 H_h2cha=value(H_h2cha); %储氢罐充气变量 H_h2dis=value(H_h2dis); %储氢罐放气变量 P_H2total=value(P_H2total); E_h=value(E_h); %储氢罐放气功率 P_CHPe=value(P_CHPe); %CHP发电功率 P_CHPh=value(P_CHPh); %CHP发热功率 P_CHPg=value(P_CHPg); %CHP耗气功率 Q_co2out=value(Q_co2out); %CHP燃烧产生的co2 P_be=value(P_be); %购电功率 P_bg=value(P_bg); %购气功率 V_CB=value(V_CB); %燃气锅炉进气量 Q_CB=value(Q_CB); %燃气锅炉发热量 P_MRh=value(P_MRh); %甲烷反应器输入氢功率 P_MRg=value(P_MRg); %甲烷反应器输出天然气功率 F=value(F); %% 绘图 figure(1); plot(load_e,'linewidth',1.5); hold on plot(load_h,'-*','linewidth',1.5); hold on plot(load_g,'--','linewidth',1.5); hold on plot(load_H,'-^','linewidth',1.5); hold on axis([1,24,0,1000]); title('典型日负荷曲线'); legend('电负荷','热负荷','气负荷','氢负荷'); xlabel('时刻'); ylabel('功率/kW'); grid on figure(2); plot(P_wt,'r--','linewidth',1.5); hold on plot(P_pv,'b-*','linewidth',1.5); hold on axis([1,24,0,1000]); legend('风机出力','光伏出力'); xlabel('时刻'); ylabel('功率/kW'); grid on figure(3); stairs(C_e,'r--','linewidth',1.5); hold on stairs(C_g,'b-*','linewidth',1.5); hold on axis([1,24,0,4]); legend('购电电价','购气价格'); xlabel('时刻'); ylabel('功率/kW'); grid on figure(4) F11=zeros(2,24); F11(1,:)=load_e; F11(2,:)=P_P2G; AA=bar(F11(:,:)','stack'); set(AA(1),'FaceColor',[0.1 0.5 0.9]); set(AA(2),'FaceColor',[0.9 0.1 0.5]); hold on FF=zeros(4,24); FF(1,:)=P_wt; FF(2,:)=P_pv; FF(3,:)=P_CHPe; FF(4,:)=P_be; bar(-FF(:,:)','stack') axis([0.5,24.5,-1500,1500]); xlabel('时刻/h'); ylabel('电功率/kW'); legend('电负荷','P2G消耗电量','风机出力','光伏出力','CHP发电功率','购电功率'); title('电负荷平衡'); grid on figure(5) F22=zeros(3,24); F22(1,:)=load_g; F22(2,:)=P_CHPg; F22(3,:)=V_CB; BB=bar(F22(:,:)','stack'); set(BB(1),'FaceColor',[0.1 0.5 0.9]); set(BB(2),'FaceColor',[0.9 0.1 0.5]); set(BB(2),'FaceColor',[0.2 0.4 0.2]); hold on FFF=zeros(2,24); FFF(1,:)=P_MRg; FFF(2,:)=P_bg; bar(-FFF(:,:)','stack') axis([0.5,24.5,-600,600]); xlabel('时刻/h'); ylabel('气功率/kW'); legend('气负荷','CHP耗气量','燃气锅炉进气量','甲烷反应器输出天然气','购气功率'); title('气负荷平衡'); grid on %% 合作博弈处理 F1=sc1(); F2=sc2(); F3=sc3(); F12=sc12(); F13=sc13(); F23=sc23(); s=[1,0,0 0,1,0 0,0,1 1,1,0 1,0,1 0,1,1 1,1,1 ]; v=[F1 F2 F3 F12 F13 F23 F]; abcd=shapley1(s,v); clc disp(['********************']); disp(['shapley分配前的收益为：']); disp([' 新能源电厂收益为 ',num2str(F1),'元']); disp([' 燃气热电厂收益为 ',num2str(F2),'元']); disp([' 电制气厂收益为 ',num2str(F3),'元']); disp(['********************']); disp(['********************']); disp(['shapley分配后的收益为：']); disp([' 新能源电厂收益为 ',num2str(abcd(1)),'元']); disp([' 燃气热电厂收益为 ',num2str(abcd(2)),'元']); disp([' 电制气厂收益为 ',num2str(abcd(3)),'元']); disp(['********************']);