Can MATLAB solve a Lyapunov equation symbolically?
Mostra commenti meno recenti
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 )
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.
9 Commenti
Sam Chak
il 28 Feb 2022
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 
?
so that the solution of X is in terms of ?
econogist
il 28 Feb 2022
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?
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
il 28 Feb 2022
Your wikipedia link directly says that the discrete case can be solved as an infinite sum.
econogist
il 28 Feb 2022
Walter Roberson
il 28 Feb 2022
"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
il 9 Mar 2022
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)
Notice the infinity in F(3,3)
assume( pe == 0 )
simplify(F)
assume(pe ~= 0)
sF = simplify(F)
c113 = children(sF(1,1),3)
c223 = children(sF(2,2),3)
assume(~c113);
assumeAlso(~c223);
simplify(sF)
But of course the infinity is still there
Sam Chak
il 10 Mar 2022
Hi @Matt Woolf
Can you check the state matrix J if it is possible for 
such that 
?
such that ?
econogist
il 10 Mar 2022
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.
. 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. Risposta accettata
Sam Chak
il 9 Mar 2022
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]
12 Commenti
econogist
il 9 Mar 2022
Walter Roberson
il 9 Mar 2022
Modificato: Walter Roberson
il 9 Mar 2022
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)
S1 = solve(eqns(1:5), [P(1,1),P(1,2),P(1,3),P(2,2),P(2,3)])
E2 = subs(eqns(6:end), S1)
symvar(E2)
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.
Walter Roberson
il 9 Mar 2022
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
il 9 Mar 2022
econogist
il 9 Mar 2022
Walter Roberson
il 10 Mar 2022
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
il 10 Mar 2022
Torsten
il 10 Mar 2022
So J is complex-valued, X is arbitrary complex (need not be Hermitian) and Q is supposed to be Hermitian ?
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)
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))
E2 = subs(eqns(N+1:end), S1).'
symvar(E2)
%}
%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
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 ?
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.
, and marginally stable if it has at least one eigenvalue on the unit circle.
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,
, 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
Più risposte (1)
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.
, , , 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.
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.
. 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.3 Commenti
econogist
il 10 Mar 2022
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
. 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.
produces the desired system eigenvalues that can be chosen with appropriate choice of the feedback gain matrix K through the pole placement technique.
econogist
il 10 Mar 2022
Categorie
Scopri di più su Matrix Computations in Centro assistenza e File Exchange
Prodotti
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
