isstable

Determine if dynamic system model is stable

Syntax

B = isstable(sys)

B = isstable(sys,'elem')

B = isstable(___,Name=Value)

Description

B = isstable(sys) returns a logical value of 1 (true) if the dynamic system model sys has stable dynamics, and a logical value of 0 (false) otherwise. If sys is a model array, then the function returns 1 only if all the models in sys are stable.

isstable returns a logical value of 1 (true) for stability of a dynamic system if:

• In continuous-time systems, all the poles lie in the open left half of the complex plane.

• In discrete-time systems, all the poles lie inside the open unit disk.

isstable is supported only for analytical models with a finite number of poles.

B = isstable(sys,'elem') returns a logical array of the same dimensions as the model array sys. The logical array indicates which models in sys are stable.

B = isstable(___,Name=Value) determines stability of a sparse model sys by computing a subset of poles based on options specified using one or more specified name-value arguments. The function ignores name-value arguments when sys is a nonsparse model. (since R2025a)

Examples

Determine Stability of Discrete-Time Transfer Function Model

Determine the stability of this discrete-time SISO transfer function model with a sample time of 0.1 seconds.

$$\mathrm{sys}\left(\mathit{z}\right)=\frac{2\mathit{z}}{4{\mathit{z}}^3+3\mathit{z}-1}$$

Create the discrete-time transfer function model.

sys = tf([2,0],[4,0,3,-1],0.1);

Examine the poles of the system.

P = abs(pole(sys))

P = 3×1

    0.9159
    0.9159
    0.2980

All the poles of the transfer function model have a magnitude less than 1, so all the poles lie within the open unit disk and the system is stable.

Confirm the stability of the model using isstable.

B = isstable(sys)

B = logical
   1

The system sys is stable.

Determine Stability of Continuous-Time Zero-Pole-Gain Model

Determine the stability of this continuous-time zero-pole-gain model.

$$\mathrm{sys}\left(\mathit{s}\right)=\frac{2}{\left(\mathit{s}+2+3\mathit{j}\right)\left(\mathit{s}+2-3\mathit{j}\right)\left(\mathit{s}-0.5\right)}$$

Create the model as a zpk model object by specifying the zeros, poles, and gain.

sys = zpk([],[-2-3*j,-2+3*j,0.5],2);

Because one pole of the model lies in the right half of the complex plane, the system is unstable.

Confirm the instability of the model using isstable.

B = isstable(sys)

B = logical
   0

The system sys is unstable.

Determine Stability of Models in Model Array

Determine the stability of an array of SISO transfer function models with poles varying from -2 to 2.

$$\left[\frac{1}{\mathit{s}+2},\frac{1}{\mathit{s}+1},\frac{1}{\mathit{s}},\frac{1}{\mathit{s}-1},\frac{1}{\mathit{s}-\mathrm{2}}\right]$$

To create the array, first initialize an array of dimension [length(a),1] with zero-valued SISO transfer functions.

a = [-2:2];
sys = tf(zeros(1,1,length(a)));

Populate the array with transfer functions of the form 1/(s-a).

for j = 1:length(a)
    sys(1,1,j) = tf(1,[1 -a(j)]);
end

isstable can tell you whether all the models in model array are stable or each individual model is stable.

Examine the stability of the model array.

B_all = isstable(sys)

B_all = logical
   0

By default, isstable returns a single logical value that is 1 (true) only if all models in the array are stable. sys contains some models with nonnegative poles, which are not stable. Therefore, isstable returns 0 (false) for the entire array.

Examine the stability of each model in the array by using 'elem' flag.

B_elem = isstable(sys,'elem')

B_elem = 5×1 logical array

   1
   1
   0
   0
   0

The function returns an array of logical values that indicate the stability of the corresponding entry in the model array. For example, B_elem(2) is 1, which indicates that the second model in the array, sys(1,1,2) is stable. This is because sys(1,1,2) has a pole at -1.

Determine Stability of Sparse Model

Since R2025a

For this example, consider tuningForkData.mat that contains the sparse second-order model of a tuning fork. For sparse state-space models, isstable determines stability based on a subset of compute in a frequency range of interest.

Load the sparse matrices from tuningForkData.mat into the workspace and create the mechss model object.

load tuningForkData.mat;
sys = mechss(M,[],K,B,F)

Sparse continuous-time second-order model with 2 outputs, 1 inputs, and 3024 degrees of freedom.
Model Properties

Use "spy" and "showStateInfo" to inspect model structure. 
Type "help mechssOptions" for available solver options for this model.

Check the stability in the frequency range up to 10000 rad/s.

b = isstable(sys,Focus=[0 1e4])

b = logical
   0

The tuning fork model sys is an undamped model and isstable returns 0 as expected. You can verify this using the damp function.

damp(sys,Focus=[0 1e4])

                                                                       
         Pole              Damping       Frequency      Time Constant  
                                       (rad/seconds)      (seconds)    
                                                                       
  0.00e+00 + 2.96e+03i     0.00e+00       2.96e+03           Inf       
  0.00e+00 - 2.96e+03i     0.00e+00       2.96e+03           Inf       
  0.00e+00 + 2.96e+03i     0.00e+00       2.96e+03           Inf       
  0.00e+00 - 2.96e+03i     0.00e+00       2.96e+03           Inf       
  0.00e+00 + 3.01e+03i     0.00e+00       3.01e+03           Inf       
  0.00e+00 - 3.01e+03i     0.00e+00       3.01e+03           Inf       
  0.00e+00 + 3.02e+03i     0.00e+00       3.02e+03           Inf       
  0.00e+00 - 3.02e+03i     0.00e+00       3.02e+03           Inf       

