FrequencyResponseFitting

Specification for fitting low-order model to frequency response

Since R2025a

Description

The FrequencyResponseFitting object uses frequency response fitting to obtain low-order approximations of sparse LTI models in the frequency band of interest. The frequency response fitting method is applicable to all types of sparse models with fewer limitations than Balanced Truncation or POD. In particular, it is effective at reducing models with unstable or undamped poles. Frequency response fitting is also applicable to non-sparse models but is usually less effective than Balanced Truncation for such models.

Additionally, this method uses the AAA interpolation algorithm rather than a least-squared fitting algorithm. Accordingly, you must not use this method to fit models to noisy data as the algorithm will tend to interpolate noise and produce poor results.

Creation

The reducespec function creates a FrequencyResponseFitting model order reduction object when you use this syntax.

R = reducespec(sys,"frfit")

Here, sys is a ordinary (nonsparse), sparse, or frequency-response data (frd) LTI model. The workflow uses this object to set up MOR tasks and store results. For the full workflow, see Task-Based Model Order Reduction Workflow.

Properties

FRData — Frequency response data used for fitting

Read-only: [] (default) | frd object

This property is read-only.

Frequency response data of the original model, returned as an frd object

Fit — Rational fit in barycentric form

Read-only: [] (default) | Barycentric object

This property is read-only.

Rational fit in barycentric form, returned as a Barycentric object.

The rational fit in the barycentric form is defined as:

$$R(s)=\frac{{\displaystyle \sum \frac{F_i{W}_i}{s-z_i}}}{{\displaystyle \sum \frac{W_i}{s-z_i}}}$$

R.Fit is a Barycentric object with properties:

• z — Vector of interpolation points z_i (complex values).

• F — 3D array of frequency response values F(z_i).

• W — 3D array of weight values W_i.

Additionally, the rational fit is real valued when the original model has real coefficients and complex valued when the original model has complex coefficients.

FitError — Worst normalized error at each iteration

Read-only: [] (default) | positive scalar

This property is read-only.

Worst normalized error at each iteration, returned as a positive scalar. The fit error is the maximum value of the normalized error E_k over the data points (z_k, F_k), where E_k is:

$$E_k=\frac{{\Vert F_k-Fit(z_k)\Vert }_F}{(AbsTol+RelTol\times {\Vert F_k\Vert }_F)}$$

The algorithm achieves the desired accuracy when this worst-case ratio drops below 1.

Options — Options for frequency response fitting

FrequencyResponseFittingOptions object

Options for frequency response fitting, specified as a FrequencyResponseFittingOptions object. Use dot notation to configure options for R. For example R.Options.AbsTol = 0.05.

For more information about available options, see FrequencyResponseFittingOptions.

Object Functions

process	Run model order reduction algorithm
view (frfit)	Plot relative fit error between original and fitted model
getrom (frfit)	Obtain reduced-order models when using frequency-response fitting method

Examples

Reduce Sparse Model Order Using Frequency Response Fitting Method

This example shows how to reduce model order of a sparss model using the frequency response fitting method. The sparse state-space model used in this example is obtained from linearizing a thermal model of heat distribution in a circular cylindrical rod.

Load the model data.

load cylindricalRod.mat
sys = sparss(A,B,C,D,E);
size(sys)

Sparse state-space model with 3 outputs, 1 inputs, and 7522 states.

The thermal model contains 7522 states.

To reduce the model order, first create a frequency response fitting specification object for sys.

R = reducespec(sys,'frfit');

To customize the specification object, you can specify additional options such as a frequency grid and the feedthrough values. Since the feedthrough values for the original model is zero, you can set it to zero in the options.

w = logspace(-7,-3,20);
R.Options.Frequency = w;
R.Options.Feedthrough = zeros(3,1);
R.Options.RelTol = 1e-3;
R.Options.Display = 'off';

Run the model reduction algorithm and visualize the frequency response of the original model and the fit error.

R = process(R);
view(R)

Obtain the reduced model from the rational fit. You can specify additional arguments of getrom such that the software returns a model in descriptor form and does not enforce any stability on the fitted model.

