modalsep

Compute modal decomposition

Since R2023b

Syntax

[H,H0] = modalsep(G)

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

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

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

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.

[H,H0] = modalsep(G,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.

[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,info] = modalsep(G,___) 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.

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 \mathit{z}<0.033$

• Region 2: $0.033\le \mathit{z}<0.040$

• Region 3: $\mathit{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

Input Arguments

G — Linear time-invariant model

ss | tf | zpk

Linear time-invariant model, specified as an ss, tf, or zpk model.

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])

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.

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.

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 modal components.
DCGain	Normalized DC contribution of modal components, returned as a vector of length nc-by-1, where nc is the number of modal components.
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. 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 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, 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.

Version History

Introduced in R2023b