How to sort eigenvalues when performing generalized nyquist criterion on 2x2 matrix?

Hi everyone,

I am currently using the GNC method to determine the stability of my system.

However, the function eig() of MATLAB switches the order of the eigenvalues. Thus, my plots are showing wrong behaviour. I have tried to sort them using sort() and also using eigs(). A plot showing the error is attached as pdf.

I would be very grateful for any lead on this problem!

Thanks in advance,

Merve

Answer 1

[Bruno Luong]

Replace statement

E1(:,k) = eig(Z_nyquist1);

by

    % track eigen values -----------------------------------------------
    Ek = eig(Z_nyquist1);
    if k > 1
        Ekp = E1(:,k-1);
        Ekf = flip(Ek);
        if norm(Ek-Ekp, 'Inf') > norm(Ekf-Ekp, 'Inf')
            Ek = Ekf;
        end
    end
    E1(:,k) = Ek;
    % track eigen values ----------------------------------------------- 
    

[William Rose]

@Merve Can,

@Bruno Luong is smart! He compares the current vector of eigenvalues to the vector of eigenvalues at the previous frequency. If the flipped current one is closer to the prior than the unflipped, used the flipped. Otherwise, use the unflipped. Nice!

Answer 2

[William Rose]

@Merve Can,

Since you did not supply code, I am not exactly sure where the problem arises. I wonder if you are not sorting all the things you need to sort. If you use the eigenvectors, then you must sort them the same way you sort the eigenvalues. The Help for eig() explains this nicely.

A=rand(3,3);
[V,D]=eig(A);           %unsorted
[d,ind]=sort(diag(D));  %sort the eigenvalues, ascending
Ds=D(ind,ind);          %sorted eigenvalue matrix
Vs=V(:,ind);            %sorted eigenvector matrix

Try it. Good luck.

Answer 3

[Merve Can]

Hi @William Rose,

the code is rather long and I don't know how to show the problem in a better way. I have now attached it.

In line 224 I am calculating the eigenvalues to generate the GNC plot for a specified frequency range. I have tried the sorting code but it did not work for me. Could you have a look?

Thanks a lot!

Merve

Answer 4

[William Rose]

@Merve Can

I downloaded your code and ran it. It runs fine. The image generated is too large to fit on my small screen.

You compute eigenvalues and eigenvectors every time through the loop over k. This loop has 100036 passes. That is a lot of eigenvalue calculations! No wonder the code takes a little while to run.

The relavant code is

Z_nyquist1      =  - Z_PFQV/Z_grid;
E1(:,k) = eig(Z_nyquist1);

where z_PFQV and Z_grid and Z_nyquist1 are each complex 2x2 matrices. I assume that these matrices, and the use of matrix right division as you are doing, is correct for this situation. I do not know if that is true but I assume it is. I notice that Z_nyquist, at least for the last loop pass, is NOT Hermitian:

>> disp(Z_nyquist1*1e3)
 -10.8234 - 0.2136i  -0.0034 + 0.1016i
   0.0034 - 0.1016i -10.8234 - 0.2136i

Since Z_nyquist1 is not Hermitian, the eigenvalues are not real - they are complex. Is that what you expect? I am not familiar with the GNC method, so I don't know. How do you want to sort the complex eigenvalues? By magnitude, by the real part, by the imaginary part? sort(compex vector) sorts by the magnitude. I have tried 4 ways: no sort, sort by magnitude, sort by real, sort by imaginary. See below. The plots are different each time. Are any of them right? I don't know, since I don't know what they're supposed to look like.

Unsorted, sort by magntitude:

Sort by real part, sort by imaginary part:

See lines 225-228 of attached code.

[Merve Can]

Hi @William Rose,

thank you for your answer!

Yes the code is correct and I have verified it with a time domain simulation. Unfortunately, the impedance needs to be recalculated for every step.

All four plots are not correct, each eigenvalue draws is supposed to draw its own contour.

Best,

Merve

[Merve Can]

Hi @Bruno Luong,

thanks a lot for this code!

I have tried to compare the vectors but my code did not work as well as yours.

Cheers!

Merve

[Bruno Luong]

Hi @Merve Can

"my code did not work as well as yours."

if you reply intended to me (your comment is on @William Rose answer)

Really? result run on TMW server (it looks fine to me):

