Inverse Lyapunov Procedure - Revision history

Schuetz2 at 05:32, 17 June 2025

2025-06-17T05:32:31Z

← Older revision		Revision as of 05:32, 17 June 2025
Line 131:		Line 131:
	title = {Inverse Lyapunov Procedure},		title = {Inverse Lyapunov Procedure},
	howpublished = {{MORwiki} -- Model Order Reduction Wiki},		howpublished = {{MORwiki} -- Model Order Reduction Wiki},
	url = <nowiki>{~~http~~://modelreduction.org/~~index.php~~/Inverse_Lyapunov_Procedure}</nowiki>,		url = <nowiki>{https://modelreduction.org/morwiki/Inverse_Lyapunov_Procedure}</nowiki>,
	year = 2018		year = 2018
	}		}

Mitchell: typo

2022-03-15T10:13:42Z

typo

← Older revision		Revision as of 10:13, 15 March 2022
Line 20:		Line 20:
	# Generate random orthogonal matrices (ie.: SVD of random matrix).		# Generate random orthogonal matrices (ie.: SVD of random matrix).
	# Compute balancing transformation for these Gramians.		# Compute balancing transformation for these Gramians.
	# Sample random input and ~~out~~ matrices.		# Sample random input and output matrices.
	# Scale output matrix to input matrix.		# Scale output matrix to input matrix.
	# Solve Lyapunov equation for system matrix.		# Solve Lyapunov equation for system matrix.

Himpe: added procedural cat

2018-08-22T12:27:00Z

added procedural cat

← Older revision		Revision as of 12:27, 22 August 2018
Line 1:		Line 1:
	[[Category:benchmark]]		[[Category:benchmark]]
			[[Category:procedural]]
	[[Category:linear]]		[[Category:linear]]
	[[Category:first differential order]]		[[Category:first differential order]]

Werner: Updated BibTeX entry.

2018-04-24T09:01:19Z

Updated BibTeX entry.

← Older revision		Revision as of 09:01, 24 April 2018
Line 124:		Line 124:

	* For the benchmark itself and its data:		* For the benchmark itself and its data:
	::The MORwiki Community. '''Inverse Lyapunov Procedure'''. MORwiki - Model Order Reduction Wiki, 2018. http://modelreduction.org/index.php/Inverse_Lyapunov_Procedure		::The MORwiki Community, '''Inverse Lyapunov Procedure'''. MORwiki - Model Order Reduction Wiki, 2018. http://modelreduction.org/index.php/Inverse_Lyapunov_Procedure

	@MISC{morwiki-invlyapproc,		@MISC{morwiki-invlyapproc,
	author = {The {MORwiki} ~~Community~~},		author = <nowiki>{{The MORwiki Community}}</nowiki>,
	title = {Inverse Lyapunov Procedure},		title = {Inverse Lyapunov Procedure},
	howpublished = {{MORwiki} -- Model Order Reduction Wiki},		howpublished = {{MORwiki} -- Model Order Reduction Wiki},
	url = {<nowiki>http://modelreduction.org/index.php/Inverse_Lyapunov_Procedure</nowiki>},		url = <nowiki>{http://modelreduction.org/index.php/Inverse_Lyapunov_Procedure}</nowiki>,
	year = {2018}		year = 2018
	}		}

	* For the background on the benchmark:		* For the background on the benchmark:

	@INPROCEEDINGS{SmiF03,		@INPROCEEDINGS{SmiF03,
	author = {Smith, S.~C. and Fisher, J.},		author = {Smith, S.~C. and Fisher, J.},
	title = {On generating random systems: a gramian approach},		title = {On generating random systems: a gramian approach},
	booktitle = {Proc. Am. Control. Conf.},		booktitle = {Proc. Am. Control. Conf.},
	volume = 3,		volume = 3,
	pages = {2743--2748},		pages = {2743--2748},
	year = 2003,		year = 2003,
	doi = {10.1109/ACC.2003.1243494}		doi = {10.1109/ACC.2003.1243494}
	}		}

	==References==		==References==

Himpe: more description and some fixes.

2018-03-06T15:12:23Z

more description and some fixes.

← Older revision		Revision as of 15:12, 6 March 2018
Line 6:		Line 6:

	==Description==		==Description==
	The '''Inverse Lyapunov Procedure''' (~~ilp~~) is a synthetic random linear system generator.		The '''Inverse Lyapunov Procedure''' (ILP) is a synthetic random linear system generator.
	It is based on reversing the [[Balanced_Truncation\|Balanced Truncation]] procedure and was developed in <ref name="smith03"/>,		It is based on reversing the [[Balanced_Truncation\|Balanced Truncation]] procedure and was developed in <ref name="smith03"/>,
	where a description of the algorithm is given.		where a description of the algorithm is given.
Line 14:		Line 14:
	A central point is the solution of the Lyapunov equations for the system matrix instead of the gramian matrix.		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.		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.
			Following, the steps of the ILP are listed:

			# Sample eigenvalues of controllability and observability Gramians.
			# Generate random orthogonal matrices (ie.: SVD of random matrix).
			# Compute balancing transformation for these Gramians.
			# Sample random input and out matrices.
			# Scale output matrix to input matrix.
			# Solve Lyapunov equation for system matrix.
			# Unbalance system.

	===Inverse Sylvester Procedure===		===Inverse Sylvester Procedure===
	A variant of the '''Inverse Lyapunov Procedure''' is the inverse Sylvester procedure <ref name="himpe17"/>.		A variant of the '''Inverse Lyapunov Procedure''' is the inverse Sylvester procedure (ILS) <ref name="himpe17"/>,
			which generates only state-space symmetric systems.
			Instead of balanced truncation, the [[wikipedia:Cross_Gramian\|cross Gramian]] is utilized for the random system generation, and hence a [[wikipedia:Sylvester_equation\|Sylvester equation]] is needs to be solved.
			The steps for the ILS are listed below:

			# Sample cross Gramian eigenvalues.
			# Sample random input matrix, and set output matrix as its transpose.
			# Solve Sylvester equation for system matrix.
			# Sample orthogonal unbalancing transformation (QR of random matrix).
			# Unbalance system.

			Even though the ILS is more limited than the ILP, for large systems it can be more efficient.

	==Data==		==Data==
Line 99:		Line 119:
	<math>B \in \mathbb{R}^{N \times M}</math>,		<math>B \in \mathbb{R}^{N \times M}</math>,
	<math>C \in \mathbb{R}^{Q \times N}</math>.		<math>C \in \mathbb{R}^{Q \times N}</math>.


	==Citation==		==Citation==


	To cite this benchmark, use the following references:		To cite this benchmark, use the following references:

Line 113:		Line 130:
	title = {Inverse Lyapunov Procedure},		title = {Inverse Lyapunov Procedure},
	howpublished = {{MORwiki} -- Model Order Reduction Wiki},		howpublished = {{MORwiki} -- Model Order Reduction Wiki},
	url = {<nowiki>http://modelreduction.org/index.php/~~Inverse Lyapunov Procedure~~</nowiki>},		url = {<nowiki>http://modelreduction.org/index.php/Inverse_Lyapunov_Procedure</nowiki>},
	year = {2018}		year = {2018}
	}		}
