modalsep

Compute modal decomposition

Since R2023b

Syntax

[H,H0] = modalsep(G)

[H,H0] = modalsep(G,N,regionFcn)

[H,H0] = modalsep(___,Name=Value)

[H,H0,info] = modalsep(___)

Description

[H,H0] = modalsep(G) computes the modal decomposition for a linear time-invariant (LTI) system G and returns the modal components as a state-space array H and the static gain H0.

$$G(s)=H_0+{\displaystyle \sum _{j=1}^m{H}_j(s)}$$

Each modal component in H_j(s) is associated with a single real pole, a pair of complex conjugate poles, or a cluster repeated poles.

For sparse models, this syntax returns a truncated modal form. By default, the function computes up to 1000 modal components associated with the poles of smallest magnitude. (since R2026a)

[H,H0] = modalsep(G,N,regionFcn) computes the region-based modal decomposition

$$G(s)=H_0+{\displaystyle \sum _{j=1}^N{H}_j(s)}$$

Here, the modal components H_j(s) have their poles in disjoint regions R_j of the complex plane. N specifies the number of regions and regionFcn is the name or a handle to the function that specifies the partition into N regions.

[H,H0] = modalsep(___,Name=Value) returns the modal decomposition based on the options specified by one or more name-value arguments. Use these options to control the granularity and accuracy of the decomposition.

For sparse models, this syntax returns a subset of modal components based on the options specified by one or more name-value arguments. This subset is controlled by the Focus and MaxOrder options. (since R2026a)

[H,H0,info] = modalsep(___) also returns a structure info containing the block-diagonalizing transformation matrices for the model and mode information for each modal component.

Examples

Compute Modal Decomposition of LTI Models

This example shows how to obtain a modal decomposition for a linear time-invariant (LTI) model using modalsep.

Create a random MIMO state-space model with 20 states.

rng(0)
G = rss(20,2,3);

Obtain the modal decomposition of this model.

[H,H0] = modalsep(G);

Examine the size of H.

size(H)

16x1 array of state-space models.
Each model has 2 outputs, 3 inputs, and between 1 and 2 states.

Typically, the modal components are of order 1 or 2, but can have higher orders in case of cluster of repeated poles.

For this model, modalsep returns the static gain H0 for each I/O pair.

H0

H0 =
 
  D = 
            u1       u2       u3
   y1        0   -1.237  -0.3334
   y2   0.1554   -2.193        0
 
Static gain.
Model Properties

Additionally, you can obtain the original model back from the modal decomposition using modalsum.

G2 = modalsum(H,H0);
sigma(G,G2,'r--')

Region-Based Modal Decomposition

This example shows how to perform region-based modal decomposition of a state-space model. In this example, you perform the modal decomposition of a high-order model based on damping ratio of the poles.

Load a model and examine the damping ratio and natural frequency of the poles.

load('highOrderModel.mat','G')
[wn,zeta] = damp(G);

Visualize the damping ratios and natural frequencies.

semilogy(zeta,wn,"x");
grid on
title("Mode Damping and Natural Frequency")
ylabel("Natural frequency")
xlabel("Damping")

Based on this, you can define regions function that decomposes the model based on the damping ratio values. For example, this model has damping ratios between 0.023 and 0.05, you can divide three regions as follows:

• Region 1: $0.023\le z<0.033$

• Region 2: $0.033\le z<0.040$

• Region 3: $z\ge 0.040$

The region function to assign region index to these values is defined at the end of this example. See Region Function Definition.

Perform the modal decomposition.

[h,h0,info] = modalsep(G,3,@regionFcn);

Examine the size of each region.

size(h(:,:,1))

State-space model with 1 outputs, 1 inputs, and 24 states.

size(h(:,:,2))

State-space model with 1 outputs, 1 inputs, and 12 states.

size(h(:,:,3))

State-space model with 1 outputs, 1 inputs, and 12 states.

For this decomposition, region 1 has 24 poles in the specified range. Region 2 and 3 have 12 poles each in their specified ranges.

Region Function Definition

function IR = regionFcn(p)
    [~,z] = damp(p);
    IR = zeros(size(z));
    IR(z<0.033 & z>=0.023) = 1;
    IR(z<0.04 & z>=0.033) = 2;
    IR(z>=0.04) = 3;
end

Compute Truncated Modal Decomposition of Sparse Models

Since R2026a

This example shows how to obtain a truncated modal decomposition for a sparse model using modalsep.

