电力系统网架规划MATLAB程序分享

资源地址：

电力系统网架规划MATLAB程序_网架规划matlab程序-电子商务文档类资源-CSDN文库

网架数据展示：

完整程序：

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

for i=1:Line

if Line_cs(i,5)==1

E1=[E1,Line_cs(i,1)];

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

%% 线路投资成本矩阵长度 * 单位造价

% 线路建设成本 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

%% 电阻电抗化简

% 电阻

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

%% 主网、光电接入(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

%% 光伏预测出力

% 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

for t=1:T

for j=1:Num_PVG

P_PVG_yc(j,t)=int_pvg(j)*by_P_PVG_timeline(t);

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

for j=1:Node

sjdmjd=1;

for i=1:Node

if JD_gl(i,j)~=0

MJD(j,sjdmjd)=i;

sjdmjd=sjdmjd+1;

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

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

C=[C,C_3];

%% 二阶锥约束

C_4=[];

for t=1:T

for j=1:Node

for i=1:Node

if (JD_gl(i,j)~=0)

C_4=[C_4,norm([2*P_ij(JD_gl(i,j),t),2*Q_ij(JD_gl(i,j),t),I_ij(JD_gl(i,j),t)-U_j(j,t)],2)<=I_ij(JD_gl(i,j),t)+U_j(j,t)]

end

C=[C,C_4];

%% 节点电压、线路电流约束

C_5=[];

C_5=[C_5,Umin<=U_j<=Umax]

for i=1:Node

for j=1:Node

if JD_gl(i,j)~=0

if ismember(JD_gl(i,j),E0)==1

C_5=[C_5,0<=I_ij(JD_gl(i,j),:)<=(SLmax/12.66)^2]

else

C_5=[C_5,0<=I_ij(JD_gl(i,j),:)<=(SLmax/12.66)^2*x_E1_line(find(E1==JD_gl(i,j)))]

end

C=[C,C_5];

%% 变压器节点功率约束

C_6=[];

C_6=[C_6,P_TRmin<=P_TR<=P_TRmax];

C_6=[C_6,P_TRmin<=Q_TR<=P_TRmax];

C=[C,C_6];

%% 风光功率约束考虑是否建设

C_7=[];

% for t=1:T

% for j=1:Node

% if ismember(j,PVG_jr)==1

% C_7=[C_7,(1-a_max_PVG)*x_pvg(find(PVG_jr==j))*P_PVG_yc(find(PVG_jr==j),t)<=P_PVG(find(PVG_jr==j),t)<=x_pvg(find(PVG_jr==j))*P_PVG_yc(find(PVG_jr==j),t)];

% end

for t=1:T

for j=1:Num_PVG

C_7=[C_7,(1-a_max_PVG)*P_PVG_yc(j,t)<=P_PVG(find(PVG_jr==j),t)<=P_PVG_yc(j,t)];

end

C=[C,C_7];

%% 线路功率约束

C_8=[];

for i=1:Node

for j=1:Node

if JD_gl(i,j)~=0

if ismember(JD_gl(i,j),E0)==1

C_8=[C_8,-SLmax<=P_ij(JD_gl(i,j),:)<=SLmax];

C_8=[C_8,-SLmax<=Q_ij(JD_gl(i,j),:)<=SLmax];

else

C_8=[C_8,x_E1_line(find(E1==JD_gl(i,j)))*(-SLmax)<=P_ij(JD_gl(i,j),:)<=x_E1_line(find(E1==JD_gl(i,j)))*SLmax];

C_8=[C_8,x_E1_line(find(E1==JD_gl(i,j)))*(-SLmax)<=Q_ij(JD_gl(i,j),:)<=x_E1_line(find(E1==JD_gl(i,j)))*SLmax];

end

C=[C,C_8];

%% 目标函数及求解

obj_OPE1=0;

obj_OPE2=0;

obj_OPE3=0;

obj_OPE4=0;

for i=1:Num_Line

obj_OPE1=obj_OPE1+x_E1_line(i)*Cline(E1(i));

end

for t=1:T

% 线路损耗

obj_OPE2=obj_OPE2+Line_closs*(Rline*I_ij(:,t)/1000);

for i=1:Node

% 从主网成本

if ismember(i,TR_jr)==1

obj_OPE3= obj_OPE3+zw_buy_TR(t)*P_TR(find(TR_jr==i),t);

end

% if ismember(i,PVG_jr)==1

% obj_OPE4= obj_OPE4+ca_PVG*(P_PVG_yc(find(PVG_jr==i),t)-P_PVG(find(PVG_jr==i),t));

% end

end

for t=1:T

for j=1:Num_PVG

obj_OPE4= obj_OPE4+ca_PVG*(P_PVG_yc(j,t)-P_PVG(j,t));

end

% obj=obj_OPE1+obj_OPE2+obj_OPE3;

obj=obj_OPE2+obj_OPE3+obj_OPE4;

ops =sdpsettings('solver','Gurobi','verbose', 2, 'debug', 1);

% ops =sdpsettings('solver','Cplex');

result = optimize(C,obj,ops);

%% 以 m_ 开头保存计算结果

m_x_E1_line=value(x_E1_line);

m_x_pvg=value(x_pvg);

m_P_TR=value(P_TR);

m_Q_TR=value(Q_TR);

m_P_LOAD=value(P_LOAD);

m_P_PVG_yc=value(P_PVG_yc);

m_P_PVG=value(P_PVG);

%% 运行时间

toc; %tic用来保存当前时间，而后使用toc来记录程序完成时间