forked from eranschweitzer/distflow
-
Notifications
You must be signed in to change notification settings - Fork 0
/
distflow_lossy.m
177 lines (162 loc) · 5.91 KB
/
distflow_lossy.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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
function [v, Pf, Qf] = distflow_lossy(mpc, opt)
%%% single phase DistFlow solution
%%% [v, Pf, Qf] = distflow_lossy(mpc)
%%% [v, Pf, QF] = distflow_lossy(mpc, opt)
%%% result = distflow_lossy(mpc, opt)
%%%
%%% INPUTS:
%%% mpc: matpower case. Must be radial, with consecutively numbered
%%% branches, and the source node should be bus 1.
%%% opt: options structure
%%% .alpha: lossless parameter either scalar or
%%% range [alphamin alphamax].
%%% .alpha_method: parametrization method (1-7):
%%% 1: scalar alpha
%%% 2: linear interpolation based on branch impedance
%%% 3: quadradic interpolation based on branch impedance
%%% 4: linear interpolation based on downstream power
%%% 5: quadratic interpolation based on downstream power
%%% 6: linear interpolation based on branch impedance TIMES downstream
%%% power
%%% 7: quadratic interpolation based on branch impedance TIMES downstream
%%% power
%%% 8: alpha = 1/2 - opt.alpha*(power^2)*(r^2+x^2)
%%% OUTPUTS:
%%% v: vector of bus voltage magnitudes (per unit)
%%% Pf: vector of sending end branch real power flows (MW)
%%% Qf: vector of sending end branch reactive flows (MVAr)
%%% result: Copy of mpc with,
%%% result.bus(:, VM) = v
%%% result.branch(:, PF) = Pf
%%% result.branch(:, QF) = Qf
%% input check
define_constants;
nb = size(mpc.bus,1);
nl = size(mpc.branch,1);
if ~all(mpc.bus(:,BUS_I) == (1:nb)')
error('distflow_lossy: bus number must be consequtive starting at 1.')
end
if (nb - (nl + 1)) ~= 0
error('distflow_lossy: case must be radial but nb=%d and nl=%d', nb, nl)
end
if mpc.bus(1,BUS_TYPE) ~= REF
error('distflow_lossy: source bus must be first bus.')
end
%% options structure
if nargin < 2
opt = [];
end
if isscalar(opt) && ~isstruct(opt)
opt = struct('alpha', opt, 'alpha_method', 1);
end
opt = optdefaults(opt);
%% operator setup
NumberOfNodes=size(mpc.bus,1);
FromNode=mpc.branch(:,F_BUS);
ToNode=mpc.branch(:,T_BUS);
Branch = 1:length(FromNode);
NumberOfBranch=length(Branch);
F = sparse(Branch,FromNode,1,NumberOfBranch,NumberOfNodes);
T = sparse(Branch,ToNode,1,NumberOfBranch,NumberOfNodes);
M0 = F-T;
M = M0(:,2:end);
r = mpc.branch(:,BR_R);
x = mpc.branch(:,BR_X);
Rline = sparse(1:NumberOfBranch, 1:NumberOfBranch,r);
Xline = sparse(1:NumberOfBranch, 1:NumberOfBranch,x);
TFT = T*F';
I = speye(NumberOfBranch);
slack_voltage = mpc.bus(mpc.bus(:,BUS_TYPE) == REF, VM);
pc = mpc.bus(:,PD)/mpc.baseMVA;
qc = mpc.bus(:,QD)/mpc.baseMVA;
B = (I-TFT)^(-1);
C = 2*(Rline*B*Rline + Xline*B*Xline) - (Rline*Rline + Xline*Xline);
%% setup alpha
mtd = opt.alpha_method;
switch mtd
case 1 %isscalar(opt.alpha)
if ~isscalar(opt.alpha)%mtd ~= 1
% warning('distflow_lossy: scalar alpha (%0.4f) was passed but alpha_method=%d.\n\tScalar alpha only makes sense with method 1. Input method selection ignored.',...
% opt.alpha, mtd)
warning('distflow_lossy: alpha_method 1 selected but a non scalar alpha was passed.\n Using first entry alpha=%0.4f', opt.alpha(1))
opt.alpha = opt.alpha(1);
end
alpha = opt.alpha*I;
% tau =(Rline*(B - I)*Rline+Xline*(B - I)*Xline)*(Rline*Rline+Xline*Xline)^(-1);
% K = (opt.alpha*I - (1-2*opt.alpha)*tau);
% Kinv = K^(-1);
case num2cell(2:7)
% methods:
% 2: linear interpolation based on branch impedance
% 3: quadradic interpolation based on branch impedance
% 4: linear interpolation based on downstream power
% 5: quadratic interpolation based on downstream power
% 6: linear interpolation based on branch impedance TIMES downstream
% power
% 7: quadratic interpolation based on branch impedance TIMES downstream
% power
switch mtd
case {2,3}
v = sqrt(r.^2 + x.^2);
case {4,5}
v = abs(B*T*(pc + 1i*qc));
case {6,7}
v = sqrt(r.^2 + x.^2).*abs(B*T*(pc + 1i*qc));
end
if ismember(mtd, [3,5,7])
v = v.^2;
end
vmax = max(v); vmin = min(v);
slope = (opt.alpha(2) - opt.alpha(1))/(vmin - vmax);
intercept = opt.alpha(2) - slope*vmin;
alpha = sparse(1:nl,1:nl,slope*v + intercept, nl, nl);
case 8
% 8: alpha = 1/2 - opt.alpha*(power^2)*(r^2+x^2)
if ~isscalar(opt.alpha)%mtd ~= 1
warning('distflow_lossy: alpha_method 8 selected but a non scalar alpha was passed.\n Using first entry alpha=%0.4f', opt.alpha(1))
opt.alpha = opt.alpha(1);
end
v = (r.^2 + x.^2).*abs(B*T*(pc + 1i*qc)).^2;
atmp = 0.5 - opt.alpha*v;
% catch extreme values
atmp(atmp <= 0) = 0.5;
atmp(atmp >= 1) = 0.5;
alpha = sparse(1:nl,1:nl,atmp, nl, nl);
otherwise
error('distflow_lossy: method %d not implemented.', mtd)
end
beta = I - C*(Rline*Rline+Xline*Xline)^(-1)*(I - 2*alpha);
%% solution
Rlossy = M\(beta\(Rline*B*T)*2);
Xlossy = M\(beta\(Xline*B*T)*2);
v = [slack_voltage^2 ;slack_voltage^2+Rlossy*pc+Xlossy*qc];
c = (Rline*Rline+Xline*Xline)^(-1)*(I - 2*alpha)*M0*v;
v = sqrt(v);
Pf = mpc.baseMVA*B*(T*pc + Rline*c);
Qf = mpc.baseMVA*B*(T*qc + Xline*c);
if nargout == 1
result = mpc;
result.bus(:,VM) = v;
result.branch(:,PF) = Pf;
result.branch(:,QF) = Qf;
v = result;
end
function opt = optdefaults(opt)
optd = struct('alpha', 0.5, 'alpha_method', 1);
if isempty(opt)
opt = optd;
else
opt = struct_compare(optd, opt);
end
function b = struct_compare(a, b)
% compares structure b to structure a.
% if b lacks a field in a it is added
% this is performed recursively, so if if a.x is a structure
% and b has field x, the the function is called on (a.x, b.x)
for f = fieldnames(a).'
if ~isfield(b, f{:})
b.(f{:}) = a.(f{:});
elseif isstruct(a.(f{:}))
b.(f{:}) = struct_compare(a.(f{:}), b.(f{:}));
end
end