bstmr

(Not recommended) Balanced stochastic model truncation (BST) via Schur method

bstmr is not recommended. Use reducespec instead. (since R2023b) For more information on updating your code, see Version History.

Syntax

GRED = bstmr(G)
GRED = bstmr(G,order)
[GRED,redinfo] = bstmr(G,key1,value1,...)
[GRED,redinfo] = bstmr(G,order,key1,value1,...)

Description

bstmr returns a reduced order model GRED of G and a struct array redinfo containing the error bound of the reduced model and Hankel singular values of the phase matrix of the original system [2].

The error bound is computed based on Hankel singular values of the phase matrix of G. For a stable system these values indicate the respective state energy of the system. Hence, reduced order can be directly determined by examining these values.

With only one input argument G, the function will show a Hankel singular value plot of the phase matrix of G and prompt for model order number to reduce.

This method guarantees an error bound on the infinity norm of the multiplicative ∥ GRED–1(G-GRED) ∥ ∞ or relative error ∥ G^-–1(G-GRED) ∥ ∞ for well-conditioned model reduction problems [1]:

$${\Vert G^{-1}(G-Gred)\Vert }_{\infty }\le {\displaystyle \prod _{k+1}^n\left(1+2{\sigma }_i(\sqrt{1+{\sigma }_i^2}+{\sigma }_i)\right)}-1$$

This table describes input arguments for bstmr.

Argument	Description
G	LTI model to be reduced (without any other inputs will plot its Hankel singular values and prompt for reduced order)
ORDER	(Optional) an integer for the desired order of the reduced model, or a vector of desired orders for batch runs

A batch run of a serial of different reduced order models can be generated by specifying order = x:y, or a vector of integers. By default, all the anti-stable part of a system is kept, because from control stability point of view, getting rid of unstable state(s) is dangerous to model a system.

'MaxError' can be specified in the same fashion as an alternative for 'ORDER'. In this case, reduced order will be determined when the accumulated product of Hankel singular values shown in the above equation reaches the 'MaxError'.

Argument	Value	Description
'MaxError'	Real number or vector of different errors	Reduce to achieve H∞ error. When present, 'MaxError' overrides ORDER input.
'Display'	'on' or 'off'	Display Hankel singular plots (default 'off').
'Order'	Integer, vector or cell array	Order of reduced model. Use only if not specified as 2nd argument.

This table describes output arguments.

Argument	Description
GRED	LTI reduced order model. Become multi-dimension array when input is a serial of different model order array.
REDINFO	A STRUCT array with three fields: REDINFO.ErrorBound (bound on ∥G–1(G-GRED) ∥∞) REDINFO.StabSV (Hankel SV of stable part of G) REDINFO.UnstabSV (Hankel SV of unstable part of G)

G can be stable or unstable, continuous or discrete.

Examples

Given a continuous or discrete, stable or unstable system, G, the following commands can get a set of reduced order models based on your selections:

rng(1234,'twister'); 
G = rss(30,5,4);
G.D = zeros(5,4);
[g1, redinfo1] = bstmr(G); % display Hankel SV plot 
                           % and prompt for order (try 15:20)
[g2, redinfo2] = bstmr(G,20); 
[g3, redinfo3] = bstmr(G,[10:2:18]);
[g4, redinfo4] = bstmr(G,'MaxError',[0.01, 0.05]);
for i = 1:4
    figure(i)
	eval(['sigma(G,g' num2str(i) ');']);
end

Algorithms

Given a state space (A,B,C,D) of a system and k, the desired reduced order, the following steps will produce a similarity transformation to truncate the original state-space system to the k^th order reduced model.

1. Find the controllability Gramian P and observability Gramian Q of the left spectral factor Φ = Γ(σ)Γ*(–σ) = Ω*(–σ)Ω(σ) by solving the following Lyapunov and Riccati equations

   AP + PA^T + BB^T = 0

   B_W = PC^T + BD^T

   QA + A^T Q + (QB_W – C^T) (–DD^T) (QB_W – C^T)^T = 0

2. Find the Schur decomposition for PQ in both ascending and descending order, respectively,

   $$\begin{array}{l}{V}_A^{T}PQ{V}_A=\left[\begin{array}{ccc}{\lambda }_1& \cdots & \cdots \\ 0& \cdots & \cdots \\ 0& 0& {\lambda }_n\end{array}\right]\\ V_D^{T}PQ{V}_D=\left[\begin{array}{ccc}{\lambda }_n& \cdots & \cdots \\ 0& \cdots & \cdots \\ 0& 0& {\lambda }_1\end{array}\right]\end{array}$$

3. Find the left/right orthonormal eigen-bases of PQ associated with the k^th big Hankel singular values of the all-pass phase matrix (W^*(s))^–1G(s).

                                         k

   $$\begin{array}{l}{V}_A=[V_{R,SMALL},\stackrel{k}{\overbrace{V_{L,BIG}}}]\\ V_D=[\stackrel{}{\overbrace{V_{R,BIG}}},V_{L,SMALL}]\end{array}$$

4. Find the SVD of (V^T _L,BIGV_R,BIG) = U Σ ςΤ

5. Form the left/right transformation for the final k^th order reduced model

   S_L,BIG = V_L,BIG U Σ(1:k,1:k)^–½

   S_R,BIG = V_R,BIG V Σ(1:k,1:k)^–½

6. Finally,

   $$\left[\begin{array}{cc}\widehat{A}& \widehat{B}\\ \widehat{C}& \widehat{D}\end{array}\right]=\left[\begin{array}{cc}{S}_{L,BIG}^{T}A{S}_{R,BIG}& S_{L,BIG}^{T}B\\ C{S}_{R,BIG}& D\end{array}\right]$$

The proof of the Schur BST algorithm can be found in [1].

Note

The BST model reduction theory requires that the original model D matrix be full rank, for otherwise the Riccati solver fails. For any problem with strictly proper model, you can shift the jω-axis via bilin such that BST/REM approximation can be achieved up to a particular frequency range of interests. Alternatively, you can attach a small but full rank D matrix to the original problem but remove the D matrix of the reduced order model afterwards. As long as the size of D matrix is insignificant inside the control bandwidth, the reduced order model should be fairly close to the true model. By default, the bstmr program will assign a full rank D matrix scaled by 0.001 of the minimum eigenvalue of the original model, if its D matrix is not full rank to begin with. This serves the purpose for most problems if user does not want to go through the trouble of model pretransformation.

References

[1] Zhou, K., “Frequency-weighted model reduction with L∞ error bounds,” Syst. Contr. Lett., Vol. 21, 115-125, 1993.

[2] Safonov, M.G., and R.Y. Chiang, “Model Reduction for Robust Control: A Schur Relative Error Method,” International J. of Adaptive Control and Signal Processing, Vol. 2, p. 259-272, 1988.

Version History

Introduced before R2006a

R2023b: Not recommended

bstmr is not recommended. To perform balanced truncation model order reduction, use the reducespec function with the following syntax.

R = reducespec(sys,"balanced");

For the full workflow, see Task-Based Model Order Reduction Workflow.