11-25-2015, 05:16 PM
Code:
linedata=[1 2 0.01008, 0.05040, 3.815629, -19.078144, 10.25, 0.05125;
1 3 0.00744, 0.03720, 5.169561, -25.847809, 7.75, 0.03875;
2 4 0.00744, 0.03720, 5.169561, -25.847809, 7.75, 0.03875;
3 4 0.01272, 0.06360, 3.023705, -15.118528, 12.75, 0.06375];
busdata=[1 0, 0, 50, 30.99, 1.00, 0 1;
2 0, 0, 170, 105.35, 1.00, 0 2;
3 0, 0, 200, 123.94, 1.00, 0 2;
4 318, 0, 80, 49.58, 1.02, 0 3]
% Bus Type: 1.Slack Bus 2.PQ Bus 3.PV Bus
ss=1i*linedata(:,8);
y=linedata(:,5)+1i*linedata(:,6);
totalbuses = max(max(linedata(:,1)),max(linedata(:,2))); % total buses
totalbranches = length(linedata(:,1)); % no. of branches
ybus = zeros(totalbuses,totalbuses);
for b=1:totalbranches
ybus((linedata(b,1)),(linedata(b,2)))=-y(b);
ybus((linedata(b,2)),(linedata(b,1))) =ybus((linedata(b,1)),(linedata(b,2)));
end
for c=1:totalbuses
for d=1:totalbranches
if linedata(d,1) == c || linedata(d,2) == c
ybus(c,c) = ybus(c,c) + y(d) + ss(d);
end
end
end
disp('TABLE 9.3 PAGE # 338 BUS ADMITTANCE MATRIX FOR EXAMPLE 9.2')
ybus
z=zeros(totalbuses,4);
busnumber=busdata(:,1);
PG=busdata(:,2);
QG=busdata(:,3);
PL=busdata(:,4);
QL=busdata(:,5);
V=busdata(:,6);
VV=V;
ANG=busdata(:,7);
type = busdata(:,8);
P = (PG-PL)./100; % per unit active power at buses
Q = (QG-QL)./100; % per unit reactive power at buses
tol=1;
iter=0;
kk=0.0001;
while tol > kk
for i = 2:totalbuses
YV = 0;
for k = 1:totalbuses
if i~=k
YV = YV + ybus(i,k)* V(k); % multiplying admittance & voltage
YV;
end
YV;
end
if busdata(i,8) == 3 %Calculating Qi for PV bus
%Q(i) = -imag(conj(V(i))*(YV + ybus(i,i)*V(i)));
Q(i) = -imag(conj(V(i))*(YV + ybus(i,i)*V(i)));
busdata(i,3)=Q(i);
end
% end
V(i) = (1/ybus(i,i))*((P(i)-j*Q(i))/conj(V(i)) - YV); % Compute Bus Voltages.
% Calculating Corrected Voltage for PV bus
if busdata(i,8) == 3
vc(i)=abs(VV(i))*(V(i)/abs(V(i)));
busdata(i,6)=vc(i);
V(i)=vc(i);
end
% Calculating Accelerated Voltage for PQ bus
if busdata(i,8) == 2
VACC(i)= VV(i)+(V(i)-VV(i));
busdata(i,6)=VACC(i);
V(i)=VACC(i);
end
%V(i)=V;
end
iter = iter + 1; % Increment iteration count.
tol = max(abs(abs(V) - abs(VV))); % Calculate tolerance.
VV = V;
end
Q;
iter;
YV;
V;
%real(VACC')
z(1:totalbuses,1)=busdata(:,1);
z(1:totalbuses,2)=busdata(:,8);
z(1:totalbuses,3)=abs(busdata(:,6));
z(1:totalbuses,4)=radtodeg(angle(V));
disp(' Bus No. Bus Type Voltage Angle ');
z
hp prime porting code ...tried but got different results
Code:
#cas
gausss(M1,M3):=
BEGIN
//M1=linedata, M2=ybus M3= busdata
ss:= .* col(M1,5);
y:= col(M1,3);
y:=1 ./ y;
totalbuses:= MAX(MAX(col(M1,1)),MAX(col(M1,2)));
totalbranches:= length(M1);
M2:=makemat(totalbuses,totalbuses);
FOR b FROM 1 TO totalbranches DO
M2(M1(b,1),M1(b,2)):=-y(b);
M2((M1(b,2)),(M1(b,1))):=M2((M1(b,1)),(M1(b,2)));
END;
FOR c FROM 1 TO totalbuses DO
FOR d FROM 1 TO totalbranches DO
IF M1(d,1)==c OR M1(d,2) == c THEN
M2(c,c) := M2(c,c) + y(d)+ss(d);
END;
END;
END;
M2;
M9:=makemat(totalbuses,4);
busnumber:=col(M3,1);
pg:=col(M3,2);
qg:=col(M3,3);
pl:=col(M3,4);
ql:=col(M3,5);
v:= col(M3,6);
vv:=v;
ANG:= col(M3,7);
types:= col(M3,8);
p := (pg - pl) ./100;
q := (qg - ql) ./100;
toll:=1;
iter:=0;
kk:=0.0001
alfa:=1.6;
WHILE toll > kk DO
FOR j FROM 2 TO totalbuses DO
YV := 0;
FOR k FROM 1 TO totalbuses DO
IF j≠k THEN
YV := YV + M2(j,k)* v(k);
END;
END;
//Calculating Qi for PV bus
IF M3(j,8) == 2 THEN
q(j) := -IM((CONJ(v(j)))*(YV + M2(j,j)*v(j)));
M3(j,3):=q(j);
END;
v(j) := (1/M2(j,j))*((p(j)-*q(j))/CONJ(v(j)) - YV); //Compute Bus Voltages.
//Calculating Corrected Voltage for PV bus
IF M3(j,8) == 3 THEN
vc(j):=(abs(vv(j)))*(v(j)/abs(v(j)));
M3(j,6):=vc(j);
v(j):=vc(j);
END;
//Calculating Accelerated Voltage for PQ bus
IF M3(j,8) == 2 THEN
VACC(j):= vv(j)+alfa*(v(j)-vv(j));
M3(j,6):=VACC(j);
v(j):=VACC(j);
END;
END;
iter := iter + 1;
toll :=abs(abs(MAX(v))-abs(MAX(vv)));
vv := v;
END;
q;
YV;
v;
END;
#end