HT4T

Gauss-Siedel

13 Downloads

Updated 27 Oct 2022

View License

Mainfunction file:
busdata=[1 1 1.025 0 0 0 0 0
2 0 1.0 0 400 200 0 0
3 1 1.0 0 0 0 300 0 ];
linedata=[1 2 0 0.025 0 1
1 3 0 0.05 0 1
2 3 0 0.025 0 1];
lfybus
lfgauss
busout
lineflow
Lfybus file:
% Y Bus admittance matrix
j=sqrt(-1); i = sqrt(-1);
fb = linedata(:,1); % From bus number
tb = linedata(:,2); % To bus number
R = linedata(:,3); % Resistance, R
X = linedata(:,4); % Reactance, X...
b = j*linedata(:,5); % Shunt Admittance, B/2
Z = R + j*X; % Z matrix
y= 1./Z; %branch admittance
nbranch=length(linedata(:,1)); % no. of branches
nbus = max(max(fb), max(tb)); % no. of buses
% Forming the Y Bus Matrix
for n = 1:nbranch
Ybus=zeros(nbus,nbus); % initialize Ybus to zero
% Formation of the off diagonal elements
for k=1:nbranch;
Ybus(fb(k),tb(k))=Ybus(fb(k),tb(k))-y(k);
Ybus(tb(k),fb(k))=Ybus(fb(k),tb(k));
end
end
% Formation of the diagonal elements
for n=1:nbus
for k=1:nbranch
if fb(k)==n | tb(k)==n
Ybus(n,n) = Ybus(n,n)+y(k) + b(k);
else, end
end
end
Lfgauss file:
basemva = 100; %Base MVA
tolerance = 0.001; %Tolerance
mi = 100; %Maximum Iterations
af = 1.8; %Acceleration factor
% Keys for check purposes
Vm=0; delta=0; yload=0; deltad =0;
nbus = length(busdata(:,1));
for k=1:nbus
n=busdata(k,1);
bt(n)=busdata(k,2);
Vm(n)=busdata(k,3);
delta(n)=busdata(k, 4);
Pd(n)=busdata(k,5);
Qd(n)=busdata(k,6);
Pg(n)=busdata(k,7);
Qg(n) = busdata(k,8);
if Vm(n) <= 0
Vm(n) = 1.0; V(n) = 1 + j*0;
else delta(n) = pi/180*delta(n);
V(n) = Vm(n)*(cos(delta(n)) + j*sin(delta(n)));
P(n)=(Pg(n)-Pd(n))/basemva;
Q(n)=(Qg(n)-Qd(n))/basemva;
S(n) = P(n) + j*Q(n);
end
DV(n)=0;
end
num = 0; AcurBus = 0; converge = 1;
Vc = zeros(nbus,1)+j*zeros(nbus,1); Sc = zeros(nbus,1)+j*zeros(nbus,1);
iter=0;
maxerror=10;
sumc1=0;
sumc2=0;
sumc3=0;
sumc4=0;
while maxerror >= tolerance & iter <= mi
iter=iter+1;
for n = 1:nbus;
YV = 0+j*0;
for L = 1:nbranch;
if fb(L) == n, k=tb(L);
YV = YV + Ybus(n,k)*V(k);
elseif tb(L) == n, k=fb(L);
YV = YV + Ybus(n,k)*V(k);
end
end
Sc = conj(V(n))*(Ybus(n,n)*V(n) + YV) ;
Sc = conj(Sc);
DP(n) = P(n) - real(Sc);
DQ(n) = Q(n) - imag(Sc);
if bt(n) == 1
S(n) =Sc; P(n) = real(Sc); Q(n) = imag(Sc); DP(n) =0; DQ(n)=0;
Vc(n) = V(n);
elseif bt(n) == 2
Q(n) = imag(Sc); S(n) = P(n) + j*Q(n);
Qgc = Q(n)*basemva + Qd(n) ;
end
if bt(n) ~= 1
Vc(n) = (conj(S(n))/conj(V(n)) - YV )/ Ybus(n,n);
else, end
if bt(n) == 0
V(n) = V(n) + af*(Vc(n)-V(n));
elseif bt(n) == 2
VcI = imag(Vc(n));
VcR = sqrt(Vm(n)^2 - VcI^2);
Vc(n) = VcR + j*VcI;
V(n) = V(n) + af*(Vc(n) -V(n));
end
end
maxerror=max( max(abs(real(DP))), max(abs(imag(DQ))) );
if iter == mi & maxerror > tolerance
fprintf('\nWARNING: Iterative solution did not converged after ')
fprintf('%g', iter), fprintf(' iterations.\n\n')
fprintf('Press Enter to terminate the iterations and print the results \n')
converge = 0; pause, else, end
end
if converge ~= 1
tech= (' ITERATIVE SOLUTION DID NOT CONVERGE'); else,
tech=(' Power Flow Solution - Gauss-Seidel Method');
end
k=0;
for n = 1:nbus
Vm(n) = abs(V(n)); deltad(n) = angle(V(n))*180/pi;
if bt(n) == 1
S(n)=P(n)+j*Q(n);
Pg(n) = P(n)*basemva + Pd(n);
Qg(n) = Q(n)*basemva + Qd(n) ;
k=k+1;
Pgg(k)=Pg(n);
elseif bt(n) ==2
k=k+1;
Pgg(k)=Pg(n);
S(n)=P(n)+j*Q(n);
Qg(n) = Q(n)*basemva + Qd(n) ;
end
sumc1 = sumc1 + Pg(n);
sumc2 = sumc2+ Qg(n);
sumc3 = sumc3 + Pd(n);
sumc4 = sumc4 + Qd(n);
yload(n) = (Pd(n)- j*Qd(n))/(basemva*Vm(n)^2);
end
busdata(:,3)=Vm'; busdata(:,4)=deltad';
Busout file:
disp(tech)
fprintf(' %g Iterations \n\n', iter)
head =[' Bus Voltage Angle ------Load------ ---Generation--- '
' No. Mag. Degree MW Mvar MW Mvar '
' '];
disp(head)
for n=1:nbus
fprintf(' |%5g', n), fprintf(' |%7.3f', Vm(n)),
fprintf(' |%8.3f', deltad(n)), fprintf(' |%9.3f', Pd(n)),
fprintf(' |%9.3f', Qd(n)), fprintf(' |%9.3f', Pg(n)),
fprintf(' |%9.3f \n', Qg(n)),
end
fprintf(' \n'), fprintf(' Total ')
fprintf(' %9.3f', sumc3), fprintf(' %9.3f', sumc4),
fprintf(' %9.3f', sumc1), fprintf(' %9.3f\n\n',sumc2),
Lineflow file:
SLT = 0;
fprintf('\n')
fprintf(' Line Flow and Losses \n\n')
fprintf(' --Line-- Power at bus & line flow --Line loss-- \n')
fprintf(' from to MW Mvar MVA MW Mvar \n')
for n = 1:nbus
busprt = 0;
for L = 1:nbranch;
if busprt == 0
fprintf(' \n'), fprintf('%6g', n), fprintf(' %9.3f', P(n)*basemva)
fprintf('%9.3f', Q(n)*basemva), fprintf('%9.3f\n', abs(S(n)*basemva))
busprt = 1;
else, end
if fb(L)==n k = tb(L);
In = (V(n) - V(k))*y(L) + b(L);
Ik = (V(k) - V(n))*y(L) + b(L)*V(k);
Snk = V(n)*conj(In)*basemva;
Skn = V(k)*conj(Ik)*basemva;
SL = Snk + Skn;
SLT = SLT + SL;
elseif tb(L)==n k = fb(L);
In = (V(n) - V(k))*y(L) + b(L)*V(n);
Ik = (V(k) - V(n))*y(L) + b(L);
Snk = V(n)*conj(In)*basemva;
Skn = V(k)*conj(Ik)*basemva;
SL = Snk + Skn;
SLT = SLT + SL;
else, end
if fb(L)==n | tb(L)==n
fprintf('%12g', k),
fprintf('%9.3f', real(Snk)), fprintf('%9.3f', imag(Snk))
fprintf('%9.3f', abs(Snk)),
fprintf('%9.3f', real(SL)),
if fb(L) ==n
fprintf('%9.3f \n', imag(SL))
else, fprintf('%9.3f\n', imag(SL))
end
else, end
end
end
SLT = SLT/2;
fprintf(' \n'), fprintf(' Total loss ')
fprintf('%9.3f', real(SLT)), fprintf('%9.3f\n', imag(SLT))
clear Ik In SL SLT Skn Snk

Cite As

Guru Prasaath (2023). HT4T (https://www.mathworks.com/matlabcentral/fileexchange/119663-ht4t), MATLAB Central File Exchange. Retrieved .

MATLAB Release Compatibility
Created with R2022b
Compatible with any release
Platform Compatibility
Windows macOS Linux
Tags Add Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
Version Published Release Notes
1.0.0