网架数据展示:
完整程序:
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=[];
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
end
end
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
end
end
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
% end
% 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
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
end
end
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
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
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来记录程序完成时间