How to solve the matrix P in the Lyapunov equation "A'P+PA=- Q" in Simulink

For the solution of matrix P in Lyapunov equation, I can call the ‘lyap’ function in the command line window or m file to solve. The method is as follows:

A=[0,1;-1,-2]';Q=[1,0;0,1];P=lyap(A,Q)

However, now I need to solve the matrix P in 'A'P+PA=- Q' in Simulink (matrices A and Q are known). At this point, the ‘lyap’ function cannot be called. Is there any good way to solve this problem? I look forward to your reply very much.

[Sam Chak]

Hi @LIULIYAN

Is there a practical reason why you want to solve the Lyapunov equation to find the matrix in the Simulink environment when the matrices and are fixed?

Unless the state matrix is time-varying, , and you want to prove the global asymptotic stability of a 2nd-order system via the Markus–Yamabe conjecture, I can't think of any 🤷‍♂️.

[LIULIYAN]

Hi@Sam Chak

In fact, matrix A is a Hurwitz matrix about the control gains Kp and Kv, which means it is variable. Q is a Positive-definite matrix. The control output of the controller I designed contains matrix P, so I need to solve it.

[Sam Chak]

Hi @LIULIYAN

If where and are the known control gains in Simulink, why is a variable?

If is Hurwitz and is positive-definite, then can be found, and you want to use it to compute the control gain matrix ? But this is usually found by solving the Algebraic Riccati Equation (ARE).

In other words, the Lyapunov equation is used to show stability, while the solution to the ARE is used to show optimality.

By the way, the lyap() function has not been supported for standalone code generation for many years. However, if you simply want to find via solving the Lyapunov equation, consider @Bruno Luong's two methods.

Answer 1

[Divyajyoti Nayak]

Hi LIULYAN,

You can use a MATLAB function block and call the lyap function inside it to calculate P.

Implement MATLAB Functions in Simulink with MATLAB Function Blocks - MATLAB & Simulink (mathworks.com)

[LIULIYAN]

I have attempted to use the 'MATLAB JUNCTION' module in Simulink to call the 'lyap' function, but an error will be reported during the simulation run. As shown in the picture below.

[Sam Chak]

Hi @Divyajyoti Nayak

Can you test if this simple Simulink model works on your machine?

P = lyap(-1, 1)

P = 0.5000

[LIULIYAN]

Hi@Sam Chak

When running the above Simulink model on my computer, it still reported an error. The version I am using is MATLAB 2016b.

Answer 2

[Bruno Luong]

I don't know simulink nor lyap function but just to tell you you can solve with standard algebra, so I hope it could be incorporated easier to simulink

A=[0,1;-1,-2]';
Q=[1,0;0,1];

P = (kron(eye(2),A') + kron(A.',eye(2))) \ -Q(:);
P = reshape(P,[2 2]),

P = 2×2

1.5000 -0.5000 -0.5000 0.5000

% Check
A'*P + P*A + Q

ans = 2×2

0 0 0 0

[Bruno Luong]

A variation to avoid using kron is solving by one of the linear solver where function-handle user supply is possible (instead of matrix representation)

A=[0,1;-1,-2]';
Q=[1,0;0,1];

P = lsqr(@(varargin) LyaProd(A, varargin{:}), -Q(:));

lsqr converged at iteration 3 to a solution with relative residual 2.5e-15.

P = reshape(P,[2 2]);
disp(P)

    1.5000   -0.5000
   -0.5000    0.5000

%%
function AP = LyaProd(A, P, opt)
P = reshape(P, [2 2]);
if strcmp(opt,'notransp')
    AP = A'*P + P*A;
else
    AP = A*P + P*A';
end
AP = AP(:);
end