case.zip_天然气_天然气网络_气 潮流_潮流_潮流 耦合

%The day ahead unit commitment reserve problem clear all E_D=[175.19 165.15 158.67 154.73 155.06 160.48 173.39 190.40 205.56 217.2 228.61 236.10 242.18 243.60 248.86 255.79 256.00 246.74 245.97 237.35 237.31 227.14 201.05 196.75]; %E_W=[16.39,10.53,8.69,8.29,17.92,16.85,17.64,13.34,18.2,17.84,15.86,19.64,14.93,11.76,9.24,10.14,11.08,12.7,16.45,18.62,16.08,18.04,19.05,17.68]; E_W=zeros(1,24); D_a=E_D*0.2 ; %Adjustable load=emission coefficient*Demand=0.2*Demand delta_d=256*0.03; delta_w=54*0.15; eta=0.95; %置信度 etal=0.9999; No_generators=3; No_timesteps=24; mpc=case6; a=mpc.gencost(:,5)*ones(1,No_timesteps); %成本函数中的 a*P*P b=mpc.gencost(:,6)*ones(1,No_timesteps); % b*P c=mpc.gencost(:,7)*ones(1,No_timesteps); % c Pmax=mpc.gen(:,9)*ones(1,No_timesteps); Pmin=mpc.gen(:,10)*ones(1,No_timesteps); RU=[55 40 5]'*ones(1,No_timesteps); RD=[55 40 5]'*ones(1,No_timesteps); P_indices=mpc.gen(:,1); No_buses=size(mpc.bus,1); No_branches=size(mpc.branch,1); PLmax=mpc.branch(:,6)*ones(1,No_timesteps); %{ i=mpc.load(:,2); j=mpc.load(:,3); c_ij=mpc.load(:,4); p_min=mpc.gas(:,2); p_max=mpc.gas(:,3); %} %Q_sp14_max=5000; %Q_sp14_min=1500; %Q_sp26_max=6000; %Q_sp26_min=2000; pf=[2.25 2.25 2.25]'*ones(1,No_timesteps); csu=[100 200 0]'*ones(1,No_timesteps); %Start up cost cdu=[0 0 0]'*ones(1,No_timesteps); %Shut down cost pWC=2; [Ybus, Bf, Pbusinj, Pfinj] = makeBdc(mpc.baseMVA, mpc.bus, mpc.branch); inv_B=[[inv(Ybus(1:No_buses-1,1:No_buses-1)) zeros(No_buses-1,1)];zeros(1,No_buses)]; for i=1:No_branches BL(i,mpc.branch(i,1))=-Ybus(mpc.branch(i,1),mpc.branch(i,2)); BL(i,mpc.branch(i,2))=Ybus(mpc.branch(i,1),mpc.branch(i,2)); end Q=BL*inv_B; P=sdpvar(No_generators,No_timesteps); %Power output of the thermal gererators P0=[150 0 0]'; WC=sdpvar(1,No_timesteps); %Wind curtailment E_W(1,:)=1*E_W; x=binvar(No_generators,No_timesteps); %"1"if the gererator starts up y=binvar(No_generators,No_timesteps); %"1"if the gererator shuts down u=binvar(No_generators,No_timesteps); %"1"if the gererator is on u0=[1 0 0]'; w_indices=5; Q_G=Q(:,P_indices); Q_D=Q; Q_W=Q(:,w_indices); D=zeros(No_buses,No_timesteps); D(1,:)=E_D/3; D(2,:)=E_D/3; D(6,:)=E_D/3; No_Gasnodes=size(mpc.gasnodes,1); %天然气节点数量 No_Gaspipeline=size(mpc.gaspipeline,1); %天然气管道数量 Gaspress=sdpvar(No_Gasnodes,No_timesteps); %天然气压力 Gaspress_max=mpc.gasnodes(:,3); Pipelineflow=sdpvar(No_Gaspipeline,No_timesteps); %天然气管道流量 Pipelineflow_fnode_index=mpc.gaspipeline(:,2); Pipelineflow_tnode_index=mpc.gaspipeline(:,3); No_Gassupplier=size(mpc.gassupplier,1); %天然气源节点数量 Gassupplier_indices=mpc.gassupplier(:,2); %天然气源节点编号 Gassupplier=sdpvar(No_Gassupplier,No_timesteps); %天然气源 NS=zeros(No_Gasnodes,No_Gassupplier); for i=1:No_Gassupplier NS(Gassupplier_indices(i),i)=1; %NS为系数矩阵，通过NS*Gassupplier得到与Gaspress维数相同的矩阵; end No_Gasloadfixed=size(mpc.gasloadfixed,1); %天然气负荷数量 Gasloadfixed_indices=mpc.gasloadfixed(:,2); %天然气负荷节点编号 Gasloadfixed=mpc.gasloadfixed(:,3)*ones(1,No_timesteps); %天然负荷 NFL=zeros(No_Gasnodes,No_Gasloadfixed); for i=1:No_Gasloadfixed NFL(Gasloadfixed_indices(i),i)=1; %NFL为系数矩阵，通过NFL*Gasloadfixed得到与Gaspress维数相同的矩阵; end No_Gasloadunfixed=size(mpc.gasloadunfixed,1); %燃气轮机数量 Gasloadunfixed_indices=mpc.gasloadunfixed(:,2); %燃气轮机节点编号 Gasloadunfixed=sdpvar(No_Gasloadunfixed,No_timesteps); %燃气轮机耗气量（待决策的气负荷） NUFL=zeros(No_Gasnodes,No_Gasloadunfixed); for i=1:No_Gasloadunfixed NUFL(Gasloadunfixed_indices(i),i)=1; %NUFL为系数矩阵，通过NUFL*Gasloadunfixed得到与Gaspress维数相同的矩阵; end V=zeros(No_Gasnodes,No_Gaspipeline); for i=1:No_Gasnodes for j=1:No_Gaspipeline if Pipelineflow_fnode_index(j)==i V(i,j)=1; elseif Pipelineflow_tnode_index(j)==i V(i,j)=-1; %V为系数矩阵，通过V*Pipelineflow得到各节点注入功率 end end end %约束条件 F=set(sum(P)+E_W-WC==E_D); %功率平衡约束 F=F+set(Q_G*P+Q_W*E_W-Q_D*D<=PLmax); %潮流约束 F=F+set(u.*Pmin<=P<=u.*Pmax); F=F+set(x-y==u-[u0 u(:,1:No_timesteps-1)]); F=F+set(x+y<=1); F=F+set(-RD<=P-[P0 P(:,1:No_timesteps-1)]<=RU); F=F+set(WC>=0); %F=F+set(Pipelineflow==Gaspress(Pipelineflow_fnode_index,:)-Gaspress(Pipelineflow_tnode_index,:)); %各个管道流量与各个节点压力的关系 %F=F+set(0<=Gaspress<=Gaspress_max*ones(1,No_timesteps)); %各个节点压力不能越限 %F=F+set(0<=Gaspress); F=F+set(NS*Gassupplier-NFL*Gasloadfixed-NUFL*Gasloadunfixed==V*Pipelineflow); %各节点功率守恒(类似于基尔霍夫电流定律，注入功率等于与其相邻各支路功率之和) F=F+set(Gasloadunfixed==(c+b.*P)); %F=F+set(Gasloadunfixed==(c+b.*P+a.*P.^2)); %电功率与气负荷的关系 %F_P=(c+b.*P+a.*P.^2); %SU=pf.*csu.*x; %SD=pf.*cdu.*y; %Cost_func=sum(sum(F_P+SU+SD))+sum(pWC*WC); Cost_func=sum(sum(Gassupplier)); %ops = sdpsettings('solver','cplex');sol = solvesdp(F,Cost_func,ops); sol = solvesdp(F,Cost_func); u_double=double(u); P_double=double(P); WC_double=double(WC); t=[1:24] figure; plot(t,double(sum(P)),'bo-',t,double(P(1,:)),'r^-',t,double(P(2,:)),'g^-',t,double(P(3,:)),'m^-'); title('燃气机调度结果'); legend('燃气机总出力','第一台燃气机出力','第二台燃气机出力','第三台燃气机出力'); figure; plot(t,double(Gassupplier(1,:)),t,double(Gassupplier(2,:))); legend('气源1','气源2'); Total_cost=double(Cost_func); P_line=double(Q_G*P+Q_W*E_W-Q_D*D); %plot(P_line(2,:)); %{ t=[1:24] figure; plot(t,E_D,'ko-',t,double(sum(P)),'bo-',t,E_W,'go-',t,double(DR),'ro-',t,E_D-double(DR),'mo-',t,zeros(1,24),'k'); title('不同类型机组出力情况'); legend('预测负荷','燃气机出力','风电机组出力','可转移负荷','预测负荷-可转移负荷'); figure; plot(t,double(sum(P)),'bo-',t,double(P(1,:)),'r^-',t,double(P(2,:)),'g^-',t,double(P(3,:)),'m^-'); title('燃气机组合情况'); legend('燃气机总出力','第一台燃气机出力','第二台燃气机出力','第三台燃气机出力'); figure; plot(t,double(sum(R))+double(RDR),'bs-',t,double(sum(R)),'gs-',t,double(RDR),'rs-',t,D_a-double(DR),'ms-'); title('备用容量分布情况'); legend('所需要的总备用','燃气机提供的备用','需求相应提供的备用','可作为备用容量的需求相应负荷'); %} %{ figure plot(t,E_D-double(DR),t,E_D,t,double(sum(P)),t,E_W,t,E_W-double(WC),t,double(P(1,:)),t,double(P(2,:)),t,double(P(3,:))); figure plot(double(WC)); figure plot(t,double(sum(R))+double(RDR),t,double(sum(R)),t,double(RDR)) figure plot(t,double(RDR)+double(DR),t,D_a) %}