-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathLinearDistFlow_DNP_withV2G_draft.m
149 lines (148 loc) · 6.7 KB
/
LinearDistFlow_DNP_withV2G_draft.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
%% Distribution Network Planning (DNP) based on Linear Distflow Optimal Power Flow (OPF) Model in
% M. E. Baran and F. F. Wu, “Optimal capacitor placement on radial distribution systems,” IEEE Trans. Power Delivery, vol. 4, no. 1, pp. 725–734,1989.
% M. E. Baran and F. F. Wu, “Optimal sizing of capacitors placed on a radial distribution system,” IEEE Trans. Power Delivery, vol. 4, no. 1, pp. 735–743, 1989
% A compact form is described in Section III. of
% Exact convex relaxation of optimal power flow in radial networksL Gan, N Li, U Topcu, SH LowIEEE Transactions on Automatic Control 60 (1), 72-87, 2014
clear all
close all
clc
%% ***********Parameters **********
N=16; % number of load nodes
L=33; % number of dis. lines
Sbase=1e6; % unit:VA
Ubase=10e3; % unit:V
Ibase=Sbase/Ubase/1.732; % unit: A
Zbase=Ubase/Ibase/1.732; % unit: Ω
LineInf=xlsread('case data\16bus33lines.xlsx','F3:L35');
NodeInf=xlsread('case data\16bus33lines.xlsx','A3:D18');
s=LineInf(:,2);
t=LineInf(:,3);
N_subs=[13,14,15,16]; % Subs nodes
N_loads=1:12; % Load nodes
v_min=0.95^2;
v_max=1.05^2;
S_max=12 ; % max power in any distribution line
M=1e8;
Pr_e=0.9; % electric price for V2G ??? Pr_e=sdpvar(XXX) ???
% line investment cost denoted by
Cost=LineInf(:,7).*LineInf(:,4);
% make load (MVA) become P_load Q_load assuming a power factor n.
S=NodeInf(N_loads,4);
n=0.8; % power factor
P_load=S*n;
Q_load=S*sqrt(1-n^2);
% make impedance z become r and x assuming a rx rate n_rx
z=LineInf(:,4).*LineInf(:,6)/Zbase; % line impedance unit:p.u.
n_rx=1;
r=z/sqrt(1+n_rx^2);
x=r/n_rx;
%% ***********Graph of the Distribution System***********
G=graph(s,t); % G=(V,E), G:graph V:vertex E:edge
idxOut = findedge(G,s,t); % !!!!!!*************** index of edge is not the same with that of in mpc.branch !!!!!!*******
In=myincidence(s,t); % 节支关联矩阵，node-branch incidence matrix
% Plot the Graph
figure;
p=plot(G,'Layout','force');
p.XData=NodeInf(:,2);
p.YData=NodeInf(:,3);
hold on;
labelnode(p,N_subs,{'13','14','15','16'});
highlight(p,N_loads,'NodeColor','y','Markersize',20,'NodeFontSize',20);
highlight(p,N_subs,'Marker','s','NodeColor','c','Markersize',30,'NodeFontSize',40);
highlight(p,s,t,'EdgeColor','k','LineStyle','-.','LineWidth',2,'EdgeFontSize',8);
text(p.XData, p.YData, p.NodeLabel,'HorizontalAlignment', 'center','FontSize', 15); % put nodes' label in right position.
p.NodeLabel={};
hold off;
%% ***********Variable statement**********
v_i=sdpvar(N,1, 'full'); % v=|U|^2=U·U*, "*" means conjugate
l_ij=sdpvar(L,1, 'full'); % l=|I|^2=I·I*
P_ij=sdpvar(L,1, 'full');
Q_ij=sdpvar(L,1, 'full');
P_v2g=sdpvar(length(N_loads),1, 'full'); % V2G active power
Q_v2g=sdpvar(length(N_loads),1, 'full'); % V2G reactive power
P_shed=sdpvar(length(N_loads),1, 'full'); % shedded load
Q_shed=sdpvar(length(N_loads),1, 'full'); % shedded load
g_subs_P=sdpvar(length(N_subs),1,'full'); % generation from substations
g_subs_Q=sdpvar(length(N_subs),1,'full'); % generation from substations
y_ij=binvar(L,1,'full'); % decision variable denoting whether this line should be invested
y_v2g=binvar(length(N_loads),1,'full'); % decision variable denoting whether this V2G station should be invested
%% ***********Constraints*************
Cons=[];
%% 1. Power balance
% S_ij=s_i+\sum_(h:h→i) S_hi for any (i,j)∈E,
% denoted by node-branch incidence matrix I
Cons_S=[];
% for i=1:N
% Cons_S=[Cons_S,In*P_ij==[-(P_load-P_shed-P_dg);g_subs_P],In*Q_ij==[-(Q_load-Q_shed-Q_dg);g_subs_Q]];
% end
Cons_S=[In*P_ij==[-(P_load-P_shed-P_v2g);g_subs_P],In*Q_ij==[-(Q_load-Q_shed-Q_v2g);g_subs_Q], P_shed>=0,Q_shed==P_shed*sqrt(1-n^2)/n];
Cons=[Cons,Cons_S];
%% 2. Voltage Calculation
% v_i-v_j=2Re(z_ij·S_ij*)=2(r·P_ij+x·Q_ij)
% → |v_i-v_j-2(r·P_ij+x·Q_ij)|≤(1-y_ij)*M
Cons_V=[v_i(N_subs)==v_max,v_i];
Cons_V=[Cons_V,abs(In'*v_i-2*r.*P_ij-2*x.*Q_ij)<=(1-y_ij)*M];
Cons=[Cons,Cons_V];
%% 3. Voltage limits
% v_min<=v<=v_max
Cons=[Cons,v_min<=v_i<=v_max];
%% 4. P_ij^2+Q_ij^2<=S_ij_max
% Cons_PQ=[abs(P_ij)<=y_ij*S_max,abs(Q_ij)<=y_ij*S_max];
Cons_PQ=[P_ij.^2+Q_ij.^2<=y_ij*S_max.^2];
Cons=[Cons, Cons_PQ];
%% 5. P_ch.*P_dis=0 (P_ch=sdpvar(T,1,'full'),P_dis=sdpvar(T,1,'full')) ??? P_es∈[-P_ch,P_dis] ???
% 0=<P_ch<=z.*P_ch_max,0=<P_dis<=(1-z).*P_dis_max ,z=binvar(T,1,'full') (Big-M)
%% 6. Energy storage energy balance
% for ti=1~T
% E_ti=\sum(t=1~ti) (n_ch*P_ch-P_dis/n_dis)+E_0
% end
%% **********Objectives*************
Obj_inv=sum(Cost.*y_ij); %investment cost, related to line investment。。。。。。。+y_v2g
Obj_ope=M*sum(P_shed); %+Pr_e*sum(P_dg); Voltage deviation; Line loss
Obj=Obj_inv+Obj_ope;
%% ********* Solve the probelm
ops=sdpsettings('solver','gurobi', 'gurobi.Heuristics',0,'gurobi.Cuts',0); %,'usex0',1,'gurobi.MIPGap',5e-2,
sol=optimize(Cons,Obj,ops)
%% Save the solution with "s_" as start
s_y_ij=value(y_ij);
s_v_i=value(v_i);
s_P_ij=value(P_ij);
s_Q_ij=value(Q_ij);
s_P_shed=value(P_shed);
s_Q_shed=value(Q_shed);
s_g_subs_P=value(g_subs_P);
s_g_subs_Q=value(g_subs_Q);
s_Obj=value(Obj);
s_Obj=value(Obj_inv);
s_Obj=value(Obj_ope);
%% Print the results in the command line
display('————————规划方案如下——————————————');
display([' 建设线路方案： ',num2str(round(s_y_ij',2))]);
display([' (1/0 表示 建设/不建设 该线路, 线路编号参见输入文件)']);
display('————————规划建设成本如下——————————————');
display([' 建设配电线路成本： ',num2str(value(Obj_inv)),' 美元']);
display([' 失负荷成本: ',num2str(round(value(Obj_ope),2)),' 美元']);
display(['>> 规划建设总成本: ',num2str(value(Obj)),' 美元']);
%% Plot the results
Gi=graph(s,t);
os=0.1; % offset value to print text in the figure
figure;
pi=plot(Gi,'Layout','force');
pi.LineWidth(idxOut) = 10*abs(s_P_ij)/max(abs(s_P_ij))+0.01;
pi.XData=NodeInf(:,2);
pi.YData=NodeInf(:,3);
labelnode(pi,N_subs,{'13','14','15','16'});
highlight(pi,N_loads,'NodeColor','y','Markersize',20,'NodeFontSize',20);
highlight(pi,N_subs,'Marker','s','NodeColor','c','Markersize',30,'NodeFontSize',40);
text(pi.XData, pi.YData, pi.NodeLabel,'HorizontalAlignment', 'center','FontSize', 15); % put nodes' label in right position.
text(pi.XData+os, pi.YData+os, num2str([round(s_P_shed,2);zeros(length(N_subs),1)]),'HorizontalAlignment', 'center','FontSize', 20, 'Color','r'); % printf the shedded load (if any).
text(pi.XData-os, pi.YData-os, num2str(sqrt(s_v_i)),'HorizontalAlignment', 'center','FontSize', 15, 'Color','g'); % printf the nodes' voltage (sqrt(v_j)).
for l=1:L
Coor1_PQ=(pi.XData(s(l))+pi.XData(t(l)))/2;
Coor2_PQ=(pi.YData(s(l))+pi.YData(t(l)))/2;
text(Coor1_PQ, Coor2_PQ, num2str(round(s_P_ij(l)+s_Q_ij(l)*i,2)),'HorizontalAlignment', 'center','FontSize', 15, 'Color','b'); % printf the complex power in distribution lines(if any).
end
pi.NodeLabel={};
%% Save the results
warning off
save('LinearDistFlow_planning_results');