027配电网网架规划及数据.zip

% @Company ： % @Version : Matlab2018a % @Author : % @Contact : % @Time : % @Desc ： close all; clear all; clc; warning off; % 去除警告 tic; % tic用来保存当前时间，而后使用toc来记录程序完成时间 %% 基本参数 T=12; % 典型日 8-19h % 8-19h 负荷各时段负荷总量 total_P_LOAD=[828,1001,1105,1105,994,1105,1105,1049,1012,810,699,626]; % 8-19h 光伏各时段出力标幺(by)值 by_P_PVG_timeline=[0.29,0.03,0.51,0.55,0.47,0.46,0.55,0.64,0.22,0.20,0.38,0.03]; Node=23; % 节点 Line=34; % 线路 o1 = 1.1198; % 负荷功率因数角 0.9 o2 = 1.2532; % DG功率因素角 0.95 % 网损成本 Line_closs=0.6; % 最大弃(abandon,a)光率 a_max_PVG=0.05; % 单位弃光成本 ca_PVG=0.6; % 主网(zw)购电分时单价 zw_buy1_TR=[0.6 0.57 0.45 0.43 0.43 0.58 0.65 0.67 0.68 0.64 0.69 0.63 0.63 0.63 0.62 0.61 0.615 0.635 0.63 0.63 0.615 0.615 0.59 0.505 ]; zw_buy_TR=zw_buy1_TR(8:19); SLmax=5300; % KVA Umin=0.95*0.95*(12.66)^2; % 20KV Umax=1.05*1.05*(12.66)^2; % 20KV % P_TRmax=7500; P_TRmax=4500; P_TRmin=0; M=99999; % 大M法处理非线性项 %% 线路节点相关信息 % 从excel中读取 % Line_cs中各项含义 1线路编号 2起、3止节点 4不允许建线路 5允许建线路 6线路长度 km Line_cs=xlsread('Lines_Nodes_canshu.xlsx','B5:H38'); %% 为节点（jd）关联(gl)矩阵 % 若不为0，则单元格数值表示为线路编号，横、纵坐标分别表示为起、止节点。 JD_gl=zeros(Node); for i=1:Line JD_gl(Line_cs(i,2),Line_cs(i,3))=i; end %% 如果线路不允许建，线路编号加入到E0当中，允许建则加入到E1当中 E0=[]; % 不允许建线路集合 E1=[]; % 允许建线路集合 for i=1:Line if Line_cs(i,4)==1 E0=[E0,Line_cs(i,1)]; end end for i=1:Line if Line_cs(i,5)==1 E1=[E1,Line_cs(i,1)]; end end Num_Line=length(E1); x_E1_line=binvar(1,Num_Line); % 待建线路建设集合 %% 每条线路长度 % A_cs中不仅仅包含了线路关联矩阵，而且保留有线路长度 横、纵坐标分别表示为起、止节点 Long_line_cs=zeros(Node,Node) for i=1:Node for j=1:Node if JD_gl(i,j)~=0 Long_line_cs(i,j)=Line_cs(JD_gl(i,j),6) % 线路长度 end end end %% 线路投资成本矩阵 长度 * 单位造价 % 线路建设成本 20万元/km cline=Long_line_cs*200000; Cline=[]; for i=1:Node for j=1:Node if JD_gl(i,j)~=0 Cline=[Cline,cline(i,j)]; end end end %% 电阻电抗化简 % 电阻 rline=Long_line_cs*0.17; % 电抗 xline=Long_line_cs*0.402; Rline=[]; %%%%%%%%线路阻抗化简 Xline=[]; for i=1:Node for j=1:Node if JD_gl(i,j)~=0 Rline=[Rline,rline(i,j)]; Xline=[Xline,xline(i,j)]; end end end %% 主网、光电接入(jr)节点 % 光伏接入节点 PVG_jr=[4,8,11,12,16,20]; % 各光伏节点接入容量 int_pvg=[435,465,345,489,564,349]; Num_PVG=length(PVG_jr); P_PVG_yc=sdpvar(Num_PVG,T,'full'); % 光伏预测（yc）出力 P_PVG=sdpvar(Num_PVG,T,'full'); x_pvg=binvar(1,Num_PVG); % 0-1变量，是否建设光伏 TR_jr=[1]; % 变压器接入节点 也就是与主网相连节点 Num_TR=length(TR_jr); P_TR=sdpvar(Num_TR,T,'full'); Q_TR=sdpvar(Num_TR,T,'full'); %% 负荷相关定义 % 各个节点基准(jz)负荷 P_load_jz=[0,188,180,136,184,160,172,164,244,252,180,204,248,160,196,144,172,188,223,145,135,268,193]/4136; % 负荷曲线 0-24 求8:19 Load_timeine1=[0.78,0.75,0.7,0.68,0.65,0.63,0.7,0.75,0.79,0.90,0.90,0.90,0.90,0.9,0.9,0.9,0.9,0.9,0.9,0.9,0.85,0.83,0.80,0.78]; Load_timeine=Load_timeine1(8:19); for t=1:T for j=1:Node P_LOAD(j,t)=total_P_LOAD(t)*P_load_jz(j)*Load_timeine(t); end end %% 光伏预测出力 % for t=1:T % for j=1:Node % if ismember(j,PVG_jr)==1 % for k=1:Num_PVG % P_PVG_yc(find(PVG_jr==j),t)=int_pvg(k)*by_P_PVG_timeline(t); % end % end % end % end for t=1:T for j=1:Num_PVG P_PVG_yc(j,t)=int_pvg(j)*by_P_PVG_timeline(t); end end %% 设决策变量 P_ij=sdpvar(Line,T); Q_ij=sdpvar(Line,T); I_ij=sdpvar(Line,T); % 电流的平方 U_j=sdpvar(Node,T); % 电压的平方 %% 约束 C=[]; % 总约束 %% 建线路 、光电 数量约束 C_1=[]; C_1=[C_1,sum(x_E1_line)==34]; % 各光伏接入节点均可正常运行 C_1=[C_1,sum(x_pvg)==Num_PVG]; C=[C,C_1]; %% sd：首节点，md：末节点 SJD=zeros(Node,5); MJD=zeros(Node,5); for j=1:Node sjdmjd=1; for i=1:Node if JD_gl(j,i)~=0 SJD(j,sjdmjd)=i; sjdmjd=sjdmjd+1; end end end for j=1:Node sjdmjd=1; for i=1:Node if JD_gl(i,j)~=0 MJD(j,sjdmjd)=i; sjdmjd=sjdmjd+1; end end end %% 有功无功集合 RIGHT1 有功 RIGHT2 无功 C_2=[]; LINE1=zeros(1,Node); % 以该节点为起始节点线路有几条 LINE2=zeros(1,Node); % 以该节点为终止节点线路有几条 for j=1:Node % 以该节点为起始节点线路有几条 LINE1(j)=sum(SJD(j,:)~=0,2); % 以该节点为终止节点线路有几条 LINE2(j)=sum(MJD(j,:)~=0,2); end for t=1:T for j=1:Node P_ALL=0; Q_ALL=0; P_ALL=-P_LOAD(j,t); Q_ALL=-P_LOAD(j,t)/tan(o1); if ismember(j,PVG_jr)==1 P_ALL=P_ALL+P_PVG(find(PVG_jr==j),t); Q_ALL=Q_ALL+P_PVG(PVG_jr==j,t)/tan(o2); end if ismember(j,TR_jr)==1 P_ALL=P_ALL+P_TR(find(TR_jr==j),t); Q_ALL=Q_ALL+Q_TR(find(TR_jr==j),t); end if LINE1(j)==0 if LINE2(j) ==0 % 没有建设线路 C_2=[C_2,P_ALL==0]; C_2=[C_2,Q_ALL==0]; else % 首节点功率分布 C_2=[C_2,-sum(P_ij(JD_gl(MJD(j,1:LINE2(j)),j),t))+sum(I_ij(JD_gl(MJD(j,1:LINE2(j)),j),t).*rline(MJD(j,1:LINE2(j)),j))/1000==P_ALL]; % KVA为单位进行换算 C_2=[C_2,-sum(Q_ij(JD_gl(MJD(j,1:LINE2(j)),j),t))+sum(I_ij(JD_gl(MJD(j,1:LINE2(j)),j),t).*xline(MJD(j,1:LINE2(j)),j))/1000==Q_ALL]; % KVA为单位进行换算 end else if LINE2(j) ==0 % 末节点功率分布 C_2=[C_2,sum(P_ij(JD_gl(j,SJD(j,1:LINE1(j))),t))==P_ALL]; % KVA为单位进行换算 C_2=[C_2,sum(Q_ij(JD_gl(j,SJD(j,1:LINE1(j))),t))==Q_ALL]; % KVA为单位进行换算 else % 正常线路功率分布 C_2=[C_2,-sum(P_ij(JD_gl(MJD(j,1:LINE2(j)),j),t))+sum(I_ij(JD_gl(MJD(j,1:LINE2(j)),j),t).*rline(MJD(j,1:LINE2(j)),j))/1000+sum(P_ij(JD_gl(j,SJD(j,1:LINE1(j))),t))==P_ALL]; % KVA为单位进行换算 C_2=[C_2,-sum(Q_ij(JD_gl(MJD(j,1:LINE2(j)),j),t))+sum(I_ij(JD_gl(MJD(j,1:LINE2(j)),j),t).*xline(MJD(j,1:LINE2(j)),j))/1000+sum(Q_ij(JD_gl(j,SJD(j,1:LINE1(j))),t))==Q_ALL]; % KVA为单位进行换算 end end end end C=[C,C_2]; %% 线路电压平衡约束 C_3=[]; for t=1:T for j=1:Node for i=1:Node if (JD_gl(i,j)~=0) if(ismember(JD_gl(i,j),E0)==1) C_3=[C_3,U_j(j,t)*1000==U_j(i,t)*1000-2*(P_ij(JD_gl(i,j),t)*rline(i,j)+Q_ij(JD_gl(i,j),t)*xline(i,j))+I_ij(JD_gl(i,j),t)*(rline(i,j)^2+xline(i,j)^2)/1000]; else C_3=[C_3,U_j(i,t)*1000-U_j(j,t)*1000-2*(P_ij(JD_gl(i,j),t)*rline(i,j)+Q_ij(JD_gl(i,j),t)*xline(i,j))+I_ij(JD_gl(i,j),t)*(rline(i,j)^2+xline(i,j)^2)/1000<= M*(1-x_E1_line(find(E1==JD_gl(i,j))))] C_3=[C_3,U_j(i,t)*1000-U_j(j,t)*1000-2*(P_ij(JD_gl(i,j),t)*rline(i,j)+Q_ij(JD_gl(i,j),t)*xline(i,j))+I_ij(JD_gl(i,j),t)*(rline(i,j)^2+xline(i,j)^2)/1000>=-M*(1-x_E1_line(find(E1==JD_gl(i,j))))] end end end end end C=[C,C_3]; %% 二阶锥约束 C_4=[]; fo