-
Notifications
You must be signed in to change notification settings - Fork 0
/
nrlfppg.m
167 lines (149 loc) · 5.27 KB
/
nrlfppg.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
% Program for Newton-Raphson Load Flow Analysis..
%nbus = 14;
Y = ybusppg(nbus); % Calling ybusppg.m to get Y-Bus Matrix..
assignin('base','Y',Y);
Z = inv(Y); % Bus Impedance Matrix
assignin('base','Z',Z);
% busd = busdatas(nbus); % Calling busdatas..
BMva = 100; % Base MVA..
bus = busd(:,1); % Bus Number..
type = busd(:,2); % Type of Bus 1-Slack, 2-PV, 3-PQ..
V = busd(:,3); % Specified Voltage..
del = busd(:,4); % Voltage Angle..
Pg = busd(:,5)/BMva; % PGi..
Qg = busd(:,6)/BMva; % QGi..
Pl = busd(:,7)/BMva; % PLi..
Ql = busd(:,8)/BMva; % QLi..
Qmin = busd(:,9)/BMva; % Minimum Reactive Power Limit..
Qmax = busd(:,10)/BMva; % Maximum Reactive Power Limit..
P = Pg - Pl; % Pi = PGi - PLi..
Q = Qg - Ql; % Qi = QGi - QLi..
Psp = P; % P Specified..
Qsp = Q; % Q Specified..
G = real(Y); % Conductance matrix..
B = imag(Y); % Susceptance matrix..
pv = find(type == 2 | type == 1); % PV Buses..
pq = find(type == 3); % PQ Buses..
npv = length(pv); % No. of PV buses..
npq = length(pq); % No. of PQ buses..
Tol = 1;
Iter = 1;
while (Tol > 1e-7) % Iteration starting..
P = zeros(nbus,1);
Q = zeros(nbus,1);
% Calculate P and Q
for i = 1:nbus
for k = 1:nbus
P(i) = P(i) + V(i)* V(k)*(G(i,k)*cos(del(i)-del(k)) + B(i,k)*sin(del(i)-del(k)));
Q(i) = Q(i) + V(i)* V(k)*(G(i,k)*sin(del(i)-del(k)) - B(i,k)*cos(del(i)-del(k)));
end
end
% Checking Q-limit violations..
if Iter <= 7 && Iter > 2 % Only checked up to 7th iterations..
for n = 2:nbus
if type(n) == 2
QG = Q(n)+Ql(n);
if QG < Qmin(n)
V(n) = V(n) + 0.01;
elseif QG > Qmax(n)
V(n) = V(n) - 0.01;
end
end
end
end
% Calculate change from specified value
dPa = Psp-P;
dQa = Qsp-Q;
k = 1;
dQ = zeros(npq,1);
for i = 1:nbus
if type(i) == 3
dQ(k,1) = dQa(i);
k = k+1;
end
end
dP = dPa(2:nbus);
M = [dP; dQ]; % Mismatch Vector
% Jacobian
% J1 - Derivative of Real Power Injections with Angles..
J1 = zeros(nbus-1,nbus-1);
for i = 1:(nbus-1)
m = i+1;
for k = 1:(nbus-1)
n = k+1;
if n == m
for n = 1:nbus
J1(i,k) = J1(i,k) + V(m)* V(n)*(-G(m,n)*sin(del(m)-del(n)) + B(m,n)*cos(del(m)-del(n)));
end
J1(i,k) = J1(i,k) - V(m)^2*B(m,m);
else
J1(i,k) = V(m)* V(n)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-del(n)));
end
end
end
% J2 - Derivative of Real Power Injections with V..
J2 = zeros(nbus-1,npq);
for i = 1:(nbus-1)
m = i+1;
for k = 1:npq
n = pq(k);
if n == m
for n = 1:nbus
J2(i,k) = J2(i,k) + V(n)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n)));
end
J2(i,k) = J2(i,k) + V(m)*G(m,m);
else
J2(i,k) = V(m)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n)));
end
end
end
% J3 - Derivative of Reactive Power Injections with Angles..
J3 = zeros(npq,nbus-1);
for i = 1:npq
m = pq(i);
for k = 1:(nbus-1)
n = k+1;
if n == m
for n = 1:nbus
J3(i,k) = J3(i,k) + V(m)* V(n)*(G(m,n)*cos(del(m)-del(n)) + B(m,n)*sin(del(m)-del(n)));
end
J3(i,k) = J3(i,k) - V(m)^2*G(m,m);
else
J3(i,k) = V(m)* V(n)*(-G(m,n)*cos(del(m)-del(n)) - B(m,n)*sin(del(m)-del(n)));
end
end
end
% J4 - Derivative of Reactive Power Injections with V..
J4 = zeros(npq,npq);
for i = 1:npq
m = pq(i);
for k = 1:npq
n = pq(k);
if n == m
for n = 1:nbus
J4(i,k) = J4(i,k) + V(n)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-del(n)));
end
J4(i,k) = J4(i,k) - V(m)*B(m,m);
else
J4(i,k) = V(m)*(G(m,n)*sin(del(m)-del(n)) - B(m,n)*cos(del(m)-del(n)));
end
end
end
J = [J1 J2; J3 J4]; % Jacobian Matrix..
X = inv(J)*M; % Correction Vector
dTh = X(1:nbus-1); % Change in Voltage Angle..
dV = X(nbus:end); % Change in Voltage Magnitude..
% Updating State Vectors..
del(2:nbus) = dTh + del(2:nbus); % Voltage Angle..
k = 1;
for i = 2:nbus
if type(i) == 3
V(i) = dV(k) + V(i); % Voltage Magnitude..
k = k+1;
end
end
Iter = Iter + 1;
Tol = max(abs(M)); % Tolerance..
end
loadflow(nbus,V,del,BMva); % Calling Loadflow.m..
assignin('base','Iter',Iter);