rsys = getrom(R,Explicit=false,Stable=false);
order(rsys)

ans = 
11

Compare the relative error between two models.

sigmaplot(sys,sys-rsys,"r--",w)

The reduced model provides a good approximation of the original sparse model.

Fit Low-Order Model to Frequency Response of Tuning Fork Sparse Model

This example shows how to fit a low-order model using the frequency response fitting method. In this example, you obtain a reduced-order model of a tuning fork model. The model is a sparse second-order (mechss) model with two outputs, one input, and 3024 degrees of freedom.

Load the model.

load tuningForkData.mat
sys = mechss(M,[],K,B,F);
w = logspace(2,6,250);
size(sys)

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

Create a model order reduction object for frequency response fitting using reducespec. Specify a frequency grid and the feedthrough value for interpolation.

R = reducespec(sys,'frfit');
R.Options.Frequency = w;
R.Options.Feedthrough = [0;0];
R.Options.Display = "off";
R = process(R);

rsys1 = getrom(R,Explicit=false,Stable=true,Minimal=true);
o1 = order(rsys1)

o1 = 
90

sigmaplot(sys,rsys1,"--",w)

This returns a model with an order which is still quite high. You can further reduce the model order by focusing on just the lower frequencies. Increasing AbsTol allows for a larger error at higher frequencies where the gain is low but the mode density increases.

R.Options.AbsTol = 1e-11;
R.Options.RelTol = 1e-3;
rsys2 = getrom(R,Explicit=false,Stable=true,Minimal=true);
o2 = order(rsys2)

o2 = 
17

The algorithm now produces a reduced model with much lower order. Compare the relative error between the original and reduced models.

sigmaplot(sys,sys-rsys1,"--",sys-rsys2,"-.",w)
legend("Original",sprintf("Order %d",o1),...
    sprintf("Order %d",o2));

Reduce MIMO Model Using Frequency Response Fitting

This example shows how to reduce a model order of a multi-input multi-output (MIMO) model using the frequency response fitting method. Typically, when you reduce MIMO models using this method, the default reduced model may not have a minimal order. This example shows how to ensure the fit is minimal.

For this example, generate a random state-space model with 50 states, two outputs, and three inputs.

rng(0);
sys = rss(50,2,3);

When using the frequency response fitting method, the software requires you to specify a grid of frequencies for evaluating the frequency response. To select a suitable grid, analyze the frequency response of the model.

sigmaplot(sys)

For this model, you can specify frequencies between 0.1 rad/s and 1000 rad/s.

w = logspace(-1,3,100);

Create a model reduction object.

R = reducespec(sys,"frfit");

Specify the desired frequency grid and run the model reduction algorithm.

R.Options.Frequency = w;
R = process(R);

Obtain the reduced model from the rational fit.

rsys = getrom(R);
order(rsys)

ans = 
30

This produces a 30th order model. By default, getrom returns a stable and explicit realization of the reduced model. However, when you use this method with MIMO models, the software may produce a nonminimal order fit. To ensure a minimal order fit, you can set the Minimal argument of getrom to true. This performs an additional balanced truncation on the initially reduced model. The software further reduces this model as long as it maintains the desired fit accuracy.

Obtain the model.

rsys2 = getrom(R,Minimal=true);
order(rsys2)

ans = 
26

Now getrom produces a model with order 26. Compare the relative error between the original and reduced models.

sigmaplot(sys,sys-rsys,":",sys-rsys2,"--")
legend("sys","sys-rsys","sys-rsys2")

The fitted models provide an accurate match for the original model.

Enforcing Stability in Fitted Reduced Models

This example shows how to enforce stability in reduced models obtained using frequency response fitting method and the effects of refining the frequency grid to avoid unstable modes.

Generate a random state-space model.

rng(1) 
sys = rss(50);  
bodeplot(sys)

First, reduce the model using a frequency grid of 30 points between 0.01 rad/s and 100 rad/s.

w = logspace(-2,2,30);
R = reducespec(sys,'frfit');
R.Options.Frequency = w;

When you use getrom with the Stable argument set to false, the software runs the standard AAA algorithm which does not enforce stability of the fitted model.