Line 128:		Line 145:
	doi = {10.1109/ACC.2003.1243494}		doi = {10.1109/ACC.2003.1243494}
	}		}


	==References==		==References==

Himpe: Added missing section and various fixes

2018-03-06T14:53:32Z

Added missing section and various fixes

@@ Line 7: / Line 7: @@
 ==Description==
 The '''Inverse Lyapunov Procedure''' (ilp) is a synthetic random linear system generator.
-It is based on reversing the [[Balanced_Truncation|Balanced Truncation]] procedure and was developed in <ref name="smith03">S.C. Smith, J. Fisher, "<span class="plainlinks">[http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1243494 On generating random systems: a gramian approach]</span>", Proceedings of the American Control Conference, 2003.</ref>, where a description of the algorithm is given.
+It is based on reversing the [[Balanced_Truncation|Balanced Truncation]] procedure and was developed in <ref name="smith03"/>,
 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 [[wikipedia:Lyapunov_equation|Lyapunov equation]] and then unbalanced.
+The balanced gramian is the basis for an associated state-space system,
 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==
+===Inverse Sylvester Procedure===
-An efficient approach to solving the Lyapunov equation is provided by empirical gramians.
+==Data==
-==Usage==
+<div class="thumbinner" style="width:20%;text-align:left;"><!--[[Media:ilp.m|-->
 <source lang="matlab">
@@ Line 82: / Line 83: @@
 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==
-<references />
+<references>
 ==Contact==
 [[User:Himpe|Christian Himpe]]

Himpe: /* Usage */ another code update

2018-01-10T12:07:05Z

Usage: another code update

← Older revision		Revision as of 12:07, 10 January 2018
Line 24:		Line 24:
	<source lang="matlab">		<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)		% ilp (inverse lyapunov procedure)
	% by Christian Himpe, 2013--2018		% by Christian Himpe, 2013--2018
Line 35:		Line 35:

	% Gramian Eigenvalues		% 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		% Gramian Eigenvectors
	[P,S,Q] = svd(randn(N));		[P,S,R] = svd(randn(N));

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

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

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

Line 72:		Line 72:
	<!--]]--></div>		<!--]]--></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>.		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>.		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">		:<source lang="matlab">
	[A,B,C] = ilp(J,N,O,s,r);		[A,B,C] = ilp(M,N,Q,s,r);
	</source>		</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==		==References==

