Chetan Shidling Staff answered 4 years ago
Here is the code:
clc;
ne=4; nn=4; nref=0;
net=[ 0 1;
0 2;
1 2;
2 3];
flag=input(‘Enter \n1-System with mutual coupling \n2-System without
mutual coupling\n’);
if flag==1
z= [(0.0+0.12i) (0.0+0.0i) (0.0+0.0i) (0.0+0.0i);
(0.0+0.0i) (0.0+0.2i) (0.0+0.1i) (0.0+0.0i);
(0.0+0.0i) (0.0+0.1i) (0.0+0.1i) (0.0+0.0i);
(0.0+0.0i) (0.0+0.0i) (0.0+0.0i) (0.0+0.24i)];
else
z= [(0.0+0.12i) (0.0+0.0i) (0.0+0.0i) (0.0+0.0i);
(0.0+0.0i) (0.0+0.2i) (0.0+0.0i) (0.0+0.0i);
(0.0+0.0i) (0.0+0.0i) (0.0+0.1i) (0.0+0.0i);
(0.0+0.0i) (0.0+0.0i) (0.0+0.0i) (0.0+0.24i)];
end
ebus=[ (0.98-0.05i);
(0.95-0.1i);
(0.92-0.12i)];
a=zeros(ne,nn-1);
for i=1:ne
nf=net(i,1);
nt=net(i,2);
if nf~=nref
a(i,nf)=1;
end
if nt~=nref
a(i,nt)=-1;
end
end
y=inv(z);
ybus=transpose(a)*y*a;
ibus=ybus*ebus;
sbus=conj(ebus).*ibus;
disp(‘ BUS INCIDENCE MATRIX (A): ‘);
disp(a);
disp(‘ PRIMITIVE IMPEDENCE MATRIX (y): ‘);
disp(y);
disp(‘ BUS ADMITTANCE MATRIX [Ybus]: ‘);
disp(ybus);
disp(‘ BUS CURRENTS [Ibus]: ‘);
disp(ibus);
disp(‘ COMPLEX POWER AT BUSES [Sbus]: ‘);
disp(sbus);