rsys1 = getrom(R,Stable=false);

The low relative error shows that the reduced model is a good fit on the specified grid.

sigma(sys,sys-rsys1,w)
legend("sys","sys-rsys1",Location="best")

However, outside the grid points, the reduced model does not provide a good fit.

bode(sys,rsys1)
legend

Additionally, this fit is also unstable.

max(real(pole(rsys1)))

ans = 
0.0335

However, when Stable is set to true, the software enforces stability after the AAA algorithm has run.

rsys2 = getrom(R,Stable=true);
max(real(pole(rsys2)))

ans = 
-2.9364e-08

In this case, the stabilization works well with only a minor degradation in fit quality. This is not always the case however, post-AAA stabilization may introduce significant errors past the frequency of the first stabilized mode.

sigma(sys,sys-rsys1,sys-rsys2,"--")
legend("sys","sys-rsys1","sys-rsys2",Location="best")

In principle, stability is encoded in the frequency response. In general, unstable poles creep in where the frequency grid is too sparse to accurately capture the system dynamics. Refining the grid is often an effective, error-free way to restore stability

w = logspace(-2,2,60);
R = reducespec(sys,'frfit');
R.Options.Frequency = w;
rsys3 = getrom(R,Stable=false);
max(real(pole(rsys3)))

ans = 
3.3737e-08

This last fit rsys3 is stable and more accurate than the stabilized fit rsys2.

sigma(sys,sys-rsys2,sys-rsys3,"--")
legend("sys","sys-rsys2","sys-rsys3",Location="best")

There is a small increase in order for the refined grid as there is more data to fit.

[order(rsys2) order(rsys3)]

ans = 1×2

    17    19

Therefore, when a stable model is required and getrom returns an unstable fit, you can either enforce stability with Stable=true, or refine the grid.

Tips

• The fit can be poor between samples and outside the frequency grid used for fitting, so it is recommended to validate the fit by comparing the fitted model with the original model on a larger or denser grid. You can also provide a few points outside the range of interest to sketch the asymptotes. It is also recommended to specify the DC and feedthrough values to get the desired behavior toward w = 0 and w = Inf.

Algorithms

The software performs frequency response fitting as follows:

1. Computes the frequency response of the model at specified frequency grid.

2. Uses the AAA algorithm with the barycentric form ([1], [2]) to fit the computed data to a rational approximation

   $$R(s)=\frac{f_{\infty }{w}_{\infty }+{\displaystyle {\sum }_{j=1}^m\frac{f_j{w}_j}{s-z_j}}}{w_{\infty }+{\displaystyle {\sum }_{j=1}^m\frac{w_j}{s-z_j}}}=\frac{N(s)}{D(s)}$$

   Here, z_j are the interpolation points, f_j are the corresponding values for f(z_j), w_j are the scalar-valued weights, and f_∞ is the feedthrough value (if specified).

   The barycentric form for r-by-s MIMO models with matrix-weights is given by:

   $$R(s)=\left(F_{\infty }{W}_{\infty }+{\displaystyle \sum _{j=1}^m\frac{F_j{W}_j}{s-z_j}}\right){\left(W_{\infty }+{\displaystyle \sum _{j=1}^m\frac{W_j}{s-z_j}}\right)}^{-1}=N(s)D{(s)}^{-1}$$

   Here, W_j are square matrices of size s. This software uses this form when R.Options.Algorithm is set to "blockAAA" when fitting a MIMO model.

3. Obtains a state-space realization from the rational approximation.

References

[1] Nakatsukasa, Yuji, Olivier Sète, and Lloyd N. Trefethen. “The AAA Algorithm for Rational Approximation.” SIAM Journal on Scientific Computing 40, no. 3 (January 2018): A1494–1522. https://doi.org/10.1137/16M1106122.

[2] Gosea, Ion Victor, and Stefan Güttel. “Algorithms for the Rational Approximation of Matrix-Valued Functions.” SIAM Journal on Scientific Computing 43, no. 5 (January 2021): A3033–54. https://doi.org/10.1137/20M1324727.

Version History

Introduced in R2025a