BalancedTruncation

Balanced truncation model order reduction

Since R2023b

Description

The BalancedTruncation object stores model order reduction (MOR) specifications for the balanced truncation of ordinary (nonsparse) linear time-invariant (LTI) models.

Creation

The reducespec function creates a balanced truncation MOR object when you use this syntax.

R = reducespec(sys,"balanced")

Here, sys is any nonsparse LTI model. The workflow uses this object to set up MOR tasks and store results. For the full workflow, see Task-Based Model Order Reduction Workflow.

Properties

Sigma — Hankel singular values

[] (default) | vector

This property is read-only.

Hankel singular values $${\sigma }_j,j=1\dots N$$, returned as a vector of size N-by-1. Here, N is the number of states in the model.

In state coordinates that equalize the input-to-state and state-to-output energy transfers, the Hankel singular values (HSVs) measure the contribution of each state to the input/output behavior. Hankel singular values relate to model order as singular values relate to matrix rank. In particular, small HSVs indicate states that you can discard to simplify the model.

Energy — Normalized state energy

[] (default) | vector

This property is read-only.

Normalized state energy, returned as a vector of size N-by-1. Here, N is the number of states in the model.

These values measure the energy of each state relative to the state with maximum energy. The normalized state energy for the kth HSV is given by $$\frac{{\sigma }_k}{{\sigma }_1}$$.

In balanced coordinates, the Hankel singular values σ are the eigenvalues of the equalized Gramians X_r = X_o. From

$$X_r={\displaystyle \underset{0}{\overset{\infty }{\int }}x(t)x{(t)}^T}dt$$

where $$x(t)=e^{tA}B$$ is the impulse response, you obtain the total state energy as follows:

$${\displaystyle \sum {\sigma }_i}=\text{tr}(X_r)={\displaystyle \underset{0}{\overset{\infty }{\int }}{\Vert x(t)\Vert }^{2}dt}$$

Hence, σ_k is the energy of kth principal state and $$\frac{{\sigma }_k}{{\sigma }_1}$$ is the normalized energy of kth principal state.

Error — Approximation error bound

[] (default) | vector

This property is read-only.

Bound on absolute or relative approximation error, returned as a vector of size N-by-1. Here, N is the number of states in the model. For more information, see Algorithms.

Loss — Neglected fraction of state energy

[] (default) | vector

Since R2024b

This property is read-only.

Neglected fraction of state energy, returned as a vector of size N-by-1. Here, N is the number of states in the model.

For example, when you select an order r, the neglected fraction of total energy is given by:

$$\frac{{\displaystyle \sum _{i=r+1}^N{\sigma }_i}}{{\displaystyle \sum _{j=1}^N{\sigma }_j}}$$

TL — Left-side matrix of transformation

[] (default) | matrix

This property is read-only.

Left-side matrix of transformation, returned as a matrix of size N-by-N. Here, N is the number of states in the model.

The balanced truncation algorithm transforms the state-space realization (A, B, C) of a model to

$$\begin{array}{ccc}{T}_L^{T}A{T}_R=\left(\begin{array}{cc}{A}_1& 0\\ 0& A_2\end{array}\right),& T_L^{T}B=\left(\begin{array}{c}{B}_1\\ B_2\end{array}\right),& C{T}_R=\left(\begin{array}{cc}{C}_1& C_2\end{array}\right).\end{array}$$

The transformation scales the system and separates the stable and unstable parts. All modes of A₁ are unstable and all modes of A₂ are stable.

TR — Right-side matrix of transformation

[] (default) | matrix

This property is read-only.

Right-side matrix of transformation, returned as a matrix of size N-by-N. Here, N is the number of states in the model.

The algorithm transforms the state-space realization (A, B, C) of a model to

$$\begin{array}{ccc}{T}_L^{T}A{T}_R=\left(\begin{array}{cc}{A}_1& 0\\ 0& A_2\end{array}\right),& T_L^{T}B=\left(\begin{array}{c}{B}_1\\ B_2\end{array}\right),& C{T}_R=\left(\begin{array}{cc}{C}_1& C_2\end{array}\right).\end{array}$$

The transformation scales the system and separates the stable and unstable parts. All modes of A₁ are unstable and all modes of A₂ are stable.

Lr — Cholesky factor of controllability Gramian

[] (default) | matrix

