# balreal

Balanced state-space realization

## Syntax

``[sysb,g] = balreal(sys)``
``[sysb,g] = balreal(sys,Name=Value)``
``[sysb,g,TL,TR] = balreal(sys,___)``

## Description

````[sysb,g] = balreal(sys)` computes a balanced state-space realization of the LTI model `sys`. For stable models, `sys` is an equivalent realization for which the controllability and observability Gramians are equal and diagonal, their diagonal entries forming the vector `g` of Hankel singular values. This balances the input-to-state and state-to-output energy transfers, and small entries of `g` indicate states that you can remove with `xelim`.If `sys` has unstable poles, the function isolates the stable part, balances it, and adds back to the unstable part to form `sysb`. The entries of `g` corresponding to unstable modes are set to `Inf`.```

example

````[sysb,g] = balreal(sys,Name=Value)` computes the balanced realization using options specified as one or more name-value arguments. (since R2023b)```
````[sysb,g,TL,TR] = balreal(sys,___)` also returns the balancing state transformation matrices `TL` and `TR`. (since R2023b)```

example

## Examples

collapse all

Consider the following zero-pole-gain model, with near-canceling pole-zero pairs:

`sys = zpk([-10 -20.01],[-5 -9.9 -20.1],1)`
```sys = (s+10) (s+20.01) ---------------------- (s+5) (s+9.9) (s+20.1) Continuous-time zero/pole/gain model. ```

Obtain a state-space realization with balanced Gramians.

`[sysb,g] = balreal(sys);`

Examine the diagonal entries of the joint Gramians.

`g'`
```ans = 1×3 0.1006 0.0001 0.0000 ```

This indicates that the last two states of `sysb` are weakly coupled to the input and output. You can then delete these states using `xelim`.

`sysr = xelim(sysb,[2 3]);`

This yields a first-order approximation of the original system.

`zpk(sysr)`
```ans = -0.00011597 (s-8638) -------------------- (s+4.981) Continuous-time zero/pole/gain model. ```

Compare the Bode responses of the original and reduced-order models.

```bp = bodeplot(sys,sysr,'r--'); bp.PhaseMatchingEnabled = "on";```

The plots shows that removing the second and third states does not have much effect on system dynamics.

For faster and more accurate results, use `reducespec` for model reduction workflows.

This example shows how to perform balanced realization of a system with unstable poles.

Consider the following system.

`$\mathit{H}\left(\mathit{s}\right)=\frac{0.1}{\left(\mathit{s}+2\right)\left(\mathit{s}+1\right)\left(\mathit{s}-1\right)\left(\mathit{s}+0.001\right)}$`

The system contains three stable and one unstable poles. One stable pole is really close to the stability boundary.

Create the model.

`sys = ss(zpk(1,[-2 -1 1 -0.001],0.1));`

Apply `balreal` to create a balanced realization.

```[sysbal,g,TL,TR] = balreal(sys); g```
```g = 4×1 Inf 25.0374 0.0391 0.0018 ```

If `sys` has unstable poles, `balreal` isolates the stable part, balances the stable part, and adds back the unstable part to form the balanced realization. The unstable pole shows up as `Inf` in the vector `g`.

You can also use the offset option to shift the stability boundary to exclude nominally stable poles.

```[sysbal,g,TL,TR] = balreal(sys,Offset=1e-2); g```
```g = 4×1 Inf Inf 0.0393 0.0018 ```

For this realization, `balreal` treats the pole at `s = –0.001` as unstable.

## Input Arguments

collapse all

Dynamic system to decompose, specified as a numeric LTI model, such as a `ss` or `tf` model.

### Name-Value Arguments

Specify optional pairs of arguments as `Name1=Value1,...,NameN=ValueN`, where `Name` is the argument name and `Value` is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose `Name` in quotes.

Example: ```[sys,g2] = balreal(sys,Algorithm="relative",FreqIntervals=[10,100]);```

Balancing algorithm, specified as either `"absolute"`' or `"relative"`. Use this option to select the balancing algorithm and control the type of reduced-order model Gr.

• `"absolute"` — Minimize the absolute error ${‖G-{G}_{r}‖}_{\infty }$.

• `"relative"` — Minimize the relative error ${‖{G}^{-1}\left(G-{G}_{r}\right)‖}_{\infty }$.

