Why does transforming transfer matrix into complex vector transfer function lead to NaN ?

Hi,

I am currently trying to transform my transfer matrix into compex vector form. The following picture shows how it is theoretically done.

My is Z_2 and has quite a high order. So when I try to use eq. (F2) to get Matlab returns Z_2_cmplx = NaN. Here is my code:

%% Calculation of Transfer Function

% Parameters

Lf_1 = 5.6e-3;

Rf_1 = 0.1;

Cf_1 = 1.61e-05;

L_1 = 1.46e-3 ;

R_1 = 0.053;

wn = 2*pi*50;

Kpc_1 = 35.839999999999996;

Kic_1 = 8.789333333333335e+02;

Kpv_1 = 0.027560488364344;

Kiv_1 = 12.664025834171640;

F_i = 0.9;

% Calculating Transfer Functions

s = tf('s');

Z_PIi = (Kpc_1 + Kic_1 * 1/s) ;

Z_CDi = (-wn*Lf_1) * [0 -1;1 0];

Z_inner = (Z_PIi + (s*Lf_1 + Rf_1))*[1 0;0 1] + (wn*Lf_1)*[0 -1;1 0] + Z_CDi;

Z_inner_inv = inv(Z_inner);

G_I = Z_PIi / Z_inner;

Z_PIv = 1/(G_I * ((Kpv_1 + Kiv_1 * 1/s) * [1 0;0 1]));

Z_PIv_inv = inv(Z_PIv); %Removes extra s

Z_CDv = -1/((wn*Cf_1*[0 -1;1 0])*G_I);

Z_CDv_inv = inv(Z_CDv);

Z_parallel_inv = Z_PIv_inv + Z_CDv_inv + [s*Cf_1 -wn*Cf_1;wn*Cf_1 s*Cf_1] + Z_inner_inv; %equal to Zinner//Zpar mentioned in the notes

Z_parallel = 1/Z_parallel_inv;

Z_1 = Z_parallel/Z_PIv;

Z_ov = [R_1, -wn*Lv_1;

wn*Lv_1, R_1];

Z_Fi = -Z_parallel / inv(G_I)*F_i;

%% Transformation into complex vector form

Z_2 = Z_parallel + Z_Fi + (Z_ov * Z_1);

Z_2_cmplx = (Z_2(1,1)+Z_2(2,2))/2 + 1i * (Z_2(2,1)+Z_2(1,2))/2 % This yields NaN :/

Can anyone please help me with this? I am really out of ideas on how to solve this issue.

Thanks a lot!!!!

Merve

[Paul]

Can't run the code because, variable Lv_1 is undefined.

Also, the expression for Z_2_cmplx doesn't follow the equation in the picture. It should be

Z_2_cmplx = (Z_2(1,1)+Z_2(2,2))/2 + 1i * (Z_2(2,1) - Z_2(1,2))/2  

Because Z_2_cmplx is supposed to be G_dq_plus (I think).

[Merve Can]

Hi Paul,

thank you for your answer. Yes, sorry it should be a “-“.Nevertheless, I am getting NaN.

 Lv_1 = (5.6e-3)*3;

Can you run the code now?

Answer 1

[Paul]

Recreating the result

%% Calculation of Transfer Function
% Parameters
Lf_1  = 5.6e-3;
Rf_1  = 0.1;
Cf_1  = 1.61e-05;
L_1   = 1.46e-3 ;
R_1   = 0.053;
Lv_1 = (5.6e-3)*3;
wn    = 2*pi*50;
Kpc_1 = 35.839999999999996;
Kic_1 = 8.789333333333335e+02;
Kpv_1 = 0.027560488364344;
Kiv_1 = 12.664025834171640;
F_i   = 0.9; 
% Calculating Transfer Functions
s = tf('s');  
Z_PIi               = (Kpc_1 + Kic_1 * 1/s) ;
Z_CDi               = (-wn*Lf_1) * [0 -1;1 0];
Z_inner             = (Z_PIi + (s*Lf_1 + Rf_1))*[1 0;0 1] + (wn*Lf_1)*[0 -1;1 0] + Z_CDi;
Z_inner_inv         = inv(Z_inner);
G_I                 = Z_PIi / Z_inner; 
Z_PIv               = 1/(G_I * ((Kpv_1 + Kiv_1 * 1/s) * [1 0;0 1]));
Z_PIv_inv           = inv(Z_PIv);                                                               %Removes extra s
Z_CDv               = -1/((wn*Cf_1*[0 -1;1 0])*G_I);
Z_CDv_inv           = inv(Z_CDv);
Z_parallel_inv      = Z_PIv_inv + Z_CDv_inv  + [s*Cf_1 -wn*Cf_1;wn*Cf_1 s*Cf_1] + Z_inner_inv;  %equal to Zinner//Zpar mentioned in the notes
Z_parallel          = 1/Z_parallel_inv;
Z_1                 = Z_parallel/Z_PIv;
Z_ov                = [R_1, -wn*Lv_1;
                       wn*Lv_1, R_1];
Z_Fi                = -Z_parallel / inv(G_I)*F_i;
Z_2                 = Z_parallel + Z_Fi + (Z_ov * Z_1);
Z_2_cmplx           = (Z_2(1,1)+Z_2(2,2))/2 + 1i * (Z_2(2,1)-Z_2(1,2))/2 % This yields NaN :/

Z_2_cmplx =
 
  NaN
 
