Can MATLAB solve a Lyapunov equation symbolically?

I have a 3x3 matrix called J. I need to show that there exists at least one X and Q such that where Q is a Hermitian matrix and is the conjugate transpose of J (see https://en.wikipedia.org/wiki/Lyapunov_equation )

Matrix J is currently in terms of many parameters and variables which makes finding X and Q by hand difficult. Can MATLAB help with this? It seems using symbolic MATLAB would help but I'm not sure how to move forward.

[Sam Chak]

Hi @Matt Woolf

Looks like you are trying to compute the analytical solution for the discrete Lyapunov equation.

If it is, then it becomes a question of how to symbolically compute the infinite sum of products of matrices.

Can you assume so that the solution of X is in terms of ?

[econogist]

Yes, I'm trying to compute the analytical solution to the discrete Lyapunov equation. I suppose I could get the solution in terms of and then just insert the actual values I have for each . I don't see, however, why that means I want the solution to the infinite sum of products of matrices. Does the Lyapunov equation imply that? If so, why?

[Walter Roberson]

Your wikipedia link directly says that the discrete case can be solved as an infinite sum.

[econogist]

Ah, I see that now. So is the idea that MATLAB cannot symbolically solve the version of the Lyapunov equation I originally posted, but it can do so for the infinite sum of products? If so, how would I go about this?

[Walter Roberson]

"One may then solve for vec ( X ) by inverting or solving the linear equations."

which is something that you can do symbolically. The sum of infinite products representation is an alternative representation, not the only way to solve the problem.

[econogist]

I followed the advice above with the code below. I am confiused by its output though (also below). How should I interpret this? Does this mean the system is stable? Why is a piecewise function being returned?

syms k Q pe D B Ur Up Sr Sp
J=[(Q*pe-(D-B)*Ur)/(pe^2) -((D-B)*Up/(pe^2)) 0; -((D-B)*Ur/(pe^2)) (Q*pe-(D-B)*Up)/(pe^2) 0; Sr Sp 1];
Q=[1 0 0; 0 1 0; 0 0 1];
F=symsum((J.^k).*Q.*(ctranspose(J)),k,1,inf);
F


F =
 
[piecewise(pe ~= 0 & pe^2 + Ur*D == Q*pe + B*Ur, (Inf*(conj(Ur)*(conj(B) - conj(D)) + conj(Q)*conj(pe)))/conj(pe)^2, pe ~= 0 & pe^2 + Ur*D ~= Q*pe + B*Ur, ((conj(Ur)*(conj(B) - conj(D)) + conj(Q)*conj(pe))*(pe^2*limit(1/pe^(2*k)*(Q*pe + B*Ur - D*Ur)^k, k, Inf) - Q*pe - B*Ur + D*Ur))/(conj(pe)^2*(- pe^2 + Q*pe - Ur*D + B*Ur))),                                                                                                                                                                                                                                                                                                                                      0,   0]
[                                                                                                                                                                                                                                                                                                                                     0, piecewise(pe ~= 0 & pe^2 + Up*D == Q*pe + B*Up, (Inf*(conj(Up)*(conj(B) - conj(D)) + conj(Q)*conj(pe)))/conj(pe)^2, pe ~= 0 & pe^2 + Up*D ~= Q*pe + B*Up, ((conj(Up)*(conj(B) - conj(D)) + conj(Q)*conj(pe))*(pe^2*limit(1/pe^(2*k)*(Q*pe + B*Up - D*Up)^k, k, Inf) - Q*pe - B*Up + D*Up))/(conj(pe)^2*(- pe^2 + Q*pe - Up*D + B*Up))),   0]
[         

[Walter Roberson]

Are you able to put any constraints on your variables? Real-valued? Integer? positive?

The piecewise pe^2 + Ur*D == Q*pe + B*Ur part is saying that it found (pe^2 + Ur*D) - (Q*pe + B*Ur) or algebraic equivalent in a denominator, and when that condition holds, there is a division by 0 that leads to an infinite outcome.

syms k Q pe D B Ur Up Sr Sp
J=[(Q*pe-(D-B)*Ur)/(pe^2) -((D-B)*Up/(pe^2)) 0; -((D-B)*Ur/(pe^2)) (Q*pe-(D-B)*Up)/(pe^2) 0; Sr Sp 1];
Q=[1 0 0; 0 1 0; 0 0 1];
F=symsum((J.^k).*Q.*(ctranspose(J)),k,1,inf)

F = 

Notice the infinity in F(3,3)

assume( pe == 0 )
simplify(F)

ans = 

assume(pe ~= 0)
sF = simplify(F)

sF = 

c113 = children(sF(1,1),3)

c113 = 

c223 = children(sF(2,2),3)

c223 = 

assume(~c113);
assumeAlso(~c223);
simplify(sF)

ans = 

But of course the infinity is still there

[Sam Chak]

Hi @Matt Woolf

Can you check the state matrix J if it is possible for such that ?

[econogist]

It's not possible that . That's why I'm having to go the Lyapunov route in the first place really since when I tried to prove stability via eigenvalues alone (hoping all three would be less than 0), the dominant eigenvalue was equal to 1.

Answer 1

[Sam Chak]

Hi @Matt Woolf

Try this script:

syms k Q pe D B Ur Up Sr Sp p11 p12 p13 p22 p23 p33

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

A = [(Q*pe-(D-B)*Ur)/(pe^2) -((D-B)*Up/(pe^2)) sym('0');
     -((D-B)*Ur/(pe^2)) (Q*pe-(D-B)*Up)/(pe^2) sym('0');
     Sr Sp sym('1')];

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*A.' - P + Q;	% discrete-time Lyapunov equations

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]

Credits to @Torsten.

See also: https://www.mathworks.com/matlabcentral/answers/1666739-use-symbolic-variable-for-lyapunov-function

[econogist]

Thanks, this looks amazing! Does this code solve when you run it though? The output tells me "Unrecognized field name "p33"" and that Maltab is unable to find an explicit solution.

[Walter Roberson]

syms k Q pe D B Ur Up Sr Sp p11 p12 p13 p22 p23 p33

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

A = [(Q*pe-(D-B)*Ur)/(pe^2) -((D-B)*Up/(pe^2)) sym('0');
     -((D-B)*Ur/(pe^2)) (Q*pe-(D-B)*Up)/(pe^2) sym('0');
     Sr Sp sym('1')];

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*A.' - P + Q;	% discrete-time Lyapunov equations

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];
symvar(eqns)

