-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathRecovery.m
executable file
·115 lines (82 loc) · 2.95 KB
/
Recovery.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
% "Recovery.m" recovers voltage phasors from the matrix solution and
% derives the generating power.
% SDP Solver of Optimal Power Flow: beta version
% Ramtin Madani([email protected])
% Morteza Ashraphijuo([email protected])
% Javad Lavaei([email protected])
% Columbia University
% Last Modified: October 07, 2014
function [V_rec, Sg_rec, Sb_rec, sf_rec, st_rec, cost_rec] = Recovery(V2,W,Pd,Qd, ...
ng,incidentG,Qmax,Qmin,Pmax,Pmin,correction,fbusN,tbusN,activeL,busType,Yf,Yt,Ybus,cc, ...
solver_rec,c2,c1,c0)
nc = length(cc);
nb = size(V2,1);
nl = length(fbusN);
Vm_rec = sqrt(V2);
phaseL = zeros(nl,nc);
for rr = 1 : nc
in1 = nb*(rr-1);
in2 = setdiff(activeL, cc{rr});
phaseL(in2,rr) = angle( diag( W(in1 + fbusN(in2), in1 + tbusN(in2)) ) );
end
source = find(busType == 3);
cvx_begin
if strcmp(solver_rec,'sdpt3')
cvx_solver sdpt3
elseif strcmp(solver_rec,'sedumi')
cvx_solver sedumi
elseif strcmp(solver_rec,'mosek')
cvx_solver mosek
end
cvx_precision best
variable Va_rec(nb,nc)
expression phaseL1(nl,nc)
for rr = 1 : nc
in2 = setdiff(activeL, cc{rr});
phaseL1(in2,rr) = Va_rec(fbusN(in2),rr) - Va_rec(tbusN(in2),rr);
end
cost_lp = sum(sum(abs(phaseL1 - phaseL)));
minimize( cost_lp );
subject to
Va_rec(source,:) == 0;
cvx_end
V_rec = Vm_rec .* exp(Va_rec*1i);
Sb_rec = zeros(nb,nc);
sf_rec = zeros(nl,nc);
st_rec = zeros(nl,nc);
for rr = 1 : nc
ex = setdiff(1 : nl, cc{rr});
Sb_rec(:,rr) = Pd + 1i*Qd + V_rec(:,rr) .* conj(Ybus{rr}*V_rec(:,rr));
sf_rec(ex,rr) = V_rec(fbusN(ex),rr) .* conj(Yf(ex,:) * V_rec(:,rr));
st_rec(ex,rr) = V_rec(tbusN(ex),rr) .* conj(Yt(ex,:) * V_rec(:,rr));
end
gb = find((sum(incidentG) ~= 0));
cvx_begin
if strcmp(solver_rec,'sdpt3')
cvx_solver sdpt3
elseif strcmp(solver_rec,'sedumi')
cvx_solver sedumi
elseif strcmp(solver_rec,'mosek')
cvx_solver mosek
end
cvx_precision best
variable Sg_rec(ng,nc) complex
cost_rec = sum(c2.*real(Sg_rec(:,1)).^2 + c1.*real(Sg_rec(:,1)) + c0);
minimize( cost_rec + (5*10^8) * ...
max([0;max(max(real(Sg_rec) - Pmax * ones(1,nc))); ...
max(max(Pmin * ones(1,nc) - real(Sg_rec))); ...
max(max(imag(Sg_rec) - Qmax * ones(1,nc))); ...
max(max(Qmin * ones(1,nc) - imag(Sg_rec))); ...
max(max(abs(Sb_rec(gb,:) - incidentG(:,gb).' * Sg_rec)))]));
% minimize( cost_rec );
subject to
% abs(Sb_rec(gb) - incidentG(:,gb).' * Sg_rec) <= 0;
% real(Sg_rec) <= Pmax * ones(1,nc);
% real(Sg_rec) >= Pmin * ones(1,nc);
% imag(Sg_rec) <= Qmax * ones(1,nc);
% imag(Sg_rec) >= Qmin * ones(1,nc);
for rr = 2 : nc
abs(real(Sg_rec(:,rr) - Sg_rec(:,1))) <= correction;
end
cvx_end
end