Difference between revisions of "Inverse Lyapunov Procedure"

Revision as of 14:07, 10 January 2018

Description

The Inverse Lyapunov Procedure (ilp) is a synthetic random linear system generator. It is based on reversing the Balanced Truncation procedure and was developed in ^[1], where a description of the algorithm is given. In aggregate form, for randomly generated controllability and observability gramians, a balancing transformation is computed. The balanced gramian is the basis for an associated state-space system, which is determined by solving a Lyapunov equation and then unbalanced. A central point is the solution of the Lyapunov equations for the system matrix instead of the gramian matrix. This is feasable due to the symmetric (semi-)positive definiteness of the gramians and the requirement for a stable system, yet with a non-unique solution.

Implementation

An efficient approach to solving the Lyapunov equation is provided by empirical gramians.

Usage

Use the following matlab code to generate a random system as described above:

function [A,B,C] = ilp(M,N,Q,s,r)
% ilp (inverse lyapunov procedure)
% by Christian Himpe, 2013--2018
% released under BSD 2-Clause License
%*
  if(nargin==5)
    rand('seed',r);
    randn('seed',r);
  end;

% Gramian Eigenvalues
  WC = exp(0.5*rand(N,1));
  WO = exp(0.5*rand(N,1));

% Gramian Eigenvectors
  [P,S,R] = svd(randn(N));

% Balancing Transformation
  WC = P*diag(sqrt(WC))*P';
  WO = R*diag(sqrt(WO))*R';
  [U,D,V] = svd(WC*WO);

% Input and Output
  B = randn(N,M);

  if(nargin>=4 && s~=0)
    C = B';
  else
    C = randn(Q,N);
  end

% Scale Output Matrix
  BB = sum(B.*B,2);  % = diag(B*B')
  CC = sum(C.*C,1)'; % = diag(C'*C)
  C = bsxfun(@times,C,sqrt(BB./CC)');

% Solve System Matrix
  A = -sylvester(D,D,B*B');

% Unbalance System
  A = V*A*U';
  B = V*B;
  C = C*U';

end

The function call requires three parameters; the number of inputs $M$ , of states $N$ and outputs $Q$ . Optionally, a symmetric system can be enforced with the parameter $s \neq 0$ . For reproducibility, the random number generator seed can be controlled by the parameter $r \in \mathbb{N}$ . The return value consists of three matrices; the system matrix $A$ , the input matrix $B$ and the output matrix $C$ .

[A,B,C] = ilp(M,N,Q,s,r);

A variant of the above code using empirical Gramians instead of a matrix equation solution can be found at http://gramian.de/utils/ilp.m , which may yield preferable results.

References

↑ S.C. Smith, J. Fisher, "On generating random systems: a gramian approach", Proceedings of the American Control Conference, 2003.

Contact

Christian Himpe

[smith03-1] S.C. Smith, J. Fisher, "On generating random systems: a gramian approach", Proceedings of the American Control Conference, 2003.

[1]

@@ Line 24: / Line 24: @@
 <source lang="matlab">
-function [A B C] = ilp(J,N,O,s,r)
+function [A,B,C] = ilp(M,N,Q,s,r)
 % ilp (inverse lyapunov procedure)
 % by Christian Himpe, 2013--2018
@@ Line 35: / Line 35: @@
 % Gramian Eigenvalues
-  WC = exp(rand(N,1));
+  WC = exp(0.5*rand(N,1));
-  WO = exp(rand(N,1));
+  WO = exp(0.5*rand(N,1));
 % Gramian Eigenvectors
-  [P,S,Q] = svd(randn(N));
+  [P,S,R] = svd(randn(N));
 % Balancing Transformation
-  WC = P*diag(WC)*P';
+  WC = P*diag(sqrt(WC))*P';
-  WO = Q*diag(WO)*Q';
+  WO = R*diag(sqrt(WO))*R';
   [U,D,V] = svd(WC*WO);
 % Input and Output
-  B = randn(N,J);
+  B = randn(N,M);
   if(nargin>=4 && s~=0)
     C = B';
   else
-    C = randn(O,N);
+    C = randn(Q,N);
   end
@@ Line 72: / Line 72: @@
 <!--]]--></div>
-The function call requires three parameters; the number of inputs <math>J</math>, of states <math>N</math> and outputs <math>O</math>.
+The function call requires three parameters; the number of inputs <math>M</math>, of states <math>N</math> and outputs <math>Q</math>.
 Optionally, a symmetric system can be enforced with the parameter <math>s \neq 0</math>.
 For reproducibility, the random number generator seed can be controlled by the parameter <math>r \in \mathbb{N}</math>.
@@ Line 78: / Line 78: @@
 :<source lang="matlab">
-[A,B,C] = ilp(J,N,O,s,r);
+[A,B,C] = ilp(M,N,Q,s,r);
 </source>
+A variant of the above code using empirical Gramians instead of a matrix equation solution can be found at http://gramian.de/utils/ilp.m , which may yield preferable results.
 ==References==