Load the sparse model.

load flowmeterSparse.mat
size(sys)

Sparse state-space model with 5 outputs, 1 inputs, and 9669 states.

The sparse model contains 9669 states. By default, modalsep computes the first 1000 modal components associated with the poles of smallest magnitude for sparse models, which may be time consuming. Computing the modal decomposition, you can see that modalsep returns an array containing 1000 modal components.

[H,H0,info] = modalsep(sys);
size(H)

1000x1 array of state-space models.
Each model has 5 outputs, 1 inputs, and 1 states.

You can limit the subset of computed modal components using the Focus and MaxOrder options. Compute the modal components in the frequency range 0 rad/s to 250 rad/s and perform region-based modal decomposition to partition modal components in the complex plane. For example, split modal components into two regions:

• Region 1 — Poles with absolute value greater than 100.

• Region 2 — Poles with absolute value less than 100.

regionFcn = @(p) 1+(abs(p)<100);
N = 2;
[Hr,H0r,infor] = modalsep(sys,N,regionFcn,Focus=[0 250]);
size(Hr)

2x1 array of state-space models.
Each model has 5 outputs, 1 inputs, and between 13 and 20 states.

The function now returns a decomposition with two modal components, with 33 modes split into the specified regions.

Additionally, you can use modalsum to combine the components and obtain a lower-order approximation of the sparse model. The accuracy depends on the number of computed modal components.

sysm = modalsum(H,H0);
sysm2 = modalsum(Hr,H0r);
sigmaplot(sys,sysm,".-",sysm2,"--",w)
legend("Original model","Approximation (1000 poles)",...
    "Approximation (33 poles)")

Use this workflow to quickly obtain a truncated modal decomposition of a sparse model. For more flexibility, first use reducespec to obtain a reduced-order model, then apply modalsep to the result.

Input Arguments

G — Linear time-invariant model

ss | tf | zpk | sparss | mechss

Dynamic system model, specified as:

• Ordinary LTI model (ss, tf, or zpk).

• Sparse LTI model (sparss or mechss). (since R2026a)

  For sparse models, modalsep is helpful when you want to quickly obtain a truncated modal decomposition. For more flexibility, first use reducespec to obtain a reduced-order model, then apply modalsep to the result.

N — Number of regions

positive scalar

Number of regions, specified as a positive scalar.

Specify this value as the number of partitions performed by the regionFcn.

regionFcn — Region function

function name | function handle

User-defined function for partitioning the complex plane, specified as a function name (string) or a function handle.

Specify a region function of the form:

IR = regionFcn(p);

The function assigns a region index IR between 1 and N to a given pole p. To specify additional inputs to the region function, use an anonymous function.

regionFcn = @(p) myFcn(p,param1,...,paramK);

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: [H,H0] = modalsep(G,InputScaling=[1,1e-3])

Focus — Frequency range of interest

[0 inf] (default) | vector

Since R2026a

Frequency range of interest, specified as a vector of form [fmin,fmax]. When you specify a value for this property, the algorithm computes only modes in this frequency range.

By default, the focus is unspecified ([0 inf]) and the algorithm computes the MaxOrder modes with smallest magnitude.

MaxOrder — Maximum order of modal approximation

1000 (default) | positive integer

Since R2026a

Maximum order of the modal approximation, specified as a positive integer. This value limits the number of eigenvalues computed by the Krylov-Schur iterations ([1]) and the order of the modal approximation of the original sparse model.

ModeTol — Accuracy of computed modes

1e-12 (default) | positive scalar

Since R2026a

Tolerance for identifying converged eigenvalues in Krylov-Schur iterations, specified as a positive scalar.

DCFrequency — Frequency for evaluating DC contribution

0 (default) | nonnegative scalar

Frequency for evaluating and matching DC contributions, specified as a nonnegative scalar.

For models with integrators, you cannot evaluate modal contributions at DC since the DC gain is infinite. To evaluate modal contributions and match gains at a different frequency, set the property to a positive value. The default value of this property corresponds to the true DC value.

InputScaling — Input scaling

[] (default) | vector

Input scaling factors, specified as a vector of length Nu, where Nu is the number of inputs in the original model sys.

Use this option to emphasize specific input channels in sys. The software evaluates the modal contributions for the scaled system (see the info argument).

OutputScaling — Output scaling

[] (default) | vector

Output scaling factors, specified as a vector of length Ny, where Ny is the number of outputs in the original model sys.