%clc; %not a fan of this
%clear;
%% Simulation time-step
Sref    = 1000;
% Sampling and switching
sf_f    = 1;
fs      = 3200*sf_f; % switching frequency
% Network
fn      = 50;
wn      = 2*pi*fn;
%% Power scaling
sf_i    = sqrt(0.15/(5.6*250));%6/500;  % current scaling factor
sf_v    = 1/(250*sf_i);%1/3;    % voltage scaling factor
Roni    = 1e-3 / sf_i * sf_v;
%% Circuit parameters
% LC filter
Lf_1    = 5.6e-3;
Rf_1    = 0.1;
Cf_1    = 600e-6*sf_i/sf_v/sf_f;
Rc      = 0;
% transmission line(grid coupling impedance)
SCR     = 32;
Zg      = (3 * (181*sf_v) * (181*sf_v)) / (SCR * Sref);
L_1     = Zg / wn;
R_1     = 0.053;
% Virtual impedance
Lv_1    = (5.6e-3)*2;
Rv_1    = 0.1;
%% Control parameters
% power measurement LPF
wc      = 20*pi;
wz_vir  = 500*pi; %differentiation filter
% inner current loop
Tsf     = 1/fs;
T_i     = 2*Tsf;
Kpc_1   = Lf_1 / T_i;
Kic_1   = (Rf_1+Roni)/T_i;
% outer voltage loop
wp      = 1/(T_i);
% Phase margin
dm      = 60*pi/180;      %phase margin, typical values between 50 and 80 degrees
Ti_u    = 1/((1-sin(dm))/(T_i+T_i*sin(dm))); % time constant of the controller
% cross frequency
wcv     = sqrt(1/(T_i*Ti_u));                 % controller bandwidth
Kpv_1   = wcv*Cf_1;
Kiv_1   = Kpv_1/Ti_u;
F_i     = 0.75;
% ?? reference voltage
Vg      = 181*sqrt(2)*sf_v;
Udc     = 800*sf_v;
% Droop Control
mp      = 0.01 / 1000 * 2 * pi * 50;
nq      = 0.005 * Vg / Sref;
LPF     = 15;
%% Steady State Values
% At least 5 decimal places
% mp = 0.01, nq = 0.005, wc = 20*pi,wz_vir = 500*pi;
% Controller Reference Frame
Ud0_c   = 95.92601;
Uq0_c   = -23.74511;
Id0_c   = 6.76003;
Iq0_c   = -0.76663;
% Theta Difference
theta_0 = 0.27442;
cos0    = cos(theta_0);
sin0    = sin(theta_0);
T_theta = [cos0 sin0;-sin0 cos0];
% Grid reference frame
Ud0_s = 98.77137;
Uq0_s = 3.13789;
Id0_s = 6.71484;
Iq0_s = 1.09392;
%% Impedance Model - MC
f_sweep1                 = 1:0.1:10000;
theta_sweep             = 0:0.01*pi:pi/2;
f_sweep2                = 1*exp(1i*theta_sweep)/(1i);
f_sweep                 = [f_sweep2 f_sweep1];
f_div                   = abs(f_sweep-50);
index                   = find(f_div<0.3);
f_sweep(index)          = [];
no_of_f                 = length(f_sweep);
E1                      = zeros(2,no_of_f);
s                       = 1i*2*pi*f_sweep;
for k = 1:no_of_f
    G_ov        = [Rv_1 -wn*Lv_1; wn*Lv_1 Rv_1];
    G_ffv       = [F_i 0;0 F_i];
    G_delv      = [0 -wn*Cf_1;wn*Cf_1 0];
    G_PIi       = [(Kpc_1 + Kic_1 * 1/s(k)) 0; 0 (Kpc_1 + Kic_1 * 1/s(k))];
    G_deli      = [0 -wn*Lf_1;wn*Lf_1 0];
    G_PIv       = [(Kpv_1 + Kiv_1 * 1/s(k)) 0; 0 (Kpv_1 + Kiv_1 * 1/s(k))];
    Z_l         = [s(k)*Lf_1+Rf_1 -wn*Lf_1;wn*Lf_1 s(k)*Lf_1+Rf_1];
    Z_c         = [s(k)*Cf_1 -wn*Cf_1;wn*Cf_1 s(k)*Cf_1];
    G_vodq      = Z_l*Z_c +[1 0;0 1] + G_PIi*G_PIv ...
        - G_PIi*G_delv + G_PIi*Z_c ...
        - G_deli*Z_c;
    G_iodq      =  -G_PIi*G_PIv*G_ov + G_PIi*G_ffv ...
        - G_PIi + G_deli - Z_l;
    Z_innerloop = G_iodq/G_vodq;
    G_voqdref   = G_PIi*G_PIv;
    % PF Droop Impedance
    Icss        = 3/2*[Id0_c Iq0_c];
    Vcss        = 3/2*[Ud0_c Uq0_c];
    Tu          = [-sin0*Ud0_s+cos0*Uq0_s;-sin0*Uq0_s-cos0*Ud0_s];
    Ti          = [-sin0*Id0_s+cos0*Iq0_s;-sin0*Iq0_s-cos0*Id0_s];
    G_LPF       = LPF / (s(k) + LPF);
    % X_Theta step by step derivation
    X_theta_1   = (1/s(k))*mp*G_LPF*(Icss*Tu+Vcss*Ti);
    X_theta     = 1+X_theta_1;
    % G_vodqs step by step derivation
    G_vodqS_1   = G_vodq*T_theta;
    G_vodqS_2   = G_vodq*Tu*mp*G_LPF*Icss*T_theta*1/(X_theta*s(k));
    G_vodqS_3   = G_iodq*Ti*mp*G_LPF*Icss*T_theta*1/(X_theta*s(k));
    G_vodqS     = G_vodqS_1 - G_vodqS_2 + G_vodqS_3;
    % G_iodqs step by step derivation
    G_iodqS_1   = G_iodq * T_theta;
    G_iodqS_2   = G_iodq*Ti*mp*G_LPF*Vcss*T_theta*1/(X_theta*s(k));
    G_iodqS_3   = G_vodq*Tu*mp*G_LPF*Vcss*T_theta*1/(X_theta*s(k));
    G_iodqS     = G_iodqS_1 - G_iodqS_2 + G_iodqS_3;
    Z_2_pv         = G_vodqS\G_iodqS;
    Z_1_pv         = G_vodqS\G_voqdref;
    %--------------------------------------------------------------------------
    % % QV Droop Impedance
    IcQss = 3/2*[-Iq0_c Id0_c;0 0];
    VcQss = 3/2*[Uq0_c -Ud0_c;0 0];
    G_Qvodq_1   = (mp*Tu*G_LPF*Icss*T_theta)/(X_theta*s(k));
    G_Qvodq_2   = IcQss*(T_theta - G_Qvodq_1);
    G_Qvodq_3   = VcQss*(mp*Ti*G_LPF*Icss*T_theta)/(X_theta*s(k));
    G_Qvodq     = -nq*G_LPF*(G_Qvodq_2 - G_Qvodq_3);
    %--------------------------------------------------------------------------
    G_Qiodq_1   = (mp*Ti*G_LPF*Vcss*T_theta)/(X_theta*s(k));
    G_Qiodq_2   = VcQss*(T_theta - G_Qvodq_1);
    G_Qiodq_3   = IcQss*(mp*Tu*G_LPF*Vcss*T_theta)/(X_theta*s(k));
    G_Qiodq     = -nq*G_LPF*(G_Qiodq_2-G_Qiodq_3);
    % % % % % % % % % % Whole System Impedance % % % % % % % % % %
    Z_PFQV_1       = G_PIi*G_PIv*G_Qiodq;
    Z_PFQV_2       = G_iodqS + Z_PFQV_1;
    Z_PFQV_3       = G_PIi*G_PIv*G_Qvodq;
    Z_PFQV_4       = G_vodqS - Z_PFQV_3;
    Z_PFQV         = Z_PFQV_4\Z_PFQV_2;
    % Grid Impedance
    Z_grid   = [s(k)*L_1+R_1 -wn*L_1;wn*L_1 s(k)*L_1+R_1];
    % Minor Loop Gain defined in the literature for grid forming converter
    Z_nyquist1      =  - Z_PFQV/Z_grid;
    %E1(:,k) = eig(Z_nyquist1);         %no sort
    % track eigen values -----------------------------------------------
    Ek = eig(Z_nyquist1);
    if k > 1
        Ekp = E1(:,k-1);
        Ekf = flip(Ek);
        if norm(Ek-Ekp, 'Inf') > norm(Ekf-Ekp, 'Inf')
            Ek = Ekf;
        end
    end
    E1(:,k) = Ek;
    % track eigen values ----------------------------------------------- 
    