Relative error gives a better match across frequency while absolute error emphasizes areas with most gain.

Regularization level, specified as `"auto"` or a nonnegative scalar value. When you set `Algorithm` to `"relative"`, `balreal` balances the model `[sys,α*I]` instead of `sys`. Set this option to `"auto"` to let `balreal` pick a suitable regularization level `α` value. Specify a value α ≥ 0 to override this default.

This regularization value helps provide a well-defined relative error at all frequencies.

Frequency intervals for computing frequency-limited Hankel singular values, specified as a matrix with two columns. Each row specifies a frequency interval `[fmin,fmax]`, where `fmin` and `fmax` are nonnegative frequencies, expressed in the frequency unit of the model. When identifying low-energy states to truncate, the software computes state contributions to system behavior in these frequency ranges only. For example:

• To restrict the computation to the range between 3 rad/s and 15 rad/s, assuming the frequency unit of the model is rad/s, set `FreqIntervals` to `[3,15]`.

• To restrict the computation to two frequency intervals, 3-15 rad/s and 40-60 rad/s, use `[3,15;40,60]`.

• To specify all frequencies below a cutoff frequency `fcut`, use `[0,fcut]`.

• To specify all frequencies above the cutoff, use `[fcut,Inf]` in continuous time, or `[fcut,pi/Ts]` in discrete time, where `Ts` is the sample time of the model.

The default value, `[]`, imposes no frequency limitation and is equivalent to `[0,Inf]` in continuous time or `[0,pi/Ts]` in discrete time. However, if you specify a `TimeIntervals` value other than `[]`, then this limit overrides `FreqIntervals = []`. If you specify both a `TimeIntervals` value and a `FreqIntervals` value, then the computation uses the union of these intervals.

Time intervals for computing time-limited Hankel singular values, specified as a matrix with two columns. Each row specifies a time interval `[tmin,tmax]`, where `tmin` and `tmax` are nonnegative times, expressed in the time unit of the model. When identifying low-energy states to truncate, the software computes state contributions to the system’s impulse response in these time intervals only. For example:

• To restrict the computation to the range between 3 s and 15 s, assuming the time unit of the model is seconds, set `TimeIntervals` to `[3,15]`.

• To restrict the computation to two time intervals, 3-15 s and 40-60 s, use `[3,15; 40,60]`.

• To specify all times from zero up to a cutoff time `tcut`, use `[0,tcut]`. To specify all times after the cutoff, use `[tcut,Inf]`.

The default value, `[]`, imposes no time limitation and is equivalent to `[0,Inf]`. However, if you specify a `FreqIntervals` value other than `[]`, then this limit overrides `Timeintervals = []`. If you specify both a `TimeIntervals` value and a `FreqIntervals` value, then the computation uses the union of these intervals.

Maximum accuracy loss factor in stable and unstable decomposition, specified as a positive scalar. For models with unstable poles, `balreal` first extracts the stable dynamics using the transformations `TL`, `TR`. Use `SepTol` to control the decomposition accuracy. Increasing `SepTol` helps separate nearby stable and unstable modes at the expense of accuracy. For more information, see `stabsep`.

Offset for the stable/unstable boundary, specified as a positive scalar value. In the stable/unstable decomposition, the stable term includes only poles satisfying

• `Re(s) < -Offset × max(1,|Im(s)|)` (continuous time)

• `|z| < 1 - Offset` (discrete time)

Increase the value of `Offset` to treat poles close to the stability boundary as unstable.

## Output Arguments

collapse all

Balanced state-space realization, returned as a state-space model.

Hankel singular values, returned as a vector.

The entries of `g` corresponding to unstable modes are set to `Inf`.

Since R2023b

Left-side matrix of the balancing state transformation, returned as a matrix with dimensions Nx-by-Nx, where Nx is the number of states in the model `sys`.

The algorithm transforms the state-space realization (A, B, C, D, E) of a model to (Ab, Bb, Cb, Db, Eb) given by:

• For explicit state-space models

`$\begin{array}{l}\begin{array}{cc}{A}_{b}={T}_{L}A{T}_{R},& {B}_{b}={T}_{L}B,\end{array}\\ \begin{array}{ccc}{C}_{b}=C{T}_{R},& {D}_{b}=D,& {E}_{b}={T}_{L}{T}_{R}=I\end{array}\end{array}$`
• For descriptor state-space models

