Re: Lyapunov Equation
- To: mathgroup at christensen.cybernetics.net
- Subject: [mg713] Re: Lyapunov Equation
- From: Bernd.Cebulski at e-technik.tu-chemnitz.de (Bernd Cebulski)
- Date: Tue, 11 Apr 1995 10:29:50
- Organization: TU Chemnitz
In article <3md5uo$9fp at news0.cybernetics.net> ggb at bigdog.engr.arizona.edu (Gene Bergmeier) writes:
>From: ggb at bigdog.engr.arizona.edu (Gene Bergmeier)
>Subject: Lyapunov Equation
>Date: 11 Apr 1995 06:00:56 GMT
>Pardon the silly question, but my roommate and I are just beginning to
>use Mathematica. We are attempting to solve the Lyapunov equation (from
>nonlinear controls). The equation is
> PA+Trans(A)P=-Q where
> Q=I
> P is symmetric
> and A is general.
>We have tried the following procedure. Set Q equal to a 2D identity
>matrix, P as a matrix of p11, p12 and p22 and A as a matrix of known
>coefficients. Then we stated the Lyapunov eq and asked Mathematica to
>solve it for p11, p12 and p22.
>Any thoughts or suggestions?
..............................
Mathematica is in general not able to solve such equations directly.
For the special case of the Ljapunow equation there are at least two
ways to go. The first one is the calculation of the matrix exponentiation
e^A. For further hints please search in the specific literature. The second
way (which I prefer) is the use of the kronecker-product.
The basic relation here is: vec(F.G.H)= kronecker(Trans(H),F) . vec(G)
where 'vec' is the vectorize-operator (writes the matrices in columns, each
over the other).
So getting values of the matrix p would follow the equation:
vec(P) = - Inverse[ kronecker(Trans(A),In) + kronecker(Im,A)] . vec(Q)
(not standard mma-notation)
where 'In' and 'Im' is the n- or m-row identity-matrix.
Thats the (very rough) theory behind ...
Now my mma-solution for that:
----------------------- begin mma-code ----------------------------------
Needs["LinearAlgebra`MatrixManipulation`"];
rows[a_]:=Part[Dimensions[a],1]; (* get the rows of a matrix *)
columns[a_]:=Part[Dimensions[a],2];(* get the columns " *)
kron[a_,b_]:=BlockMatrix[Array[Part[Outer[Times,a,b],#1,#2]&,Dimensions[a]]];
tovector[a_]:=Flatten[Transpose[a]]; (* vectorize-operator)
tomatrix[a_]:=Transpose[Partition[a, Sqrt[rows[a]]]]; (* the inverse ...*)
ljapunow[a_,c_]:=tomatrix[ Inverse[kron[Transpose[a],
IdentityMatrix[rows[a]]] +
kron[IdentityMatrix[columns[a]],Transpose[a]]].
tovector[-Transpose[c]]]
(* gets the matrix P from the system matrix A and the symmetric positive
definite matrix C (perhaps identity-matrix, but not necessary) *)
------------------------ end mma-code ---------------------------------------
BTW, I used this for evaluation of stability for nonlinear systems with the
method of 'M. A. AISERMANN'. The matrix Q (in the above it is C) will be
optimized for a maximum range of (global) stability. For further
work please contact me per e-mail.
Cheerio, Bernd.
.-----------------------------------------------.
| Bernd Cebulski |
| Dep. of Electrical Machines and Drives |
| Chemnitz, University of Technology |
| Germany |
|--- bernd.cebulski at e-technik.tu-chemnitz.de ---|
| Phone: +49 (371) 5313318 |
| Fax: +49 (371) 5313361 |
()_________________\ /_________________()
________________ | | ________________
()____________| | | | | |____________()
______________\| | | |/______________
()_________________| |_________________()