Use this option to emphasize specific output channels in sys. The software evaluates the modal contributions for the scaled system (see the info argument).

SepTol — Relative accuracy of modal decomposition

1e-12 (default) | scalar between 0 and 1

Relative accuracy of modal decomposition, specified as a scalar between 0 and 1.

This option limits the condition number of the block diagonalizing transformation to roughly SepTol/eps. Increasing SepTol helps yield smaller modal components at the expense of accuracy.

For region-based decompositions, the function ignores this option because the block sized are predefined.

Display — Show or hide progress report

"on" (default) | "off"

Since R2026a

Show or hide progress report, specified as either "off" or "on".

Output Arguments

H — Modal components

model array

Modal components, returned as an array of models of the same type as G. For example, if G is a zpk model, modalsep returns the modal components as an array of zpk models.

For sparse state-space models, modalsep returns the modal components as an array of ss models.

Use H(:,:,j) to retrieve the submodel H_j(s).

H0 — Static gain

matrix | ss object

Static gain, returned as a matrix or a static ss object.

info — Additional information

structure | []

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

Field	Description
Mode	Average mode (pole) value in the modal component Hj, returned as a vector of length nc-by-1, where nc is the number of computed modal components.
DCGain	Normalized DC contribution of modal components, returned as a vector of length nc-by-1, where nc is the number of computed modal components.
TL	Left-side matrix of the block-diagonalizing transformation, returned as a matrix. For ordinary LTI models, TL is a matrix of size Nx-by-Nx, where Nx is the number of states in the model G For sparse models, TL is a matrix of size nc-by-Nx, where nc is the number of computed modal components. For mechss models, TL corresponds to the equivalent first-order representation sparss(G). (since R2026a) The algorithm transforms the state-space realization (A, B, C, D, E) of a model to block diagonal matrices (Am, Bm, Cm, Dm, Em) given by: For explicit state-space models $$\begin{array}{l}\begin{array}{cc}{A}_m=T_{L}A{T}_R,& B_m=T_{L}B,\end{array}\\ \begin{array}{ccc}{C}_m=C{T}_R,& D_m=D,& E=T_L{T}_R=I\end{array}\end{array}$$ For descriptor state-space models $$\begin{array}{l}\begin{array}{cc}{A}_m=T_{L}A{T}_R,& B_m=T_{L}B,\end{array}\\ \begin{array}{ccc}{C}_m=C{T}_R,& D_m=D,& E_m=T_{L}E{T}_R\end{array}\end{array}$$
TR	Right-side matrix of the block-diagonalizing transformation, returned as a matrix. For ordinary LTI models, TR is a matrix of size Nx-by-Nx, where Nx is the number of states in the model G For sparse models, TR is a matrix of size Nx-by-nc, where nc is the number of computed modal components. For mechss models, TR corresponds to the equivalent first-order representation sparss(G). (since R2026a) The algorithm transforms the state-space realization (A, B, C, D, E) of a model to block diagonal matrices (Am, Bm, Cm, Dm, Em) given by: For explicit state-space models $$\begin{array}{l}\begin{array}{cc}{A}_m=T_{L}A{T}_R,& B_m=T_{L}B,\end{array}\\ \begin{array}{ccc}{C}_m=C{T}_R,& D_m=D,& E=T_L{T}_R=I\end{array}\end{array}$$ For descriptor state-space models $$\begin{array}{l}\begin{array}{cc}{A}_m=T_{L}A{T}_R,& B_m=T_{L}B,\end{array}\\ \begin{array}{ccc}{C}_m=C{T}_R,& D_m=D,& E_m=T_{L}E{T}_R\end{array}\end{array}$$

For zpk and tf models, the function returns an empty value [] for this argument.

References

[1] Stewart, G. W. “A Krylov--Schur Algorithm for Large Eigenproblems.” SIAM Journal on Matrix Analysis and Applications 23, no. 3 (January 2002): 601–14. https://doi.org/10.1137/S0895479800371529.

Version History

Introduced in R2023b

R2026a: Support for Sparse Models: Compute truncated modal form

You can now use modalsep to compute a truncated modal form of sparse state-space models. This is helpful when you want to quickly obtain a truncated modal form of sparse models. For an example, see Compute Truncated Modal Decomposition of Sparse Models.

For more flexibility, first use reducespec to first obtain a reduced-order model, then apply modalsep to the result.