end
%% Unity Circle
theta = 0 : 0.01 : 2*pi;
radius =1;
x = radius * cos(theta);
y = radius * sin(theta);
%%
h=figure;
Fontsize = 24;
%set(h,'position',[100 100 900 700]); %not a fan of this
plot(real(E1(2,:)),imag(E1(2,:)),'m')
set(gca,'FontSize',Fontsize)
hold on
plot(real(E1(1,:)),imag(E1(1,:)),'b')
hold on
plot(real(E1(2,:)),-imag(E1(2,:)),'m')
hold on
plot(real(E1(1,:)),-imag(E1(1,:)),'b')
hold on
plot(-1,0,'kx','MarkerSize',10,'LineWidth',1)
ax = gca;
ax.TickLabelInterpreter = 'latex';
grid on
title('Nyquist Plot','interpreter','latex');
xlabel('Real Axis','interpreter','latex')
ylabel('Imaginary Axis','interpreter','latex')
legend('$\lambda_\mathrm{1}$','$\lambda_\mathrm{2}$','interpreter','latex')
axes('position',[.65 .2 .25 .25])
box on % put box around new pair of axes
plot(real(E1(2,:)),imag(E1(2,:)),'m')
set(gca,'FontSize',Fontsize)
xlim([-2 1])
ylim([-2 2])
hold on
plot(real(E1(2,:)),-imag(E1(2,:)),'m')
xlim([-2 1])
ylim([-2 2])
hold on
plot(real(E1(1,:)),imag(E1(1,:)),'b')
xlim([-2 1])
ylim([-2 2])
grid on
hold on
plot(real(E1(1,:)),-imag(E1(1,:)),'b')
xlim([-2 1])
ylim([-2 2])
hold on
plot(-1,0,'kx')
hold on
plot(x,y,'b','MarkerSize',10,'Linewidth',1)
ax = gca;
ax.TickLabelInterpreter = 'latex';

[Merve Can]

Hi,

I meant my own code did not work as well as yours. Sorry!