主要内容
该程序以33节点系统为例实现了考虑源-荷-储协同互动的主动配电网优化调度模型，程序采用配电网二阶锥约束、储能约束、分布式电源约束、可平移、可削减负荷约束等，以负荷调用成本、储能调用成本、储能补贴和收益、购电成本、网损成本之和作为目标函数，程序参考下图文献，虽然给出了riqian和rinei两个matlab程序，但是实际上未对日内进行滚动优化，程序采用matlab+cplex求解，程序注释清楚，方便学习！
部分代码
mpc = IEEE33BW; pload = mpc.Pload*100;%节点有功负荷,原数据100MVA，改为1MVA PPload=[1.77184023 1.760691425 1.985343246 2.164759614 2.525216573 2.836282757 2.717289933 2.912703959 2.938170633 3.200578131 3.93626314 3.680872774 3.53168293 3.215685929 2.903666471 2.626416093 2.33835944 2.582209352 3.063339069 3.465731913 3.01938202 2.591419695 2.254491338 1.784180138]; qload = mpc.Qload*100;%节点无功负荷 branch = mpc.branch; branch(:,3) = branch(:,3)*1/(10^2);%求阻抗标幺值 R = real(branch(:,3)); X = imag(branch(:,3)); T = 24;%时段数为24小时 nb = 33;%节点数 nl = 32;%支路数 nwt = 2;%风机 npv = 1;%光伏 S_pv=1.5; RR=xlsread('光照强度数据.xlsx','B2:B25'); pp_pv=xlsread('光照强度数据.xlsx','D3:AA3'); C_e=xlsread('电价.xlsx','A1:X1'); C_ae=xlsread('激励电价.xlsx','A1:X1'); C_de=xlsread('需求响应激励价格.xlsx','A1:X1'); R_STD=1000;%标准太阳辐射 R_C=150;%辐射点 PP_pv=zeros(npv,T); % %%%%光伏出力模型 for i=1:1:npv for t=1:1:T if(0<=RR(t))&(RR(t)<=R_C) PP_pv(i,t)=S_pv*((RR(t))^2)/(R_STD*R_C); else if(R_C<=RR(t))&(RR(t)<=R_STD) PP_pv(i,t)=S_pv*RR(t)/R_STD; else if(R_STD<=RR(t)) PP_pv(i,t)=S_pv; end end end end end v=xlsread('风速数据.xlsx','B2:B25'); pp_wt=xlsread('风速数据.xlsx','G2:AD2'); v_in=3;%切入风速 v_r=11.3;%额定风速 v_out=25;%切出风速 S_wt=[1;1.5]; PP_wt=zeros(nwt,T); %%%风电出力模型 for i=1:1:nwt for t=1:1:T if (0<=v(t)) &(v(t)<=v_in)|(v_out<=v(t)) PP_wt(i,t)=0; else if(v_in<=v(t))&(v(t)else if(v_r<=v(t))&(v(t)2;%ESS数 v_in=3;%切入风速 v_r=11.3;%额定风速 v_out=25;%切出风速 upstream = zeros(nb,nl); dnstream = zeros(nb,nl); for t = 1:nl upstream(t,t)=1; end for t = [1:16,18:20,22:23,25:31] dnstream(t,t+1)=1; end dnstream(1,18) = 1; dnstream(2,22) = 1; dnstream(5,25) = 1; dnstream(33,1) = 1; Vmax = [1.05*1.05*ones(nb-1,T) 1.05*1.05*ones(1,T)]; Vmin = [0.95*0.95*ones(nb-1,T) 0.95*0.95*ones(1,T)]; Pgmax = [zeros(nb-1,T) 5*ones(1,T)]; Pgmin = [zeros(nb-1,T) 0*ones(1,T)]; Qgmax = [zeros(nb-1,T) 3*ones(1,T)]; Qgmin = [zeros(nb-1,T) -1*ones(1,T)]; %% 2.设变量 V = sdpvar(nb,T);%电压的平方 I = sdpvar(nl,T);%电流的平方 P = sdpvar(nl,T);%线路有功 Q = sdpvar(nl,T);%线路无功 Pg = sdpvar(nb,T);%发电机有功 Qg = sdpvar(nb,T);%发电机无功 p_wt = sdpvar(nwt,T);%风机有功 p_pv = sdpvar(npv,T);%光伏有功 s_IL=binvar(1,T); Temp_shift=binvar(3,T); Lshift=sdpvar(3,T); Lshift_old=xlsread('可转移负荷.xlsx','A1:X3'); PLOAD=sdpvar(nb,T); Pload=sum(PLOAD); p_dch = sdpvar(ness,T);%ESS放电功率 p_ch = sdpvar(ness,T);%ESS充电功率 u_dch = binvar(ness,T);%ESS放电状态 u_ch = binvar(ness,T);%ESS充电状态 E_ess = sdpvar(ness,T);%ESS的电量 S_IL1=sdpvar(1,T,'full'); S_IL2=sdpvar(1,T,'full'); %% 3.设约束 C = []; %% 需求响应约束 %可平移负荷1约束,9-14为可平移负荷区间 C = [C,sum(Temp_shift(1,1:8))==0,sum(Temp_shift(1,9:14))==1,sum(Temp_shift(1,15:24))==0]; %可平移负荷2约束,19-22为可平移负荷区间 C = [C,sum(Temp_shift(2,1:18))==0,sum(Temp_shift(2,19:22))==1,sum(Temp_shift(2,23:24))==0]; %可平移负荷3约束,8-21为可平移负荷区间 C = [C,sum(Temp_shift(3,1:7))==0,sum(Temp_shift(3,8:21))==1,sum(Temp_shift(3,22:24))==0]; %可平移负荷1&可平移负荷2&可平移负荷3平移 for k = 1:24 C = [C, implies(Temp_shift(1,k),Lshift(1,:)==circshift(Lshift_old(1,:),[0,k-17]))]; C = [C, implies(Temp_shift(2,k),Lshift(2,:)==circshift(Lshift_old(2,:),[0,k-3]))]; C = [C, implies(Temp_shift(3,k),Lshift(3,:)==circshift(Lshift_old(3,:),[0,k-23]))]; end TL1=[zeros(3,T);Lshift_old(1,:)-Lshift(1,:);zeros(29,T)]; TL2=[zeros(14,T);Lshift_old(2,:)-Lshift(2,:);zeros(18,T)]; TL3=[zeros(29,T);Lshift_old(3,:)-Lshift(3,:);zeros(3,T)]; %%%%可削减负荷%%%%% C=[C,0<=S_IL1<=0.8*pload(10,:)]; C=[C,0<=S_IL2<=0.8*pload(26,:)]; IL1=[zeros(9,T);S_IL1(1,:);zeros(23,T)]; IL2=[zeros(25,T);S_IL2(1,:);zeros(7,T)]; C=[C,pload+TL1+TL2+TL3-IL1-IL2==PLOAD]; %% 储能装置（ESS）约束 %充放电状态约束 C = [C, u_dch(1,:) + u_ch(1,:) <= 1];%表示充电，放电，不充不放三种状态 C = [C, u_dch(2,:) + u_ch(2,:) <= 1]; C = [C, u_dch(1,:) + u_dch(2,:) <= 1]; C = [C, u_ch(1,:) + u_ch(2,:) <= 1]; %充放电时刻约束 C = [C, [sum(u_dch(1,8:15))==0,sum(u_dch(1,18:21))==0]]; C = [C, [sum(u_ch(1,1:7))==0,sum(u_ch(1,14:18))==0,sum(u_ch(1,22:24))==0]]; C = [C, [sum(u_dch(2,8:15))==0,sum(u_dch(2,18:21))==0]]; %功率约束 C = [C, 0 <= p_dch <= u_dch*0.3]; C = [C, 0 <= p_ch<= u_ch*0.3]; %容量约束 for t = 1:23 C = [C, E_ess(1,t+1) == E_ess(1,t) + 0.9*p_ch(1,t) - 0.9*p_dch(1,t)]; end for t = 1:23 C = [C, E_ess(2,t+1) == E_ess(2,t) + 0.9*p_ch(2,t) - 0.9*p_dch(2,t)]; end C = [C, 0.1<=E_ess(1,:)/1.03<=0.9]; C = [C, 0.1<=E_ess(2,:)/1.21<=0.9]; %投入节点选择 P_dch1 = [zeros(16,T);p_dch(1,:);zeros(16,T)]; P_ch1 = [zeros(16,T);p_ch(1,:);zeros(16,T)]; P_dch2 = [zeros(16,T);p_dch(2,:);zeros(16,T)]; P_ch2 = [zeros(16,T);p_ch(2,:);zeros(16,T)]; %% 风机+光伏约束 C = [C, 0 <= p_wt,p_wt <= 1.2*PP_wt]; P_wt = [zeros(15,T);p_wt(1,:);zeros(10,T);p_wt(2,:);zeros(6,T)]; C = [C, 0 <= p_pv,p_pv <=1.2*PP_pv]; p_pv =1.2*PP_pv; P_pv = [zeros(15,T);p_pv;zeros(17,T)]; %% 潮流约束 Pin = -upstream*P + upstream*(I.*(R*ones(1,T))) + dnstream*P;%节点注入有功 Qin = -upstream*Q + upstream*(I.*(X*ones(1,T))) + dnstream*Q;%节点注入无功 C = [C, Pin+pload+TL1+TL2+TL3-IL1-IL2- P_wt- P_pv- Pg-P_dch1+ P_ch1-P_dch2+ P_ch2==0]; C = [C, Qin + qload - Qg == 0]; %欧姆定律约束 C = [C, V(branch(:,2),:) == V(branch(:,1),:) - 2*(R*ones(1,24)).*P - 2*(X*ones(1,24)).*Q + ((R.^2 + X.^2)*ones(1,24)).*I]; %二阶锥约束 C = [C, V(branch(:,1),:).*I >= P.^2 + Q.^2]; %% 通用约束 %节点电压约束 C = [C, Vmin <= V,V <= Vmax]; %发电机功率约束 C = [C, Pgmin <= Pg,Pg <= Pgmax,Qgmin <= Qg,Qg <= Qgmax]; %支路电流约束 C = [C, 0 <= I,I <= 12]; %% 4.设目标函数 Pload2=PPload+sum(p_dch-p_ch); Pload3=Pload+sum(p_dch-p_ch); %%%%%%%%% for t=1:1:23 f1_0=0.25*[(max(PPload)-min(PPload))/mean(PPload(:))]+0.75*[(max(abs(PPload(t)-Pload(t+1))))/mean(PPload(:))]; f1_1=0.25*[(max(Pload)-min(Pload))/mean(Pload(:))]+0.75*[(max(abs(Pload(t)-Pload(t+1))))/mean(Pload(:))]; f1_2=0.25*[(max(Pload2)-min(Pload2))/mean(Pload2)]+0.75*[(max(abs(Pload2(t)-Pload2(t+1))))/mean(Pload2(:))]; f1_3=0.25*[(max(Pload3)-min(Pload3))/mean(Pload3)]+0.75*[(max(abs(Pload3(t)-Pload3(t+1))))/mean(Pload3(:))]; end %%%%%%%电压差最小%%%%%%%%% f1=0; for t=1:1:T for n=1:1:nb-1 f1=f1+abs(V(n,t)-V(n+1,t)); end end %%%%%%%%% B1=sum(sum(p_dch-p_ch).*C_ae); B2=sum(sum(p_dch))*200; R_ess=B1+B2; %储能补贴和收益%%%%% C_LA=sum([0.062,0.043,0.048]*Temp_shift.*C_de)+sum((S_IL1+S_IL2).*C_de);%%%负荷调用成本%%%% C_buy=sum(sum(Pg).*C_e); C_loss=sum(sum(I.*(R*ones(1,T))))*400; %%%网损成本%%%% C_ess=250*sum(sum(p_dch)+sum(p_ch)); %%%储能调用成本%%%% f2=C_LA+C_loss+C_ess-(B1+B2)+C_buy; %%%%储能寿命损耗成本%%%% n=sum(sum(u_dch + u_ch)); L=0; for D=0:0.1:1 L=L+n/(71470*D^4-170100*D^3+146400*D^2-56500*D+12230); end C_ESS1=819000*2+2953000*0.375; C_ESS2=0; for t=1:1:10 C_ESS2=C_ESS2+69000*2*((1+0.015)/(1+0.09))^t; end C_ESS=C_ESS1+C_ESS2; C_day=C_ESS*L; %%%%%储能寿命损耗成本%%%%%% %%%%%%%%%%%%%%% %%%%%%配电网综合收益%%%%%% peak_0=PPload(9)+PPload(10)+PPload(11)+PPload(12)+PPload(13)+PPload(14)+PPload(19)+PPload(20); peak_3=Pload3(9)+Pload3(10)+Pload3(11)+Pload3(12)+Pload3(13)+Pload3(14)+Pload3(19)+Pload3(20); lamda3=(peak_0-peak_3)/peak_0; %%%%场景三的削峰率 C_gridup=1; delta_n3=((log10(1+lamda3))/log10(1+0.015)); %%%%%%延缓改造年限%%%% F1=C_gridup*[1-((1+0.015)/(1+0.09))^delta_n3]; %%%%%%减少配电网升级改造费用%%%%
结果一览
下载链接
