How to calculate minimum singular value of I+L and I+inv(L) for MIMO system?

Hello,

I want to calculate minimum singular value for MIMO system, and here is an example

%%
A = [0.1873	 0.8403	-0.6003	-2.1384	 0.1240;
    -0.0824	-0.8880	 0.4899	-0.8395  1.4366;
    -1.9330  0.1001  0.7394  1.3546 -1.9609;
    -0.4390 -0.5445  1.7119 -1.0722 -0.1977;
    -1.7947  0.3035 -0.1941  0.9610 -1.2078];
B = [ 2.9080 -0.2725;
    0.8252  1.0984;
    1.3790 -0.2779;
    -1.0582  0.7015;
    -0.4686 -2.0518];
C = [1 0 0 0 0; 0 0 1 0 0];
D = zeros(2,1);
%%
N = 100;
a = logspace(-2,2,N);
Q = diag([1 1 10 10 2]);
R = eye(2);
w = logspace(-3,3,1000);
for i = 1:N
    Qi = a(i)*Q;
    K = lqr(A,B,Qi,R);
    [re,im] = nyquist(A,B,K,0,1,w);
    l = (re+sqrt(-1)*im);
    sr = ones(size(l)) + ones(size(l))./l;
    rd = ones(size(l)) + l;
    srmin2(i,:) = min(abs(sr));
    rdmin2(i,:) = min(abs(rd));
end
plot(a, srmin2);

Is it correct code for calculation of minimum singular value of I+L and I+inv(L) for MIMO system?

and if it is correct, does the plot means minimum singular value for each input?

Answer 1

[Ayush]

Hello Matt,

Based on the code provided by you for calculating minimum singular value of I+L and I+inv(L) for a MIMO system, I can figure out some issues that might result in a different outcome than what you are intending.

Here are some areas in your code that can be looked at for proper calculation and results for the MIMO system.

1. The “nyquist” function is incorrectly used with parameters (A,B,K,0,1,w) because in MATLAB this function does not accept state-space matrices and thus requires a system model object created by the “ss” function. Please refer to the below documentations to know more about the “nyquist” and “ss” function respectively:

• https://www.mathworks.com/help/control/ref/dynamicsystem.nyquist.html
• https://www.mathworks.com/help/control/ref/ss.html

2. The computation of the transfer function (L) was also incorrect, since the “nyquist” function returns a complex frequency which include a real part (re), imaginary part (im) and a vector of frequencies (wout). The code does not extract the frequency required for the MIMO system and should have used the “squeeze” function to remove dimensions of length 1.

• https://www.mathworks.com/help/matlab/ref/squeeze.html

3. The calculation of the minimum singular value is not as desired since the use of “min” function in the above code for min(abs(sr)) and min(abs(rd)) will give minimum over all frequencies and all singular values, not the minimum singular values at each frequency. Hence, the “svd” function can be used to decompose into the singular value and then find its minimum. Please refer to documentation to know more about the “svd” function:

• https://www.mathworks.com/help/matlab/ref/double.svd.html

You can also refer to the below example code snippet, that incorporates the above mentioned changes for the issues found:

A = [0.1873 0.8403 -0.6003 -2.1384 0.1240;
    -0.0824 -0.8880 0.4899 -0.8395 1.4366;
    -1.9330 0.1001 0.7394 1.3546 -1.9609;
    -0.4390 -0.5445 1.7119 -1.0722 -0.1977;
    -1.7947 0.3035 -0.1941 0.9610 -1.2078];
B = [2.9080 -0.2725;
    0.8252 1.0984;
    1.3790 -0.2779;
    -1.0582 0.7015;
    -0.4686 -2.0518];
C = [1 0 0 0 0; 0 0 1 0 0];
D = zeros(2, 2); % D should be 2x2 because B has 2 columns

N = 100;
a = logspace(-2, 2, N);
Q = diag([1 1 10 10 2]);
R = eye(2);
w = logspace(-3, 3, 1000);
srmin2 = zeros(N, length(w));
rdmin2 = zeros(N, length(w));

for i = 1:N
    Qi = a(i) * Q;
    K = lqr(A, B, Qi, R);
    sys = ss(A - B * K, B, C, D); 

    for j = 1:length(w)
        [re, im, ~] = nyquist(sys, w(j));
        L = squeeze(re(1,1,:)) + 1i * squeeze(im(1,1,:));

        % Calculate singular values of I+L and I+inv(L)
        sv_IL = svd(eye(size(L)) + L);
        sv_I_invL = svd(eye(size(L)) + inv(L));

        % Store the minimum singular values
        srmin2(i, j) = min(sv_IL);
        rdmin2(i, j) = min(sv_I_invL);
    end
end

% Plot the minimum singular value across all frequencies for each alpha
[min_srmin2, ~] = min(srmin2, [], 2); % Minimum over frequencies
[min_rdmin2, ~] = min(rdmin2, [], 2); % Minimum over frequencies

plot(a, min_srmin2);
xlabel('Alpha');
ylabel('Minimum Singular Value of I+L');
title('Minimum Singular Value of I+L vs Alpha');

Hope it helps!