ans = 

S1 = solve(eqns(1:5), [P(1,1),P(1,2),P(1,3),P(2,2),P(2,3)])

S1 = struct with fields:
    p11: -(pe^2*(2*B^3*Q*Up^3 + 2*Ur*B^3*Q*Up^2 + 4*B^2*Q^2*Up^2*pe + 2*Ur*B^2*Q^2*Up*pe - 6*D*B^2*Q*Up^3 - 6*Ur*D*B^2*Q*Up^2 - 2*Ur*B^2*Up*pe^3 + 3*B*Q^3*Up*pe^2 + Ur*B*Q^3*pe^2 - 8*D*B*Q^2*Up^2*pe - 4*Ur*D*B*Q^2*Up*pe + 6*D^2*B*Q*Up^3 + 6*Ur*D^2*B*Q…
    p12: (pe^2*(B - D)*(2*B^2*Q*Up^2*Ur + 2*B^2*Q*Up*Ur^2 + B*Q^2*Up^2*pe + 4*B*Q^2*Up*Ur*pe + B*Q^2*Ur^2*pe - 4*D*B*Q*Up^2*Ur - 4*D*B*Q*Up*Ur^2 - B*Up^2*pe^3 - B*Ur^2*pe^3 + Q^3*Up*pe^2 + Q^3*Ur*pe^2 - D*Q^2*Up^2*pe - 4*D*Q^2*Up*Ur*pe - D*Q^2*Ur^2*pe…
    p13: (pe^2*(Q^2*Sr*pe^8 + 2*Q^3*Sr*pe^7 - 2*Q^4*Sr*pe^6 - Q^5*Sr*pe^5 + Q^6*Sr*pe^4 - Q*Sr*pe^9 - B*Sp*Up*pe^8 - B*Sr*Ur*pe^8 + D*Sp*Up*pe^8 + D*Sr*Ur*pe^8 + 2*B^4*Q^2*Sr*Up^4 + 2*D^4*Q^2*Sr*Up^4 - B^3*Sp*Ur^3*pe^4 - B^3*Sr*Up^3*pe^4 + D^3*Sp*Ur^3…
    p22: -(pe^2*(2*B^3*Q*Ur^3 + 2*Up*B^3*Q*Ur^2 + 4*B^2*Q^2*Ur^2*pe + 2*Up*B^2*Q^2*Ur*pe - 6*D*B^2*Q*Ur^3 - 6*Up*D*B^2*Q*Ur^2 - 2*Up*B^2*Ur*pe^3 + 3*B*Q^3*Ur*pe^2 + Up*B*Q^3*pe^2 - 8*D*B*Q^2*Ur^2*pe - 4*Up*D*B*Q^2*Ur*pe + 6*D^2*B*Q*Ur^3 + 6*Up*D^2*B*Q…
    p23: (pe^2*(Q^2*Sp*pe^8 + 2*Q^3*Sp*pe^7 - 2*Q^4*Sp*pe^6 - Q^5*Sp*pe^5 + Q^6*Sp*pe^4 - Q*Sp*pe^9 - B*Sp*Up*pe^8 - B*Sr*Ur*pe^8 + D*Sp*Up*pe^8 + D*Sr*Ur*pe^8 + 2*B^4*Q^2*Sp*Ur^4 + 2*D^4*Q^2*Sp*Ur^4 - B^3*Sp*Ur^3*pe^4 - B^3*Sr*Up^3*pe^4 + D^3*Sp*Ur^3…

E2 = subs(eqns(6:end), S1)

E2 = 

symvar(E2)

ans = 

Notice that p33 is not in the variables for the last remaining equation. Indeed, p33 is not in eqns at all.

But also notice that Q is a variable . And if you look at your equations you use Q and then later set it to a value, but you do not substitute it into the existing uses. If you do substitute it into existing uses, then even more variables vanish.

And remember that you first use Q as-if it is a scalar, but later you drop a matrix in on it. That invalidates assumptions: scalar Q * somthing has a different meaning that matrix Q * something.

[econogist]

How can that be the answer though? The matrix dimensions clearly won't work to solve the Lyanpunov equation since you'd have to multiply a 3x3 matrix by a 1x8 matrix which is impossible.

[Walter Roberson]

Not sure where you find a 1 x 8 in there ?

But it is true that if Q starts out 3 x 3 then you cannot build A correctly -- and if you duplicate elements to match dimensions, the result would certainly not be 3 x 3.

That suggests to me that the Q that appears in A is not the same as the Q that is assigned eye(3) . Or perhaps the Q that appears in A needs to be indexed.

[econogist]

The 1 x 8 I'm seeing is ans=(B D Q Sp Sr Up Ur pe). Is this supposed to be value of P that will solve the Lyapunov equation?

I see what you mean with Q. You're correct that the Q that appears in A is not the same as the Q assigned eye(3). Sorry for that confusion. Perhaps the Q that appears in A could be called R.

[econogist]

I should also mention that Q does not have to be equal to eye(3) and just has to be any Hermitian matrix...

[Walter Roberson]

No, you have a system of 6 equations in 13 variables. E2 is the 6th equation after solutions for 5 of the variables have been found, so E2 is the equation that in your code you would ideally have solved for p33.

However, when you look at the list of variables involved in that last equation, E2, you see that list of 8 variables, and you see that the last equation does not involve p33 at all, so you cannot solve for it. And if you track backwards, you will find that the whole set of equations in eqns does not involve p33 -- so there was never any chance to solve for p33 (even if you had asked to solve for it as the very first thing.)

That list of 8 values is not the solution to the equation: it is the list of variables involved in the final equation after we have solved the first 5 equations.

Is the Q that is involved in creating A the same as the Q that is added in creating L ? If so then how is the size difference to be solved?

[econogist]

No, the Q that is involved in creating A is not the same as that added in creating L. I should have used a different letter there.

So are you saying a solution does not exist?

[Torsten]

So J is complex-valued, X is arbitrary complex (need not be Hermitian) and Q is supposed to be Hermitian ?

[Walter Roberson]

syms k q pe D B Ur Up Sr Sp p11 p12 p13 p22 p23 p33

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

A = [(q*pe-(D-B)*Ur)/(pe^2) -((D-B)*Up/(pe^2)) sym('0');
     -((D-B)*Ur/(pe^2)) (q*pe-(D-B)*Up)/(pe^2) sym('0');
     Sr Sp sym('1')];

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

N = sym(zeros(3));	% zero matrix
L = A*P*A.' - P + Q;	% discrete-time Lyapunov equations

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];
symvar(eqns)

ans = 

vars = [P(1,1),P(1,2),P(1,3),P(2,2),P(2,3)];
N = 4;
S1 = solve(eqns(1:N), vars(1:N))

S1 = struct with fields:
    p11: -(pe^2*(Q1_1*pe^7 + Q1_1*pe^3*q^4 - 2*Q1_1*pe^5*q^2 + B^3*Q1_1*Up^3*q - 2*B^3*Q1_2*Up^3*q + B^3*Q2_2*Up^3*q - D^3*Q1_1*Up^3*q + 2*D^3*Q1_2*Up^3*q - D^3*Q2_2*Up^3*q - B^2*Q1_1*Up^2*pe^3 + B^2*Q2_2*Up^2*pe^3 - D^2*Q1_1*Up^2*pe^3 + D^2*Q2_2*Up^2…
    p12: -(pe^2*(Q1_2*pe^7 + Q1_2*pe^3*q^4 - 2*Q1_2*pe^5*q^2 - B^2*Q1_2*Up^2*pe^3 + B^2*Q2_2*Up^2*pe^3 + B^2*Q1_1*Ur^2*pe^3 - B^2*Q1_2*Ur^2*pe^3 - D^2*Q1_2*Up^2*pe^3 + D^2*Q2_2*Up^2*pe^3 + D^2*Q1_1*Ur^2*pe^3 - D^2*Q1_2*Ur^2*pe^3 + B^2*Q1_2*Up^2*pe*q^2…
    p13: -(- Q1_3*pe^11 + Q1_3*pe^5*q^6 - 3*Q1_3*pe^7*q^4 + 3*Q1_3*pe^9*q^2 - B*Up*p23*pe^9 + D*Up*p23*pe^9 - Q1_2*Sp*pe^10*q - Q1_1*Sr*pe^10*q - Q1_2*Sp*pe^6*q^5 + 2*Q1_2*Sp*pe^8*q^3 - Q1_1*Sr*pe^6*q^5 + 2*Q1_1*Sr*pe^8*q^3 + B^2*Q1_3*Up^2*pe^7 + B^2*…
    p22: -(pe^2*(Q2_2*pe^7 + Q2_2*pe^3*q^4 - 2*Q2_2*pe^5*q^2 + B^3*Q1_1*Ur^3*q - 2*B^3*Q1_2*Ur^3*q + B^3*Q2_2*Ur^3*q - D^3*Q1_1*Ur^3*q + 2*D^3*Q1_2*Ur^3*q - D^3*Q2_2*Ur^3*q + B^2*Q1_1*Ur^2*pe^3 - B^2*Q2_2*Ur^2*pe^3 + D^2*Q1_1*Ur^2*pe^3 - D^2*Q2_2*Ur^2…

E2 = subs(eqns(N+1:end), S1).'

E2 = 

symvar(E2)

ans = 

%}
%S2 = solve(E2, P(3,3))
%results = [S1.p11; S1.p12; S1.p13; S1.p22; S1.p23; P(3,3)];
%results = subs(results, P(3,3), S2)

You can solve E2(1) for p23 and substitute that in, but you can see that the final result would not have p33 in it, so no matter what Q matrix you use, you canot solve for p33

[Torsten]

@Matt Woolf

You write:

I have a 3x3 matrix called J. I need to show that there exists at least one X and Q such that

J*X*J^H - X + Q = 0 where Q is a Hermitian matrix and J^H is the conjugate transpose of J.

I don't understand why you don't take an arbitrary Hermitian matrix X (e.g. X = 0) and set Q = X - J*X*J^H which is Hermitian, so fulfills the requirement. Or where do I go wrong here ?

[Sam Chak]

Thanks @Walter Roberson for showing how to solve this rigorously.

The stability conditions for asymptotic stability says that a discrete-time linear system is asymptotically stable if its eigenvalues are inside the unit circle , and marginally stable if it has at least one eigenvalue on the unit circle.

I have double-checked the state matrix supplied by @Matt Woolf, that

will always produce at least one eigenvalue on the unit circle.

Naturally, since , when trying to solve the discrete-time Lyapunov equation, , then will be eliminated in the process. For example,

A = [-0.1 -0.2 0; 0.3 0.4 0; 0.5 0.6 1]
eig(A)
Q = eye(3);
P = dlyap(A, Q)

the chosen state matrix will produce an error message:

Here is a second example that a state matrix with all its eigenvalues are inside the unit circle.

clear all
A = [-0.2 -0.2 0.4; 0.5 0 1; 0 -0.4 -0.5]
eig(A)          % check the eigenvalues of A are inside the unit circle
Q = eye(3)      % identity matrix
P = dlyap(A, Q)
eig(P)          % check the signs of the eigenvalues of P
A*P*A' - P + Q  % check if it produces the zero matrix

Because the eigenvalues of P are all positive, the matrix is positive definite and the discrete-time system is asymptotically stable.

By the way, Q does not have be an identity matrix. It is possible to investigate the stability of the system using Q = diag([0, 0, 1]), so long as the rank of the observability matrix is full rank. 3rd example:

clear all
A = [-0.2 -0.2 0.4; 0.5 0 1; 0 -0.4 -0.5]
eig(A)
C = [0 0 1];        % output matrix
Q = C'*C            % a Hermitian matrix
rank(obsv(A, C))    % check the rank of the observability matrix
P = dlyap(A, Q)
eig(P)
A*P*A' - P + Q

Answer 2

[Sam Chak]

Hi @Matt Woolf

Let's try a new one by reducing the non-dependant variables. I have simplified the computation by letting , , , and , with and maintained. Also, I have replaced with in order to investigate whether the solution exists.

clear all; clc
syms a b c d Sr Sp p11 p12 p13 p22 p23 p33

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

A = [a b sym('0');
     c d sym('0');
     Sr Sp sym('0.5')];

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

Q = sym(eye(3));	% identity matrix

L = A*P*A.' - P + Q;	% discrete-time Lyapunov equations

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]
simplify(sol)

pretty(S.p33)

In this example, the solution exists, though is still pretty long even after simplification, not mentioning that you have to substitute k, q, pe, D, B, Ur, Up, back into the a, b, c, d.

More importantly, you need to recheck your derivation/modeling of the state matrix if . If it is, it could imply that the system is marginally stable and you can't find any P that satisfies the discrete Lyapunov equation.

[econogist]

The reason I'm tying to solve the Lyapunov equation in the first place is due to this pesky eigenvalue equal to 1. There's no way to change that value I'm afraid.

I'm fairly certain the system is stable but just can't prove it analytically. I was hoping I'd be able to get a fairly simple matrix for P in terms of the variables I've written that would solve the Lyapunov equation for certain parameter/variable values.

[Sam Chak]

Hi @Matt Woolf

The stability issue maybe not as bad as you think. Although I don't know what your system really is, if your system

is controllable, then there exists a feedback gain matrix K in the controller

that arbitrarily assigns the system eigenvalues to any set . When the two equations are combined, you can rewrite the closed-loop system as

,

where the closed-loop state matrix is

.

In other words, produces the desired system eigenvalues that can be chosen with appropriate choice of the feedback gain matrix K through the pole placement technique.

[econogist]

Thanks so much for this response, Sam. I'm going to have to sit with it for a bit as some of what you wrote is beyond my knowledge.

The full system is here, by the way: https://math.stackexchange.com/questions/4398126/is-this-3-equation-nonlinear-system-stable