electrical engineering end of semester project

% Example Power Flow in a 3-bus Test System using Gauss-Seidel method Iter = 0; Maxiter = 20; converge = 0; tolerance = 0.001; j = sqrt(-1); % ---------- Input line impedances ------------- % Z = [0 0.0047 + 0.0474i 0.0062 + 0.0632i 0.0047 + 0.0474i 0 0.0047 + 0.0474i 0.0062 + 0.0632i 0.0047 + 0.0474i 0]; % ---------- Base Values ----------- % kVLL=345; MVA3Ph=100; Zbase=kVLL^2/MVA3Ph; XL_km=0.376; % ohm/km at 60 Hz RL_km= 0.037; B_km=4.5; % B in micro-mho/km Z13_ohm=(RL_km+j*XL_km)*200; %200 km long B13_Micro_Mho=4.5*200; %200 km long Z12_ohm=(RL_km+j*XL_km)*150; %150 km long B12_Micro_Mho=4.5*150; %150 km long Z23_ohm=(RL_km+j*XL_km)*150; %150 km long B23_Micro_Mho=4.5*150; %150 km long %--------- line impedances in per unit ----------% Z13=Z13_ohm/Zbase; Z12=Z12_ohm/Zbase; Z23=Z23_ohm/Zbase; %--------- susceptances in per unit ----------% B13=B13_Micro_Mho*Zbase*10^-6; B12=B12_Micro_Mho*Zbase*10^-6; B23=B23_Micro_Mho*Zbase*10^-6; %-------Finding the admittance matrix---------% Y(1,1)=1/Z12 + 1/Z13; Y(1,2)=-1/Z12; Y(1,3)=-1/Z13; Y(2,1)=-1/Z12; Y(2,2)=1/Z12 + 1/Z23; Y(2,3)=-1/Z23; Y(3,1)=-1/Z13; Y(3,2)=-1/Z23; Y(3,3)=1/Z13 + 1/Z23; G = real(Y); B = imag(Y); % ------- Initialize line variables --------- % I = zeros(1,length(Z)); P = zeros(1,length(Z)); Q = zeros(1,length(Z)); S = zeros(1,length(Z)); % -------- Input bus parameters --------- % v = ones(1,length(Z)); %angle = zeros(length(Z),1); Pload = zeros(1,length(Z)); Qload = zeros(1,length(Z)); Psched = zeros(1,length(Z)); Qsched = zeros(1,length(Z)); Bustype(1) = 'S'; % PQ bus Bustype(2) = 'G'; % PV bus Bustype(3) = 'L'; % PQ bus v(2) = 1.05; Pload(2) = -2; Pload(3) = 5; Qload(3) = 1; % ------ Set up power schedule --------- % totload = 0.0; for i = 1:length(Z) Psched(i)= -Pload(i); Qsched(i)= -Qload(i); totload = totload + Pload(i); end % --- Calculate P and Q at each bus ----- % for i = 1:length(Z) I(i) = 0.0 + sqrt(-1)*0.0; for j = 1:length(Z) I(i) = I(i) + Y(i,j)*v(j); end S(i) = v(i) * conj(I(i)); P(i) = real(S(i)); Q(i) = imag(S(i)); end % ------- Display bus voltages ---------- % fprintf('\n %s %2d %s \n', 'Iter ',Iter, ' Vreal Vimag Vmag Vangle P Q ') for i = 1:length(Z) vitermag = abs(v(i)); viterang = (180./pi)*atan2(imag(v(i)),real(v(i))); fprintf(' %s %2d %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f \n',' Bus ',i, real(v(i)), imag(v(i)), vitermag, viterang, P(i), Q(i) ) end while converge == 0 Iter = Iter + 1; MAXDP = 0.0; MAXDPbus = 0; MAXDQ = 0.0; MAXDQbus = 0; for i = 1 : length(Z) %---------calculate net P and Q at bus i--------% I(i) = 0.0 + sqrt(-1)*0.0; for j = 1:length(Z) I(i) = I(i) + Y(i,j)*v(j); end S(i) = v(i) * conj(I(i)); P(i) = real(S(i)); Q(i) = imag(S(i)); if Bustype(i) == 'G' Qsched(i) = Q(i); end deltap(i) = abs(P(i) - Psched(i)); deltaq(i) = 0.0; if Bustype(i) == 'L' deltaq(i) = abs(Q(i) - Qsched(i)); end if Bustype(i) == 'S' deltap(i) = 0.0; deltaq(i) = 0.0; end if Bustype(i) ~= 'S' if deltap(i) > MAXDP MAXDP = deltap(i); MAXDPbus = i; end if deltaq(i) > MAXDQ MAXDQ = deltaq(i); MAXDQbus = i; end %--------- Y * V for row i of Y matrix without Yii term---------% sum = 0.0; for j = 1:length(Z) if j ~= i sum = sum + Y(i,j)*v(j); end end vnew = (1.0/Y(i,i))*( ((Psched(i) - sqrt(-1)*Qsched(i))/(conj(v(i)))) - sum); v(i) = vnew; end end %------------ Print and save result from last iteration-----------% %-----calculate net P and Q at bus i-----% for i = 1:length(Z) I(i) = 0.0 + sqrt(-1)*0.0; for j = 1:length(Z) I(i) = I(i) + Y(i,j)*v(j); end S(i) = v(i) * conj(I(i)); P(i) = real(S(i)); Q(i) = imag(S(i)); if Bustype(i) == 'G' Qsched(i) = Q(i); end end fprintf('\n %s %2d %s \n', 'Iter ',Iter, ' Vreal Vimag Vmag Vangle P Q ') for i = 1:length(Z) vitermag = abs(v(i)); viterang = (180./pi)*atan2(imag(v(i)),real(v(i))); fprintf(' %s %2d %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f \n',' Bus ',i, real(v(i)), imag(v(i)), vitermag, viterang, P(i), Q(i) ) end % -------- check for convergence ---------- % if MAXDP < tolerance if MAXDQ < tolerance converge = 1; end end if Iter > Maxiter converge = 1; end end