Now add Rayleigh damping $\mathit{C}=\alpha \mathit{M}+\beta \mathit{K}$ with minimum damping 0.01 at frequency 3000 rad/s.

wmin = 3000;
zmin = 0.01;
alpha = zmin*wmin;
beta = zmin/wmin;
C_d = alpha*M+beta*K;
sys_d = mechss(M,C_d,K,B,F)

Sparse continuous-time second-order model with 2 outputs, 1 inputs, and 3024 degrees of freedom.
Model Properties

Use "spy" and "showStateInfo" to inspect model structure. 
Type "help mechssOptions" for available solver options for this model.

Check the stability for the damped model.

isstable(sys_d,Focus=[0 1e4])

ans = logical
   1

With added damping, all the computed poles of the model lie in the left-half plane.

damp(sys_d,Focus=[0 1e4])

                                                                       
         Pole              Damping       Frequency      Time Constant  
                                       (rad/seconds)      (seconds)    
                                                                       
 -2.96e+01 + 2.96e+03i     1.00e-02       2.96e+03         3.38e-02    
 -2.96e+01 - 2.96e+03i     1.00e-02       2.96e+03         3.38e-02    
 -2.96e+01 + 2.96e+03i     1.00e-02       2.96e+03         3.38e-02    
 -2.96e+01 - 2.96e+03i     1.00e-02       2.96e+03         3.38e-02    
 -3.01e+01 + 3.01e+03i     1.00e-02       3.01e+03         3.33e-02    
 -3.01e+01 - 3.01e+03i     1.00e-02       3.01e+03         3.33e-02    
 -3.02e+01 + 3.02e+03i     1.00e-02       3.02e+03         3.31e-02    
 -3.02e+01 - 3.02e+03i     1.00e-02       3.02e+03         3.31e-02    

Input Arguments

sys — Dynamic system

dynamic system model | model array

Dynamic system, specified as a SISO or MIMO dynamic system model or an array of SISO or MIMO dynamic system models. Dynamic systems that you can use include continuous-time or discrete-time numeric LTI models such as tf, zpk, or ss models.

If sys is a generalized state-space model genss or an uncertain state-space model uss (Robust Control Toolbox), isstable checks the stability of the current or nominal value of sys.

If sys is a sparse state-space model (sparss or mechss), the software computes a subset of poles in a specified frequency band of focus to check the stability. For sparse models, use the name-value arguments to specify computation options. If you do not specify any options, the software computes the first 1000 poles with smallest magnitude. For mechss models with Rayleigh damping, the software computes poles from the eigenvalues of (K,M). Otherwise, the software computes the poles from the equivalent sparss model. (since R2025a)

If sys is an array of models, isstable checks the stability of every model in the array.

• If you use B = isstable(sys), the output is 1 (true) only if all the models in the array are stable.

• If you use B = isstable(sys,'elem'), the output is a logical array, the entries of which indicate the stability of the corresponding entry in the model array.

For more information on model arrays, see Model Arrays.

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: B = isstable(sparseSys,Focus=[0 10],Display="off")

Focus — Frequency range of interest

[0 Inf] (default) | vector

Since R2025a

Frequency range of interest, specified as a vector of form [fmin,fmax]. When you specify a frequency range of focus, the algorithm computes only the poles with natural frequency in this range. For discrete-time models, the software approximates the equivalent natural frequency through Tustin transform.

For symmetric semi-definite sparse models, you can specify any frequency range with 0 ≤ fmin < fmax ≤ ∞. For nonsymmetric sparss models, you must specify frequency range as either [0,fmax] or [fmin,inf].

Since the software computes all poles in the specified frequency range, you typically specify a low-frequency range to limit computing a large number of poles. By default, the focus is unspecified ([0 Inf]) and the algorithm computes up to MaxNumber poles.

MaxNumber — Maximum number of poles and zeros to compute

1000 (default) | positive integer

Since R2025a

Maximum number of poles and zeros to compute, specified as a positive integer. This value limits the number of poles computed by the algorithm.

Shift — Spectral shift

0 (default) | finite scalar

Since R2025a

Spectral shift, specified as a finite scalar.

The software computes poles with the natural frequency in the specified range [0,fmax] using inverse power iterations for A-sigma*E, which obtains eigenvalues closest to shift sigma. When A is singular and sigma is zero, the algorithm fails as no inverse exists. Therefore, for sparse models with integral action (s = 0 or at z = 1 for discrete-time models), you can use this option to implicitly shift poles or zeros to the value closest to this shift value. Specify a shift value that is not equal to an existing pole value of the original model.

Tolerance — Accuracy of computed poles

1e-12 (default) | positive finite scalar

Since R2025a

Tolerance for accuracy of computed poles, specified as a positive finite scalar. This value controls the convergence of computed eigenvalues in inverse power iterations.

Display — Show or hide progress report

"on" (default) | "off"

Since R2025a

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

Output Arguments

B — True or false result

1 | 0 | logical array

True or false result, returned as 1 for a stable model or 0 for an unstable model.

The 'elem' flag causes isstable to return an array of logical values with same dimensions as the model array. The values in the array indicate the stability of the corresponding entry in the model array.

Algorithms

For sparse-state space models, the software uses the Krylov--Schur algorithm [1] for inverse power iterations to compute poles in the specified frequency band.

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 R2012a

R2025a: Sparse Model Support: Determine stability based on subset of computed poles

You can now determine stability of sparse models based on a subset of computed poles in a frequency range of interest. By default, the function computes first 1000 poles with the smallest magnitude. To avoid computing a large number of poles, you can use the new syntax B = isstable(sparseSys,Name=Value) to specify options such as frequency range of interest or maximum number of poles to compute.