stabsep

Stable-unstable decomposition

Syntax

[G1,G2] = stabsep(G)

[G1,G2] = stabsep(G,Offset)

[G1,G2] = stabsep(G,___,Name=Value)

[G1,G2,info] = stabsep(___)

Description

[G1,G2] = stabsep(G) decomposes the linear time-invariant (LTI) model G into its stable and unstable parts.

$$G=G_1+G_2$$

G1 contains all stable modes (poles) that can be separated from the unstable modes in a numerically stable way, and G2 contains the remaining modes. G1 is always proper and contains the feedthrough term (if any).

[G1,G2] = stabsep(G,Offset) shifts the stability boundary so that G1 includes all poles satisfying: (since R2023b)

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

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

[G1,G2] = stabsep(G,___,Name=Value) specifies additional options for controlling accuracy or focusing on unstable modes. (since R2023b)

[G1,G2,info] = stabsep(___) returns a structure info containing the block-diagonalizing transformation matrices. (since R2023b)

Examples

Stable-Unstable Decomposition for Zero-Pole-Gain Model

This example shows how to separate a zero-pole-gain model into its stable and unstable parts.

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 a zpk model of this system and obtain the stable-unstable decomposition.

h = zpk(1,[-2 -1 1 -0.001],0.1);
[h1,h2] = stabsep(h)

h1 =
 
           0.1
  ---------------------
  (s+0.001) (s+1) (s+2)
 
Continuous-time zero/pole/gain model.
Model Properties

h2 =
 
  -2.8326e-18
  -----------
     (s-1)
 
Continuous-time zero/pole/gain model.
Model Properties

The stable part h1 contains all the stable poles. You can use the Offset argument to shift the stability boundary and exclude the nominally stable pole (s = –0.001) from the stable part of the decomposition.

Compute the decomposition with the shifted boundary.

[h1,h2] = stabsep(h,0.1)

h1 =
 
  -0.050075 (s+2.999)
  -------------------
      (s+1) (s+2)
 
Continuous-time zero/pole/gain model.
Model Properties

h2 =
 
  0.050075 (s-1)
  ---------------
  (s+0.001) (s-1)
 
Continuous-time zero/pole/gain model.
Model Properties

For this decomposition, the nominally stable pole is in the unstable part h2.

By default, the function provides the stable part in the first output argument and the unstable part in the second output argument. You can use the Focus option to reverse the order.

[h1,h2] = stabsep(h,Focus="unstable")

h1 =
 
  -2.1943e-18
  -----------
     (s-1)
 
Continuous-time zero/pole/gain model.
Model Properties

h2 =
 
           0.1
  ---------------------
  (s+0.001) (s+1) (s+2)
 
Continuous-time zero/pole/gain model.
Model Properties

For this decomposition, h1 contains the unstable pole and the nominally stable pole.

Input Arguments

G — Dynamic system to decompose

numeric LTI model

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

Offset — Offset for stable/unstable boundary

0 (default) | scalar

Since R2023b

Offset for stable/unstable boundary, specified as a scalar. 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 exclude poles close to the stability boundary from G1.

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.

Example: [G1,G2] = stabsep(G,0.1,SepTol=100)

Focus — Focus of decomposition

"stable" (default) | "unstable"

Since R2023b

Focus of decomposition, specified as either "stable" or "unstable".

• "stable" — The first argument G1 includes all poles satisfying

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

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

• "unstable" — The first argument G1 includes all poles satisfying

  • Re(s) > Offset × max(1,|Im(s)|) (continuous time)

  • |z| > 1 + Offset (discrete time)

When an exact stable/unstable split is numerically impossible, G1 includes a maximum subset of the desired modes and G2 contains the remaining mix of stable and unstable modes.

SepTol — Accuracy loss factor

10 (default) | positive scalar

Since R2023b

Accuracy loss factor for stable/unstable decomposition, specified as a positive scalar value. When decomposing a model G(s), stabsep ensures that the frequency responses of G and G₁ + G₂ differ by no more than SepTol times the absolute accuracy of the computed value of G(s). Increasing SepTol helps separate nearby stable and unstable modes at the expense of accuracy.

Output Arguments

G1 — Stable part of decomposition

numeric LTI model

Stable part of the decomposition, returned as a model of same type as G.

When Focus is set to "unstable", the function returns the unstable part of the decomposition in the G1 argument.

G2 — Unstable part of decomposition

numeric LTI model

Unstable part of the decomposition, returned as a model of same type as G.

When Focus is set to "unstable", the function returns the stable part of the decomposition in the G2 argument.

info — Additional information

structure | []

Since R2023b

Additional information about the decomposition, returned as structure with these fields.

Field	Description
TL	Left-side matrix of the block-diagonalizing transformation, returned as a matrix with dimensions Nx-by-Nx, where Nx is the number of states in the model G.
TR	Right-side matrix of the block-diagonalizing transformation, returned as a matrix with dimensions Nx-by-Nx, where Nx is the number of states in the model G.

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

$$\begin{array}{ccc}{T}_{L}A{T}_R=\left(\begin{array}{cc}{A}_1& 0\\ 0& A_2\end{array}\right),& T_{L}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.

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

References

[1] Bartels, R. H., and G. W. Stewart. “Algorithm 432 [C2]: Solution of the Matrix Equation AX + XB = C [F4].” Communications of the ACM 15, no. 9 (September 1972): 820–26. https://doi.org/10.1145/361573.361582.

[2] Bavely, Connice A., and G. W. Stewart. “An Algorithm for Computing Reducing Subspaces by Block Diagonalization.” SIAM Journal on Numerical Analysis 16, no. 2 (April 1979): 359–67. https://doi.org/10.1137/0716028.

[3] Demmel, James. “The Condition Number of Equivalence Transformations That Block Diagonalize Matrix Pencils.” SIAM Journal on Numerical Analysis 20, no. 3 (June 1983): 599–610. https://doi.org/10.1137/0720040.

Version History

Introduced before R2006a

R2023b: Use simplified syntax

Starting in R2023b, you can specify Offset directly as the second input argument and the Focus and SepTol options as name-value arguments. For example:

[G1,G2] = stabsep(G,0.1,Focus="unstable");

As a result of this change, stabsepOptions and the syntax [G1,G2] = stabsep(G,options) are not recommended.