Use symbolic variable for lyapunov function

I am trying to find a value for a lyapunov function but I do not know the numeric values. When I run the lyapunov command, I get an error that only numeric arrays can be used. Is there a way for using only symbolic variable to get the answer.

[Sam Chak]

Hi @Kashish Pilyal

Without showing your work and the relevant info about your system, it's very difficult to provide assistance.

[Kashish Pilyal]

 syms k_p k_d h integer
 
 A=[0     1    0;
   -k_p -k_d   0;
    0    1/h  -1/h];

Q=eye(3);

P=lyap(A',Q);

So, I need to find the value of matrix P in terms of symbols in matrix A. The MATLAB gives an error that lyap function will only take numeric values. Is there any other way as by hand it is getting too messy.

[Sam Chak]

Hi @Kashish Pilyal

I have tested and verified the results symbolically that holds.

clear all; clc
syms a b c

A = sym('A', [3 3]);	% state matrix
P = sym('P', [3 3]);	% positive definite matrix

A = [sym('0') sym('1') sym('0');
     -a -b sym('0');
     sym('0') c -c];

P = [((a^3 + 2*a^2*b*c + 2*a^2*c^2 + a^2 + a*b^2 + a*b*c + a*c^2 + b^3*c + b^2*c^2)/(2*a*b*(c^2 + b*c + a))) (1/(2*a)) (-a/(2*(c^2 + b*c + a)));
     (1/(2*a)) ((a^2 + 2*a*c^2 + b*a*c + a + c^2 + b*c)/(2*a*b*(c^2 + b*c + a))) (c/(2*(c^2 + b*c + a)));
     (-a/(2*(c^2 + b*c + a))) (c/(2*(c^2 + b*c + a))) (1/(2*c))];

Q = sym(eye(3));	% identity matrix
L = A.'*P + P*A + Q;	% Lyapunov equation
simplify(L)

Result:

[Kashish Pilyal]

@Sam Chak Thank you for the help.

Answer 1

[Sam Chak]

Hi @Kashish Pilyal

If you are writing for a journal paper or a thesis, the following explanation might be helpful.

Let , , and .

There are a few ways to solve this symbolically.

syms a b c p11 p12 p22 p23 p33 p31
eqns = [1 - 2*a*p12 == 0, - a*p22 - b*p12 + c*p31 + p11 == 0, 1 - 2*b*p22 + 2*c*p23 + 2*p12 == 0, - b*p23 - c*p23 + c*p33 + p31 == 0, 1 - 2*c*p33 == 0, - a*p23 - c*p31 == 0];
S = solve(eqns);
sol = [S.p11; S.p12; S.p22; S.p23; S.p33; S.p31]

Result:

The result has been verified numerically:

clear all; clc
A = [0 1 0; -1 -2 0; 0 1 -1]
Q = eye(3)
P = lyap(A', Q)
A'*P + P*A

[Kashish Pilyal]

Thank you for the answer as it will help me in my report, but I need to ask one thing. To verify the answer, I have to show that numerically? As when I use only symbolic variables, the equation does not hold true and I think that could be a problem in my report as the results do not match. Is there a way out of this?

[Sam Chak]

Improved Script (merging with @Torsten's code):

syms k_p k_d h p11 p12 p13 p22 p23 p33

A = sym('A', [3 3]);	% state matrix
P = sym('P', [3 3]);	% positive definite matrix

A = [sym('0') sym('1') sym('0');
     -k_p -k_d sym('0');
     sym('0') sym('1')/h sym('-1')/h];

P = [p11 p12 p13;
     p12 p22 p23;
     p13 p23 p33];

Q = sym(eye(3));	% identity matrix
N = sym(zeros(3));	% zero matrix
L = A.'*P + P*A + Q;	% Lyapunov equation

eqns = [L(1, 1) == 0, L(1, 2) == 0, L(1, 3) == 0, L(2, 2) == 0, L(2, 3) == 0, L(3, 3) == 0];
S = solve(eqns);
sol = [S.p11; S.p12; S.p13; S.p22; S.p23; S.p33]

[Torsten]

I don't understand why you don't have to use

p21,p31,p32

as symbolic and solve

L(1,1)==0, L(1,2)==0, L(1,3)==0, L(2,1)==0, L(2,2)==0, L(2,3)==0, L(3,1)==0, L(3,2)==0, L(3,3)==0

.

[Sam Chak]

Hi @Torsten

My apologies for failing to inform that P has to be a symmetric matrix . Allow me to quote the theorem directly from Prof. Hassan Khalil's book, "Nonlinear Control":

Theorem: A matrix A is Hurwitz if and only if for every positive definite symmetric matrix Q, there exists a positive definite symmetric matrix P that satisfies the Lyapunov equation . Moreover, if A is Hurwitz, then P is the unique solution.

From the property of symmetry, we know that , , and .

I'm still learning and not good at expressing the control law and equations in the symbolic form in MATLAB. That's why I worked out the equations manually and then used MATLAB to solve the derived set of linear equations. Thanks for your original script in solving the symbolic equations.

[Torsten]

Ah, I didn't know this.

Thank you for the info.

Answer 2

[Torsten]

syms k_p k_d h
A = sym('A', [3 3]);
X = sym('X', [3 3]);
A = [sym('0') sym('1') sym('0');
     -k_p -k_d sym('0');
     sym('0') sym('1')/h  sym('-1')/h];
Q = sym(eye(3));
N = sym(zeros(3));
B = A.'*X + X*A + Q;
F = solve(B==N)

[Kashish Pilyal]

Thank you for the answer but I have tried this method too. The matrix F in this case comes out to be empty. It is 0 by 1 symbolic. I actually managed to get the answer now. I had to write all equations seperately like this

syms P [3 3]
X= (A'*P)+(P*A);
eqnA=X(3,3)==-1;
eqnB=X(3,2)==0;
eqnC=X(3,1)==0;
eqnD=X(2,3)==0;
eqnE=X(1,3)==0;
eqnF=X(1,1)==-1;
eqnG=X(2,1)==0;
eqnH=X(2,2)==-1;
eqnI=X(1,2)==0;

Then I used solve and got the values of each element although the lyapunov function equation A'P+PA=-Q (in my case I) does not hold properly still. In other words the lyapunov equation does not give the identity matrix as output.