`$\begin{array}{l}\begin{array}{cc}{A}_{b}={T}_{L}A{T}_{R},& {B}_{b}={T}_{L}B,\end{array}\\ \begin{array}{ccc}{C}_{b}=C{T}_{R},& {D}_{b}=D,& {E}_{b}={T}_{L}E{T}_{R}\end{array}\end{array}$`

The function returns an empty value `[]` for this argument when the input model `sys` is not a state-space model.

Since R2023b

Right-side matrix of the balancing state transformation, returned as a matrix with dimensions Nx-by-Nx, where Nx is the number of states in the model `sys`.The algorithm transforms the state-space realization (A, B, C, D, E) of a model to (Ab, Bb, Cb, Db, Eb) given by:

• For explicit state-space models

`$\begin{array}{l}\begin{array}{cc}{A}_{b}={T}_{L}A{T}_{R},& {B}_{b}={T}_{L}B,\end{array}\\ \begin{array}{ccc}{C}_{b}=C{T}_{R},& {D}_{b}=D,& {E}_{b}={T}_{L}{T}_{R}=I\end{array}\end{array}$`
• For descriptor state-space models

`$\begin{array}{l}\begin{array}{cc}{A}_{b}={T}_{L}A{T}_{R},& {B}_{b}={T}_{L}B,\end{array}\\ \begin{array}{ccc}{C}_{b}=C{T}_{R},& {D}_{b}=D,& {E}_{b}={T}_{L}E{T}_{R}\end{array}\end{array}$`

The function returns an empty value `[]` for this argument when the input model `sys` is not a state-space model.

## Tips

For model order reduction purposes, use `reducespec`.

## Algorithms

Consider the stable state-space model

`$\begin{array}{l}E\stackrel{˙}{x}=Ax+Bu\\ y=Cx+Du\end{array}$`

with controllability and observability Gramians Wc and Wo.

The state coordinate transformation $x={T}_{R}\overline{x}$ produces the equivalent model

`$\begin{array}{l}{T}_{L}E{T}_{R}\stackrel{˙}{\overline{x}}={T}_{L}A{T}_{R}\overline{x}+{T}_{L}Bu\\ y=C{T}_{R}\overline{x}+Du\end{array}$`

`balreal` transforms the Gramians to

`$\begin{array}{cc}{\overline{W}}_{c}={T}_{L}{W}_{c}{T}_{L}^{T},& {\overline{W}}_{o}={T}_{R}^{T}{W}_{o}\end{array}{T}_{R}$`

such that

`${\overline{W}}_{c}={\overline{W}}_{o}=diag\left(g\right)$`

If the model has unstable poles, the function isolates the stable part, balances it, and adds back to the unstable part to form the realization. The entries of `g` corresponding to unstable modes are set to `Inf`.

See [1] and [2] for details on the algorithm.

If you use the `TimeIntervals` or `FreqIntervals` options, then `balreal` bases the balanced realization on time-limited or frequency-limited controllability and observability Gramians. For information about calculating time-limited and frequency-limited Gramians, see `gram` and [4].

## References

[1] Laub, A., M. Heath, C. Paige, and R. Ward. “Computation of System Balancing Transformations and Other Applications of Simultaneous Diagonalization Algorithms.” IEEE Transactions on Automatic Control 32, no. 2 (February 1987): 115–22. https://doi.org/10.1109/TAC.1987.1104549.

[2] Moore, B. “Principal Component Analysis in Linear Systems: Controllability, Observability, and Model Reduction.” IEEE Transactions on Automatic Control 26, no. 1 (February 1981): 17–32. https://doi.org/10.1109/TAC.1981.1102568.

[3] Laub, Alan J. “Computation of ‘Balancing’ Transformations.” Joint Automatic Control Conference, no. 17 (1980): 84. https://doi.org/10.1109/JACC.1980.4232165.

[4] Gawronski, Wodek, and Jer-Nan Juang. “Model Reduction in Limited Time and Frequency Intervals.” International Journal of Systems Science 21, no. 2 (February 1990): 349–76. https://doi.org/10.1080/00207729008910366.

## Version History

Introduced before R2006a

expand all