This property is read-only.

Cholesky factor of the controllability Gramian of the stable part, returned as a matrix of size N-by-N. Here, N is the number of states in the model.

The controllability Gramian of the stable part (A₂, B₂, C₂) is $$X_r=L_r{L}_r^T$$.

Lo — Cholesky factor of observability Gramian

[] (default) | matrix

This property is read-only.

Cholesky factor of the observability Gramian of the stable part, returned as a matrix of size N-by-N. Here, N is the number of states in the model.

The observability Gramian of the stable part (A₂, B₂, C₂) is $$X_o=L_o{L}_o^T$$.

Regularization — Regularization level

[] (default) | matrix

This property is read-only.

Regularization level when R.Options.Algorithm is set to "relative", returned as a scalar value.

Options — Options for balanced truncation

BalancedTruncationOptions object

Options for balanced truncation of LTI models, specified as a BalancedTruncationOptions object. Use dot notation to configure options for R. For example R.Options.Algorithm = "relative".

For more information about available options, see BalancedTruncationOptions.

Object Functions

process	Run model order reduction algorithm
view (balanced)	Plot state contributions when using balanced truncation method
getrom (balanced)	Obtain reduced-order models when using balanced truncation method

Examples

Balanced Truncation of LTI Model in Particular Frequency Band

This example shows how to obtain a reduced-order model for a linear time-invariant (LTI) model using the balanced truncation method. In this example, you reduce a high-order model with a focus on the dynamics in a particular frequency range.

Load a model and examine its frequency response.

load('highOrderModel.mat','G')
bodeplot(G)

G is a 48th-order model with several large peak regions around 5.2 rad/s, 13.5 rad/s, and 24.5 rad/s, and smaller peaks scattered across many frequencies. Suppose that for your application you are only interested in the dynamics near the second large peak, between 10 rad/s and 22 rad/s. Focus the model reduction on the region of interest to obtain a good match with a low-order approximation.

Create a model order reduction task and specify the desired frequency interval.

R = reducespec(G,"balanced");
R.Options.FreqIntervals = [10,22];

Analyze the model and compute the derived information.

R = process(R);

Use the getrom function to retain all the modes between 10 rad/s and 22 rad/s.

[rsys,info] = getrom(R,Order=[10,18],Method="truncate");

Examine the frequency responses of the reduced-order models. Also, examine the difference between those responses and the original response (the absolute error).

subplot(2,1,1);
bodemag(G,rsys(:,:,1,1),rsys(:,:,2,1),logspace(0.5,1.5,100))
title('Bode Magnitude Plot')
legend("Original","Order 10", "Order 18")
subplot(2,1,2);
bodemag(G-rsys(:,:,1,1),G-rsys(:,:,2,1),logspace(0.5,1.5,100));
title('Absolute Error Plot')
legend('Order 10','Order 18');

With the frequency-limited energy computation, even the 10th-order approximation is quite good in the region of interest.

Reduce Order of LTI Model Using Balanced Truncation

This example shows how to create a model order reduction specification for a dense LTI model using the balanced truncation method.

For this example, generate a random discrete-time state-space model with 40 states.

rng(0)
sys = drss(40);

Plot the Bode response.

bode(sys)

Create the specification object.

R = reducespec(sys,"balanced");

Notice that reducespec does not perform any computation and creates only the object. This allows you to set additional options before running the model order reduction algorithm, such as relative error control as the algorithm objective.

R.Options.Algorithm = "relative";

For balanced truncation, you can visualize the contributions in terms of the Hankel singular values or normalized energies. By default, the view function plots Hankel singular values.

view(R)

For this example, select an order of 15 since it is the first order with a relative error less than 1e-4. In general, you select the order based on the desired absolute or relative fidelity. Then, compute the reduced-order model.

rsys = getrom(R,Order=15);

Plot the Bode response of both models.

bode(sys,rsys,'r--')
legend("Original","Order 15")

Frequency-Weighted Balanced Truncation

This example shows how to improve the balanced truncation approximation with the help of frequency weights to emphasize accuracy in a particular frequency band.

For this example, generate a random state-space model with 30 states, two outputs, and three inputs.

rng(2)
sys = rss(30,2,3);
size(sys)

State-space model with 2 outputs, 3 inputs, and 30 states.

Create a model order reduction task.

R = reducespec(sys,"balanced");

