How to solve singular matrix?

function dx = plant(t,x,u)
dx = zeros(6,1);
x1 = [x(1);x(2);x(3)];
x2 = [x(4);x(5);x(6)];
beta3=0.33;
beta4=2;
D=0.1;
M=1;
a1=0.1;%k1
a2=0.1;%k2
a3=0.01;%k3
a4=0.1;%k4
b=2.3;
f2= [(-(D*a1+a3)*x1(1)/M)-(D*a2*b^3*x2(3)/M);1;2];
G=[(1+a4*b^3)/M 0 0;0 1 0;0 0 1];
p1=[beta3*sin(x1(1));1;2];
p2=[beta4*(x2(1)^2);0;3];
dx1 = x2  + p1;
dx2 = f2  + G*u  + p2;
dx = [dx1; dx2];

Answer 1

[Sam Chak]

Hi @YING CHIH

The Plant system itself does not have any division-by-zero singularity issues. However, the simulation is terminated due to the state variable number 4 () experiencing an exponential growth within a finite time. This is likely caused by the inability of the External Input u to stabilize the Plant system. To resolve this issue, you simply need to design a stabilizing input vector signal.

uExtIn  = ones(3, 1);       % External Input signal
tspan   = [0 0.765];
x0      = [1; zeros(5,1)];
[t, x]  = ode45(@(t, x) plant(t, x, uExtIn), tspan, x0);
plot(t, x), grid on, legend('show', 'location', 'nw')

function dx = plant(t, x, u)
    dx      = zeros(6,1);
    x1      = [x(1);
               x(2);
               x(3)];
    x2      = [x(4);
               x(5);
               x(6)];
    beta3   = 0.33;
    beta4   = 2;
    D       = 0.1;
    M       = 1;
    a1      = 0.1;  % k1
    a2      = 0.1;  % k2
    a3      = 0.01; % k3
    a4      = 0.1;  % k4
    b       = 2.3;
    f2      = [(-(D*a1+a3)*x1(1)/M)-(D*a2*b^3*x2(3)/M);
               1;
               2];
    G       = [(1+a4*b^3)/M 0 0;
               0            1 0;
               0            0 1];
    p1      = [beta3*sin(x1(1));
               1;
               2];
    p2      = [beta4*(x2(1)^2);
               0;
               3];
    
    dx1     = x2 + p1;
    dx2     = f2 + G*u + p2;
    dx      = [dx1; 
               dx2];
end

[Sam Chak]

To illustrate my point, if a stabilizing input vector signal is designed using the principles of linear algebra and satisfies the stability condition, then the system responses will be stable.

uExtIn  = [];           % No external Input signal
tspan   = [0 10];
x0      = [1; zeros(5,1)];
[t, x]  = ode45(@(t, x) plant(t, x, uExtIn), tspan, x0);
plot(t, x), grid on, legend('show', 'location', 'nw')

function dx = plant(t, x, u)
    dx      = zeros(6,1);
    x1      = [x(1);
               x(2);
               x(3)];
    x2      = [x(4);
               x(5);
               x(6)];
    beta3   = 0.33;
    beta4   = 2;
    D       = 0.1;
    M       = 1;
    a1      = 0.1;  % k1
    a2      = 0.1;  % k2
    a3      = 0.01; % k3
    a4      = 0.1;  % k4
    b       = 2.3;
    f2      = [(-(D*a1+a3)*x1(1)/M)-(D*a2*b^3*x2(3)/M);
               1;
               2];
    G       = [(1+a4*b^3)/M 0 0;
               0            1 0;
               0            0 1];
    p1      = [beta3*sin(x1(1));
               1;
               2];
    p2      = [beta4*(x2(1)^2);
               0;
               3];
    
    % ----------------------------------
    A       = [zeros(3), eye(3); zeros(3, 6)];
    B       = [zeros(3); G];
    K       = lqr(A, B, eye(6), eye(3));
    uLQR    = - K*x;
    % ----------------------------------
    
    dx1     = x2 + p1;
    dx2     = f2 + G*uLQR + p2;
    dx      = [dx1; 
               dx2];
end

[YING CHIH]

I tried using lqr but encountered "Function 'lqr' not supported for code generation." However, because it was a controller issue, I tried to change u, which is my controller u, and after making the changes, it will be the same as the initial problem.

function u_false  = fcn(t,u)
u_false=zeros(3,1);
if t<1.5
    ut1=0;
else
    ut1=0.6;
end
if t<0.5
    ut2=0;
else
    ut2=0.9;
end
if t<1
    ut3=0;
else
    ut3=0.2;
end
tran1 = [ut1  0 0; 0 ut2 0;0 0 ut3];

u_false = (eye(3)-tran1)*u;

[Sam Chak]

Since I was not aware of the specifics of your u signal, the implementation of lqr in the code was merely my demonstration for automatically designing a stabilizing input signal u in MATLAB.

Based on your code, it appears that your controllers are If-Else discontinuous functions of time. If you are encountering a persistent error message, it suggests that your discontinuous signal u is causing destabilization of the system.

As long as your input signal u is stabilizing, the error message will disappear. A better way is to manually design u that 100% guarantees the stability.

Answer 2

[Sam Chak]

Hi @YING CHIH

Hear is another stabilizing input vector signal.

uExtIn  = [];           % No external Input signal
tspan   = [0 20];
x0      = [1; zeros(5,1)];
[t, x]  = ode45(@(t, x) plant(t, x, uExtIn), tspan, x0);
subplot(211)
plot(t, x(:,1:3)), grid on, legend('x_1', 'x_2', 'x_3', 'location', 'ne')
title('Responses of X1 vector')
subplot(212)
plot(t, x(:,4:6)), grid on, legend('x_4', 'x_5', 'x_6', 'location', 'ne')
title('Responses of X2 vector')

function dx = plant(t, x, u)
    dx      = zeros(6,1);
    x1      = [x(1);
               x(2);
               x(3)];
    x2      = [x(4);
               x(5);
               x(6)];
    beta3   = 0.33;
    beta4   = 2;
    D       = 0.1;
    M       = 1;
    a1      = 0.1;  % k1
    a2      = 0.1;  % k2
    a3      = 0.01; % k3
    a4      = 0.1;  % k4
    b       = 2.3;
    f2      = [(- (D*a1 + a3)*x1(1)/M) - (D*a2*b^3*x2(3)/M);
               1;
               2];
    G       = [(1+a4*b^3)/M 0 0;
               0            1 0;
               0            0 1];
    p1      = [beta3*sin(x1(1));
               1;
               2];
    p2      = [beta4*(x2(1)^2);
               0;
               3];
    
    % ----------------------------------
    k1      = 1;
    k2      = 1;
    x2d     = - k1*x1 - p1;
    dx2d    = - k2*(x2 - x2d);
    u       = G\(dx2d - p2 - f2);
    % ----------------------------------
    
    dx1     = x2 + p1;
    dx2     = f2 + G*u + p2;
    dx      = [dx1; 
               dx2];
end