clear
close all
clc

U1  = 100;
R1  = 20;
R2  = 30;
R3  = 10;
R4  = 20;
XL1 = 10i;
XL2 = 30i;
XC1 = -40i;
II  = 2*exp(40i);   % dont use 'I' as a variable name matlab will confuse it with 'i' (imaginary part)

syms I1 I2 I3 I4 I6 Uj

% equations from second Kirchhoff law
eq1='U1-I1*R1-I3*R2-I4*R3-I1*XC1=0';
eq2='-I2*XL1+II*R4-Uj+I3*R2=0';
eq3='Uj-II*R4-I6*XL2+I4*R3=0';

% from first Kirchhoff law
eq4='I1-I3-I2=0';
eq5='I3-I4-II=0';
eq6='I4+I6-I1=0';

% symbolic solver
[I1,I2,I3,I4,I6,Uj]=solve(eq1,eq2,eq3,eq4,eq5,eq6,I1,I2,I3,I4,I6,Uj)

% get numbers from symbols
I1_var=subs(I1)
I2_var=subs(I2)
I3_var=subs(I3)
I4_var=subs(I4)
I6_var=subs(I6)
Uj_var=subs(Uj)

% power balance
I1_abs=abs(I1_var)
I2_abs=abs(I2_var)
I3_abs=abs(I3_var)
I4_abs=abs(I4_var)
I6_abs=abs(I6_var)
Uj_abs=abs(Uj_var)

Sodd = U1*conj(I1_var)+Uj_var*conj(II)
Pod=real(Sodd)
Qod=imag(Sodd)

% final verification
Pp = ((I1_abs^2)*R1)+((I3_abs^2)*R2)+((abs(II)^2)*R4)+((I4_abs^2)*R3)
Qb = ((I2_abs^2)*XL1)+((I6_abs^2)*XL2)+((I1_abs^2)*XC1)

 
I1 =
 
(R2*U1 + R3*U1 + U1*XL1 + U1*XL2 - II*R2*XL2 + II*R3*XL1)/(R1*R2 + R1*R3 + R2*XC1 + R3*XC1 + R1*XL1 + R1*XL2 + R2*XL1 + R2*XL2 + R3*XL1 + R3*XL2 + XC1*XL1 + XC1*XL2)
 
 
I2 =
 
-(II*R1*R3 - R3*U1 - R2*U1 + II*R3*XC1 + II*R1*XL2 + II*R2*XL2 + II*R3*XL2 + II*XC1*XL2)/(R1*R2 + R1*R3 + R2*XC1 + R3*XC1 + R1*XL1 + R1*XL2 + R2*XL1 + R2*XL2 + R3*XL1 + R3*XL2 + XC1*XL1 + XC1*XL2)
 
 
I3 =
 
(U1*XL1 + U1*XL2 + II*R1*R3 + II*R3*XC1 + II*R1*XL2 + II*R3*XL1 + II*R3*XL2 + II*XC1*XL2)/(R1*R2 + R1*R3 + R2*XC1 + R3*XC1 + R1*XL1 + R1*XL2 + R2*XL1 + R2*XL2 + R3*XL1 + R3*XL2 + XC1*XL1 + XC1*XL2)
 
 
I4 =
 
-(II*R1*R2 - U1*XL2 - U1*XL1 + II*R2*XC1 + II*R1*XL1 + II*R2*XL1 + II*R2*XL2 + II*XC1*XL1)/(R1*R2 + R1*R3 + R2*XC1 + R3*XC1 + R1*XL1 + R1*XL2 + R2*XL1 + R2*XL2 + R3*XL1 + R3*XL2 + XC1*XL1 + XC1*XL2)
 
 
I6 =
 
(R2*U1 + R3*U1 + II*R1*R2 + II*R2*XC1 + II*R1*XL1 + II*R2*XL1 + II*R3*XL1 + II*XC1*XL1)/(R1*R2 + R1*R3 + R2*XC1 + R3*XC1 + R1*XL1 + R1*XL2 + R2*XL1 + R2*XL2 + R3*XL1 + R3*XL2 + XC1*XL1 + XC1*XL2)
 
 
Uj =
 
(R2*U1*XL2 - R3*U1*XL1 + II*R1*R2*R3 + II*R1*R2*R4 + II*R1*R3*R4 + II*R2*R3*XC1 + II*R2*R4*XC1 + II*R3*R4*XC1 + II*R1*R2*XL2 + II*R1*R3*XL1 + II*R1*R4*XL1 + II*R2*R3*XL1 + II*R1*R4*XL2 + II*R2*R3*XL2 + II*R2*R4*XL1 + II*R2*R4*XL2 + II*R3*R4*XL1 + II*R3*R4*XL2 + II*R2*XC1*XL2 + II*R3*XC1*XL1 + II*R4*XC1*XL1 + II*R4*XC1*XL2 + II*R1*XL1*XL2 + II*R2*XL1*XL2 + II*R3*XL1*XL2 + II*XC1*XL1*XL2)/(R1*R2 + R1*R3 + R2*XC1 + R3*XC1 + R1*XL1 + R1*XL2 + R2*XL1 + R2*XL2 + R3*XL1 + R3*XL2 + XC1*XL1 + XC1*XL2)
 

I1_var =

   2.5805 + 1.2511i


I2_var =

   2.9553 - 1.0763i


I3_var =

  -0.3748 + 2.3274i


I4_var =

   0.9590 + 0.8372i


I6_var =

   1.6214 + 0.4139i


Uj_var =

 -48.6856 +70.0749i


I1_abs =

    2.8678


I2_abs =

    3.1452


I3_abs =

    2.3574


I4_abs =

    1.2731


I6_abs =

    1.6734


Uj_abs =

   85.3275


Sodd =

  4.2741e+002 -1.4603e+002i


Pod =

  427.4136


Qod =

 -146.0326


Pp =

  427.4136


Qb =

        0 -1.4603e+002i

Published with MATLAB®