Himpe: /* Usage */ remove emgr line

2018-01-10T11:32:07Z

Usage: remove emgr line

← Older revision		Revision as of 11:32, 10 January 2018
Line 80:		Line 80:
	[A,B,C] = ilp(J,N,O,s,r);		[A,B,C] = ilp(J,N,O,s,r);
	</source>		</source>

	~~'''ilp''' is compatible with [[wikipedia:MATLAB\|MATLAB]] and [[wikipedia:GNU_Octave\|OCTAVE]] and the matlab code can be downloaded from: [http://gramian.de/ilp.m ilp.m].~~
	~~The [[Emgr\|Empirical Gramian Framework]] can be obtained at [http://gramian.de/emgr.m http://gramian.de].~~

	==References==		==References==

Himpe: /* Usage */ updated Matlab code

2018-01-10T11:29:57Z

Usage: updated Matlab code

← Older revision		Revision as of 11:29, 10 January 2018
Line 26:		Line 26:
	function [A B C] = ilp(J,N,O,s,r)		function [A B C] = ilp(J,N,O,s,r)
	% ilp (inverse lyapunov procedure)		% ilp (inverse lyapunov procedure)
	% by Christian Himpe, 2013-~~2014 ( http://gramian.de )~~		% by Christian Himpe, 2013--2018
	% released under BSD 2-Clause License ~~( http://gramian.de/#license )~~		% released under BSD 2-Clause License
	%*		%*
			if(nargin==5)
	~~if(exist('emgr')~=2) disp('emgr framework is required. Download at http://gramian.de/emgr.m'); return; end~~		rand('seed',r);
			randn('seed',r);
	if(nargin==5) rand('seed',r); randn('seed',r); end;		end;

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

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

	% Balancing Transformation		% Balancing Transformation
	WC = Pdiag(WC)P';		WC = Pdiag(WC)P';
	WO = Qdiag(WO)Q';		WO = Qdiag(WO)Q';
	[U D V] = svd(WC*WO);		[U,D,V] = svd(WC*WO);

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

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

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

	% Solve System Matrix		% Solve System Matrix
	f = @(x,u,~~p) -D~~x+Bu;		A = -sylvester(D,D,B*B');
	~~g = @(x,u,p) C*x;~~
	~~A = -emgr(f,g,[J N O],[0 0.01 1],~~'~~c') - (1e-13)*eye(N~~);

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


			end
	</source>		</source>
	<!--]]--></div>		<!--]]--></div>

Himpe: /* Usage */ Code Update

2014-06-03T10:08:29Z

Usage: Code Update

← Older revision		Revision as of 10:08, 3 June 2014
Line 26:		Line 26:
	function [A B C] = ilp(J,N,O,s,r)		function [A B C] = ilp(J,N,O,s,r)
	% ilp (inverse lyapunov procedure)		% ilp (inverse lyapunov procedure)
	% by Christian Himpe, 2013 ( http://gramian.de )		% by Christian Himpe, 2013-2014 ( http://gramian.de )
	% released under BSD 2-Clause License ( http://gramian.de/#license )		% released under BSD 2-Clause License ( http://gramian.de/#license )
	%*		%*
Line 34:		Line 34:
	if(nargin==5) rand('seed',r); randn('seed',r); end;		if(nargin==5) rand('seed',r); randn('seed',r); end;

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

	%% Gramian Eigenvectors		% Gramian Eigenvectors
	X = randn(N,N);		[P S Q] = svd(randn(N));
	~~[U E V] = svd(X~~);

	%% Balancing ~~Trafo~~		% Balancing Transformation
	[P ~~D Q] = svd(~~diag(WC.*WO));		WC = Pdiag(WC)P';
	W = -D;		WO = Qdiag(WO)Q';
			[U D V] = svd(WC*WO);

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

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

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

	%% Solve System Matrix		% Solve System Matrix
	f = @(x,u,p) Wx+Bu;		f = @(x,u,p) -Dx+Bu;
	g = @(x,u,p) C*x;		g = @(x,u,p) C*x;
	A = -emgr(f,g,[J N O],[0 0.01 1],'c');		A = -emgr(f,g,[J N O],[0 0.01 1],'c') - (1e-13)*eye(N);

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

	~~%% Unbalance System~~
	~~T = U'*P';~~
	~~A = TAT';~~
	~~B = T*B;~~
	~~C = C*T';~~

	</source>		</source>