Static gain.

Look at the coefficients of the numerator and denominator of Z_2(1,1). They are enormous, here's the five largest. Same thing with Z_2(2,2)

[num,den] = tfdata(Z_2(1,1));
sortnum = sort(num{:});
sortden = sort(num{:});
sortnum(end-5:end)

ans = 1×6

1.0e+186 * 0.0271 0.1262 0.4510 1.1589 1.4996 1.9038

sortden(end-5:end)

ans = 1×6

1.0e+186 * 0.0271 0.1262 0.4510 1.1589 1.4996 1.9038

When Z_2(1,1) and Z_2(2,2) are combined, problems result with the coefficients getting too large.

[num,den] = tfdata(Z_2(1,1)+Z_2(2,2));
any(isinf(num{:}))

ans = logical
   1

any(isinf(den{:}))

ans = logical
   1

There may be a better way to do all that transfer function manipulation. Are those equations derived from a block diagram?

[Merve Can]

Yes it’s a method published in a paper that I am supposed to implement. This is the title „ Impedance Circuit Model of Grid-Forming Inverter: Visualizing Control Algorithms as Circuit Elements“.

I am out of ideas on how to go on from here. Do you have an idea?

[Paul]

Without looking at the paper, I tried to mimic the original equations, but sticking with ss form wherever possible and finding minimal realizations at each step.

Please check carefully that this code is functionally equivalent to your code.

I didn't take the last step of combining the real and imaginary parts of Z2_cmplx, because, frankly, I've never seen anything where an LTI system is multiplied by 1i

Hopefully you have some idea of what the result should look like and compare to the result achieved below.

%% Calculation of Transfer Function
% Parameters
Lf_1  = 5.6e-3;
Rf_1  = 0.1;
Cf_1  = 1.61e-05;
L_1   = 1.46e-3 ;
R_1   = 0.053;
Lv_1 = (5.6e-3)*3;
wn    = 2*pi*50;
Kpc_1 = 35.839999999999996;
Kic_1 = 8.789333333333335e+02;
Kpv_1 = 0.027560488364344;
Kiv_1 = 12.664025834171640;
F_i   = 0.9; 
% Calculating Transfer Functions
s = tf('s');  
Z_PIi               = ss(Kpc_1 + Kic_1 * 1/s) ;
Z_CDi               = ss((-wn*Lf_1) * [0 -1;1 0]);
Z_inner             = (Z_PIi + ((s*Lf_1 + Rf_1))*[1 0;0 1]) + (wn*Lf_1)*[0 -1;1 0] + Z_CDi;
Z_inner_inv         = minreal(ss(zpk(inv(Z_inner))));

3 states removed.

% G_I                 = Z_PIi / Z_inner; 
G_I                 = minreal(ss(zpk(Z_PIi * Z_inner_inv))); 

5 states removed.

% Z_PIv               = 1/(G_I * ((Kpv_1 + Kiv_1 * 1/s) * [1 0;0 1]));
% Z_PIv_inv           = inv(Z_PIv);                                                               %Removes extra s
Z_PIv_inv           = minreal(ss(zpk(G_I * ((Kpv_1 + Kiv_1 * 1/s) * [1 0;0 1]))));

3 states removed.

% Z_CDv               = -1/((wn*Cf_1*[0 -1;1 0])*G_I);
% Z_CDv_inv           = inv(Z_CDv);
Z_CDv_inv           = -((wn*Cf_1*[0 -1;1 0])*G_I);

Z_parallel_inv      = Z_PIv_inv + Z_CDv_inv  + [s*Cf_1 -wn*Cf_1;wn*Cf_1 s*Cf_1] + Z_inner_inv;  %equal to Zinner//Zpar mentioned in the notes
% Z_parallel          = 1/Z_parallel_inv;
Z_parallel          = minreal(ss(zpk(inv(Z_parallel_inv))));

18 states removed.

% Z_1                 = Z_parallel/Z_PIv;
Z_1                 = minreal(Z_parallel*Z_PIv_inv);

3 states removed.

Z_ov                = [R_1, -wn*Lv_1;
                       wn*Lv_1, R_1];

% Z_Fi                = -Z_parallel / inv(G_I)*F_i;
Z_Fi                = minreal(-Z_parallel * G_I * F_i);

1 state removed.

Z_2                 = minreal(Z_parallel + Z_Fi + (Z_ov * Z_1));

10 states removed.

% Z_2_cmplx           = (Z_2(1,1)+Z_2(2,2))/2 + 1i * (Z_2(2,1)-Z_2(1,2))/2;
Z_2_cmplx_real        = minreal(Z_2(1,1) + Z_2(2,2))/2; 

26 states removed.

Z_2_cmplx_imag        = minreal(Z_2(2,1) - Z_2(1,2))/2;

26 states removed.

bode(Z_2_cmplx_real)

bode(Z_2_cmplx_imag)

[Merve Can]

Thank you for your answer! I will try to use the minreal( command. Unfortunately, I can't use the state space function because my aim is to build an impedance model to investigate the stability of my system.

[Paul]

Why can't a state space model be used? It's just a different representation of the same system, so should be usable for any anlaysis that needs to be done. I'd be very interested if there's some type of analysis that can't be done with a state space representation but can be done with some othe representation.

[Merve Can]

Maybe I can explain it in a better way in an e-mail. Could you please contact me via [e-mail]?

[Paul]

No, I don't do email exchanges.