90号资源-源程序：论文可在知网下载《计及需求响应的区域综合能源系统双层优化调度策略》本人博客有解读

%======================================================================== %中国电力期刊论文 %论文题目---《计及需求响应的区域综合能源系统双层优化调度策略》 %双层优化：KKT条件（双层模型转化为单层模型）+ Big-M法（非线性模型转化为混合整数线性模型） %======================================================================== clc; clear; close all; warning off; tic %% 选择用哪个求解器求解 %% 常量定义 WT=xlsread('data','B2:B25'); %风机出力数据导入 Ee0=xlsread('data','C2:C25'); %原始电负荷数据导入 Eh1=xlsread('data','D2:D25'); %热负荷数据导入 WT=WT'; Ee0=Ee0'; Eh1=Eh1'; M=1E8; %Big-M法中的M Pgbuy=1.2*ones(1,24); %RIESO向电网购电电价 Pgsell=0.5*ones(1,24); %RIESO向电网售电电价 Phbuy=0.3*ones(1,24); %RIESO向热网购热价格 Phsell=0.2*ones(1,24); %RIESO向热网售热价格 Pesell=[0.4*ones(1,6),0.7*ones(1,2),1.2*ones(1,3),0.7*ones(1,5),... 1.2*ones(1,4),0.7*ones(1,2),0.4*ones(1,2)]; %RIESO向用户的售电电价 Phsell_u=0.3*ones(1,24); %RIESO向用户的售热价格 CLS=0.55*ones(1,24); %向DRA购买转移负荷的价格 CLC=0.55*ones(1,24); %向DRA购买削减负荷的价格 CLScus=0.3*ones(1,24); %向用户购买单位负荷转移量成本 CLCcus=0.5*ones(1,24); %向用户购买单位负荷削减量的补偿成本 %% 决策变量定义 ELSdown=sdpvar(1,24); %电负荷需求响应转出量 ELSup=sdpvar(1,24); %电负荷需求响应转入量 ELC=sdpvar(1,24); %电负荷需求响应削减量 Ee=sdpvar(1,24); %RIESO向电网的购/售电量 Ee_b=sdpvar(1,24); %RIESO向电网的购电量 Ee_s=sdpvar(1,24); %RIESO向电网的售电量 Eh=sdpvar(1,24); %RIESO向热网的购/售热量 Eh_b=sdpvar(1,24); %RIESO向热网的购热量 Eh_s=sdpvar(1,24); %RIESO向热网的售热量 EICE=sdpvar(1,24); %内燃发电机的输出功率 EGB=sdpvar(1,24); %燃气锅炉的输出功率 ERE=sdpvar(1,24); %余热锅炉的输出功率 Ee1=sdpvar(1,24); %用户参与需求响应后的电负荷 Temp_b_e=binvar(1,24); %RIESO向电网买电标识 Temp_s_e=binvar(1,24); %RIESO向电网售电标识 Temp_static_e=binvar(1,24); %RIESO购/售电静置标识 Temp_b_h=binvar(1,24); %RIESO向热网购热标识 Temp_s_h=binvar(1,24); %RIESO向热网售热标识 Temp_static_h=binvar(1,24); %RIESO购/售热静置标识 %% 上层RIESO的约束条件 C= []; for t=1:24 %电功率平衡约束 C=[C,Ee0(1,t)-ELSdown(1,t)+ELSup(1,t)-ELC(1,t)==WT(1,t)+EICE(1,t)+Ee(1,t)]; %热功率平衡约束 C=[C,Eh1(1,t)==((EICE(1,t)*0.65*0.83/0.35)*0.85+EGB(1,t)*0.85)+Eh(1,t)]; %需求响应后的电负荷 C=[C,Ee1(1,t)==Ee0(1,t)-ELSdown(1,t)+ELSup(1,t)-ELC(1,t)]; %余热锅炉的输出功率 C=[C,ERE(1,t)==EICE(1,t)*0.65*0.83/0.35]; %机组出力上下限约束 C=[C,0<=EICE(1,t)<=500,0<=EGB(1,t)<=800]; %与电网热网能量交互约束 C=[C,-1000<=Ee(1,t)<=1000,0<=Ee_b(1,t)<=1000,-1000<=Ee_s(1,t)<=0]; C=[C,implies(Temp_b_e(1,t),[Ee(1,t)>=0,Ee_b(1,t)==Ee(1,t),Ee_s(1,t)==0])]; C=[C,implies(Temp_s_e(1,t),[Ee(1,t)<=0,Ee_b(1,t)==0,Ee_s(1,t)==Ee(1,t)])]; C=[C,implies(Temp_static_e(1,t),[Ee(1,t)==0,Ee_b(1,t)==0,Ee_s(1,t)==0])]; C=[C,Temp_b_e(1,t)+Temp_s_e(1,t)+Temp_static_e(1,t)==1]; C=[C,-100<=Eh(1,t)<=100,0<=Eh_b(1,t)<=100,-100<=Eh_s(1,t)<=0]; C=[C,implies(Temp_b_h(1,t),[Eh(1,t)>=0,Eh_b(1,t)==Eh(1,t),Eh_s(1,t)==0])]; C=[C,implies(Temp_s_h(1,t),[Eh(1,t)<=0,Eh_b(1,t)==0,Eh_s(1,t)==Eh(1,t)])]; C=[C,implies(Temp_static_h(1,t),[Eh(1,t)==0,Eh_b(1,t)==0,Eh_s(1,t)==0])]; C=[C,Temp_b_h(1,t)+Temp_s_h(1,t)+Temp_static_h(1,t)==1]; end for t=2:24 C=[C,-100<=EICE(1,t)-EICE(1,t-1)<=100];%内燃机的上升和下降爬坡速率约束 C=[C,-300<=EGB(1,t)-EGB(1,t-1)<=300]; %燃气锅炉的上升和下降爬坡速率约束 end %% 进行下层的KKT条件处理 %下层等式约束 C=[C,sum(ELSdown)==sum(ELSup)]; %处理式（29）-（31）KKT条件 %定义等式约束的拉格朗日乘子 Delta1=sdpvar(1,24); Delta2=sdpvar(1,24); %定义不等式约束的拉格朗日乘子 u1bt=sdpvar(1,24); u2bt=sdpvar(1,24); u3at=sdpvar(1,24); u4at=sdpvar(1,24); u5t=sdpvar(1,24); u6ct=sdpvar(1,24); u7ct=sdpvar(1,24); Alpha=binvar(1,24); %时刻t第a段电负荷转出的状态变量 Beta=binvar(1,24); %时刻t第b段电负荷转入的状态变量 Chi=binvar(1,24); %时刻t第c段电负荷削减的状态变量 %依次添加式（29）-（31）KKT平衡条件约束 for t=1:24 C=[C, CLS(1,t)-CLScus(1,t)-u3at(1,t)+u4at(1,t)+Delta1(1,t)==0, -u1bt(1,t)+u2bt(1,t)-Delta1(1,t)==0, CLC(1,t)-CLCcus(1,t)-u6ct(1,t)+u7ct(1,t)==0, ]; end %处理式（39）-（45）KKT条件 %非线性不等式约束转换为线性不等式约束 %引入Big-M法所需的布尔变量 0-1变量 v39=binvar(1,24); v40=binvar(1,24); v41=binvar(1,24); v42=binvar(1,24); v43=binvar(1,24); v44=binvar(1,24); v45=binvar(1,24); %Big-M法处理割平面约束 for t=1:24 C=[C, 0<=ELSup(1,t)-0.1*Ee0(1,t)*Beta(1,t)<=M*v39(1,t), 0<=u1bt(1,t)<=M*(1-v39(1,t)), 0<=1*Ee0(1,t)*Beta(1,t)-ELSup(1,t)<=M*v40(1,t), 0<=u2bt(1,t)<=M*(1-v40(1,t)), 0<=ELSdown(1,t)-0.1*Ee0(1,t)*Alpha(1,t)<=M*v41(1,t), 0<=u3at(1,t)<=M*(1-v41(1,t)), 0<=1*Ee0(1,t)*Alpha(1,t)-ELSdown(1,t)<=M*v42(1,t), 0<=u4at(1,t)<=M*(1-v42(1,t)), 0<=1-Alpha(1,t)-Beta(1,t)<=M*v43(1,t), 0<=u5t(1,t)<=M*(1-v43(1,t)), 0<=ELC(1,t)-0.1*Ee0(1,t)*Chi(1,t)<=M*v44(1,t), 0<=u6ct(1,t)<=M*(1-v44(1,t)), 0<=1*Ee0(1,t)*Chi(1,t)-ELC(1,t)<=M*v45(1,t), 0<=u7ct(1,t)<=M*(1-v45(1,t)), ]; end %% 上层RIESO的目标函数 %售能收益 Rsell=sum((Ee0(1,:)-ELSdown(1,:)+ELSup(1,:)-ELC(1,:)).*Pesell+Eh1(1,:).*Phsell_u); %DR资源的购买成本 Cdr=sum(ELSdown(1,:).*CLS+ELC(1,:).*CLC); %购电成本 Cgrid=sum(Ee_b(1,:).*Pgbuy+Ee_s(1,:).*Pgsell); %购热成本 Chbuy=sum(Eh_b(1,:).*Phbuy+Eh_s(1,:).*Phsell); %燃料成本 Cfuel=sum((0.0013*EICE(1,:).*EICE(1,:)+0.16*EICE(1,:))+(0.0005*EGB(1,:).*EGB(1,:)+0.11*EGB(1,:))); %设备运维成本 Cope=sum(0.03*WT(1,:)+0.02*EICE(1,:)+0.025*ERE(1,:)+0.02*EGB(1,:)); %目标函数 objective=Rsell-Cdr-Cgrid-Chbuy-Cfuel-Cope; %上层优化目标为RIESO的日净利润最大 %% 下层DRA的目标函数 OF_lower=Cdr-sum(CLScus(1,:)*0.3)-sum(CLCcus(1,:)*0.5); %下层优化目标为DRA的日净利润最大 %% P_CLS=sdpvar(1,24); P_CLC=sdpvar(1,24); PCdr=sdpvar(1,24); OF_upper=sdpvar(1,24); for t=1:24 P_CLS(1,t)=ELSdown(1,t)*CLS(1,t); %DRA转移利润（非纯利润） P_CLC(1,t)=ELC(1,t)*CLC(1,t); %DRA削减利润（非纯利润） PCdr(1,t)=ELSdown(1,t)*CLS(1,t)+ELC(1,t)*CLC(1,t);%DRA利润 OF_upper(1,t)=((Ee0(1,t)-ELSdown(1,t)+ELSup(1,t)-ELC(1,t))*Pesell(1,t)+Eh1(1,t)*Phsell_u(1,t))-... PCdr(1,t)-(Ee_b(1,t)*Pgbuy(1,t)+Ee_s(1,t)*Pgsell(1,t))-(Eh_b(1,t)*Phbuy(1,t)+Eh_s(1,t)*Phsell(1,t))-... ((0.0013*EICE(1,t)*EICE(1,t)+0.16*EICE(1,t))+(0.0005*EGB(1,t)*EGB(1,t)+0.11*EGB(1,t)))-... (0.03*WT(1,t)+0.02*EICE(1,t)+0.025*ERE(1,t)+0.02*EGB(1,t)); end ops = sdpsettings('solver','cplex', 'verbose', 2);%参数指定程序用cplex求解器 result1=optimize(C,-objective,ops) if result1.problem == 0 disp('求解成功'); else disp('求解过程中出错'); end %转化字符类型，输出运行结果 value(ELSdown); value(ELSup); value(ELC); value(Ee); value(EICE); EGB=0.85*value(EGB); ERE=0.85*value(ERE); value(Ee1); value(Eh1); Eh=value(Eh); P_CLS=value(P_CLS); P_CLC=value(P_CLC); Ee_s=value(Ee_s); objective=value(objective); %RIESO的最大日净利润 disp(['RIESO的最大日净利润为: ', num2str(objective),'元']); OF_lower=value(