Obtain the reduced order model.

rsys = getrom(R,Order=19);
sigma(sys,sys-rsys,"r--")
legend("sys","sys-rsys");

This reduced-order model does not provide a great approximation in the high frequencies. To improve the response, you can use the options to specify frequency weights to emphasize accuracy in a particular frequency band. These weights should have high gain in frequency bands of interest and low gain elsewhere.

Specify the frequency weights with a high gain in the 3 rad/s to 100 rad/s band.

wi = tf([1 0 0 0],[1 6 18 27])*eye(3);
wo = tf(1,[0.01 1])*eye(2);
R.Options.InputWeight = wi;
R.Options.OutputWeight = wo;

Obtain the reduced-order model.

rsysFW = getrom(R,Order=19);
sigma(sys,sys-rsys,"r--",sys-rsysFW,"k-.")
legend("sys","rsys","rsysFW");

The model with frequency weighted reduction provides a better approximation across the higher frequencies.

Use Input-Output Scaling to Emphasize Accuracy in Particular Channels

This example shows how you can use input and output scaling to emphasize accuracy of model order reduction in a particular channel of a MIMO model.

Load the CD player state-space model.

load cdromData.mat cdrom
size(cdrom)

State-space model with 2 outputs, 2 inputs, and 120 states.

Create a model order reduction task.

R = reducespec(cdrom,"balanced");

Specify the input and output weights as static values such that the (1,2) channel is dominant.

R.Options.InputWeight = diag([1e-3 1]); 
R.Options.OutputWeight=diag([1 0.1]);

Obtain the reduced-order model.

rsys = getrom(R,MaxError=1e-2,Method="truncate");
bode(cdrom,rsys)

The reduced-order model provides a good match for the scaled I/O channel.

Model order reduction method controls the overall error of the MIMO model, so it is unlikely you get high relative accuracy on all channels when the gains are very different. Scaling an I/O channel helps you control the accuracy for that particular channel.

Algorithms

Balanced truncation first decomposes the system G into its stable and unstable parts:

$$G=G_s+G_u$$

When you specify R.Options.Algorithm as "absolute", the software uses the balanced truncation method of [1] to reduce G_s. This computes the HSVs σ_j based on the controllability and observability Gramians. For order r, the absolute error $${\Vert G_s-G_r\Vert }_{\infty }$$ is bounded by $$2{\displaystyle \sum _{j=r+1}^n{\sigma }_j}$$. Here, n is the number of states in G_s.

When you specify R.Options.Algorithm as "relative", the software uses the balanced stochastic truncation method of [2] to reduce G_s. For square G_s, this computes the HSVs σ_j of the phase matrix $$F={\left(W\text{'}\right)}^{-1}G$$, where W(s) is a stable, minimum-phase spectral factor of GG’:

$$W\text{'}(s)W(s)=G(s)G\text{'}(s)$$

For order r, the relative error $${\Vert G_s{}^{-1}(G_s-G_r)\Vert }_{\infty }$$ is bounded by

$${\displaystyle \prod _{j=r+1}^H\left(\frac{1+{\sigma }_j}{1-{\sigma }_j}\right)}-1\approx 2{\displaystyle \sum _{j=r+1}^n{\sigma }_j}$$

where $$2{\displaystyle \sum _{j=r+1}^n{\sigma }_j}\ll 1$$.

References

[1] Varga, A., "Balancing-Free Square-Root Algorithm for Computing Singular Perturbation Approximations," Proc. of 30th IEEE CDC, Brighton, UK (1991), pp. 1062-1065.

[2] Green, M., "A Relative Error Bound for Balanced Stochastic Truncation", IEEE Transactions on Automatic Control, Vol. 33, No. 10, 1988

Version History

Introduced in R2023b

R2024b: Select model order using new criteria based on total neglected state energy

You can now select model order based on the fraction of total energy you want to neglect. The BalancedTruncation and SparseBalancedTruncation objects now contain a new property Loss to store the energy loss as a fraction of total energy. You can use view(R,"loss") to visualize a bar chart of neglected fraction of total energy and getrom with the new name-value argument MaxLoss to pick the reduced order based on maximum energy loss as a fraction of total energy. For example, to obtain a reduced model that neglects 1% of the total energy, you can specify rsys = getrom(R,MaxLoss=0.01).