Synthetic parametric model

Description

Figure 1: System poles for different parameter values.

On this page you will find a synthetic parametric model with one parameter for which one can easily experiment with different system orders, values of the parameter, as well as different poles and residues (see Fig. 1). Also, the decay of the Hankel singular values can be changed indirectly through the parameter.

Model

We consider a dynamical system in the frequency domain given by its pole-residue form of the transfer function

$\begin{align} H(s,\varepsilon) & = \sum_{k=1}^{N}\frac{r_{k}}{s-p_{k}}\\ & = \sum_{k=1}^{N}\frac{r_{k}}{s-(\varepsilon a_{k} + jb_{k})}, \end{align}$

with $p_{k} = \varepsilon a_{k} + jb_{k}$ the poles of the system, $j$ the imaginary unit, and $r_{k}$ the residues. The parameter $\varepsilon$ is used to scale the real part of the system poles. We can write down the state-space realization of the system's transfer function as

$\begin{align} H(s,\varepsilon) = \widehat{C}(sI_{N} - (\varepsilon \widehat{A}_{\varepsilon} + \widehat{A}_{0}))^{-1}\widehat{B}, \end{align}$

with the corresponding system matrices $\widehat{A}_{\varepsilon} \in \mathbb{R}^{N \times N}$, $\widehat{A}_{0} \in \mathbb{C}^{N \times N}$, $\widehat{B} \in \mathbb{R}^{N}$, and $\widehat{C}^{T} \in \mathbb{C}^{N}$ given by

$\begin{align} \varepsilon\widehat{A}_{\varepsilon} + \widehat{A}_{0} & = \varepsilon \begin{bmatrix} a_{1} & & \\ & \ddots & \\ & & a_{N} \end{bmatrix} + \begin{bmatrix} jb_{1} & & \\ & \ddots & \\ & & jb_{N} \end{bmatrix},\\ \widehat{B} & = \begin{bmatrix}1, & \ldots, & 1 \end{bmatrix}^{T},\\ \widehat{C} & = \begin{bmatrix}r_{1}, & \ldots, & r_{n} \end{bmatrix}. \end{align}$

One notices that the system matrices $\widehat{A}_{0}$ and $\widehat{C}$ have complex entries. For rewriting the system with real matrices, we assume that $N$ is even, $N=2m$, and that all system poles are complex and ordered in complex conjugate pairs, i.e.,

$\begin{align} p_{1} & = \varepsilon a_{1} + jb_{1},\\ p_{2} & = \varepsilon a_{1} - jb_{1},\\ & \ldots\\ p_{N-1} & = \varepsilon a_{m} + jb_{m},\\ p_{N} & = \varepsilon a_{m} - jb_{m}. \end{align}$

Corresponding to the system poles, also the residues are written in complex conjugate pairs

$\begin{align} r_{1} & = c_{1} + jd_{1},\\ r_{2} & = c_{1} - jd_{1},\\ & \ldots\\ r_{N-1} & = c_{m} + jd_{m},\\ r_N & = c_{m} - jd_{m}. \end{align}$

Using this, the realization of the dynamical system can be written with matrices having real entries by

$\begin{align} A_{\varepsilon} & = \begin{bmatrix} A_{\varepsilon, 1} & & \\ & \ddots & \\ & & A_{\varepsilon, m} \end{bmatrix}, & A_{0} & = \begin{bmatrix} A_{0, 1} & & \\ & \ddots & \\ & & A_{0, m} \end{bmatrix}, & B & = \begin{bmatrix} B_{1} \\ \vdots \\ B_{m} \end{bmatrix}, & C & = \begin{bmatrix} C_{1}, & \cdots, & C_{m} \end{bmatrix}, \end{align}$

with $A_{\varepsilon, k} = \begin{bmatrix} a_{k} & 0 \\ 0 & a_{k} \end{bmatrix}$, $A_{0, k} = \begin{bmatrix} 0 & b_{k} \\ -b_{k} & 0 \end{bmatrix}$, $B_{k} = \begin{bmatrix} 2 \\ 0 \end{bmatrix}$, $C_{k} = \begin{bmatrix} c_{k}, & d_{k} \end{bmatrix}$.

Figure 2: Frequency response of synthetic parametric system for different parameter values.

Numerical Values

Figure 3: Hankel singular values of synthetic parametric system for different parameter values.

We construct a system of order $N = 100$. The numerical values for the different variables are

• $a_{k}$ equally spaced in the interval $[-10^3, -10]$,

• $b_{k}$ equally spaced in the interval $[10, 10^3]$,

• $c_{k} = 1$,

• $d_{k} = 0$,

• $\varepsilon \in \left[\frac{1}{50}, 1\right]$.

The frequency response of the transfer function $H(s,\varepsilon) = C(sI_{N}-(\varepsilon A_{\varepsilon} + A_{0}))^{-1}B$ is plotted for parameter values $\varepsilon \in \left[\frac{1}{50}, \frac{1}{20}, \frac{1}{10}, \frac{1}{5}, \frac{1}{2}, 1\right]$ in Fig. 2.

Other interesting plots result for small values of the parameter $\varepsilon$. For example, for $\varepsilon = \frac{1}{100}$ or $\frac{1}{1000}$, the peaks in the frequency response become more pronounced, since the poles move closer to the imaginary axis.

For $\varepsilon \in \left[\frac{1}{50}, \frac{1}{20}, \frac{1}{10}, \frac{1}{5}, \frac{1}{2}, 1\right]$, we also plotted the decay of the Hankel singular values in Fig. 3. Notice that for small values of the parameter, the decay of the Hankel singular values is very slow.

Data and Scripts

This benchmark includes one data set. The matrices can be downloaded in the MatrixMarket format:

• Synth_matrices.tar.gz (1.28 kB)

The matrix name is used as an extension of the matrix file.

System data of arbitrary even order $N$ can be generated in MATLAB or Octave by the following script:

N = 100; % Order of the resulting system.

% Set coefficients.
a = -linspace(1e1, 1e3, N/2).';
b =  linspace(1e1, 1e3, N/2).';
c = ones(N/2, 1);
d = zeros(N/2, 1);

% Build 2x2 submatrices.
aa(1:2:N-1, 1) = a;
aa(2:2:N, 1)   = a;
bb(1:2:N-1, 1) = b;
bb(2:2:N-2, 1) = 0;

% Set up system matrices.
Ae = spdiags(aa, 0, N, N);
A0 = spdiags([0; bb], 1, N, N) + spdiags(-bb, -1, N, N);
B  = 2 * sparse(mod(1:N, 2)).';
C(1:2:N-1) = c.';
C(2:2:N)   = d.';
C          = sparse(C);

or in Python:

import numpy as np
import scipy.sparse as sps

N = 100  # Order of the resulting system.

# Set coefficients.
a = -np.linspace(1e1, 1e3, N//2)
b = np.linspace(1e1, 1e3, N//2)
c = np.ones(N//2)
d = np.zeros(N//2)

# Build 2x2 submatrices.
aa = np.empty(N)
aa[::2] = a
aa[1::2] = a
bb = np.zeros(N)
bb[::2] = b

# Set up system matrices.
Ae = sps.diags(aa, format='csc')
A0 = sps.diags([bb, -bb], [1, -1], (N, N), format='csc')
B = np.zeros((N, 1))
B[::2, :] = 2
C = np.empty((1, N))
C[0, ::2] = c
C[0, 1::2] = d

Beside that, the plots in Fig. 1 and Fig. 2 can be generated in MATLAB and Octave using the following script:

% Get residues of the system.
r(1:2:N-1, 1) = c + 1j * d;
r(2:2:N, 1)   = c - 1j * d;

ep = [1/50; 1/20; 1/10; 1/5; 1/2; 1]; % Parameter epsilon.
jw = 1j * linspace(0, 1.2e3, 5000).'; % Frequency grid.

% Computations for all given parameter values.
p   = zeros(2 * length(a), length(ep));
Hjw = zeros(length(ep), 5000);
for k = 1:length(ep)
    p(:, k)   = [ep(k) * a + 1j * b; ep(k) * a - 1j * b]; % Poles.
    [jww, pp] = meshgrid(jw, p(:, k));
    Hjw(k, :) = (r.') * (1 ./ (jww - pp)); % Frequency response.
end

% Plot poles.
figure;
plot(real(p), imag(p), '.', 'MarkerSize', 12);
xlabel('Re(p)');
ylabel('Im(p)');
legend( ...
    '\epsilon = 1/50', ...
    '\epsilon = 1/20', ...
    '\epsilon = 1/10', ...
    '\epsilon = 1/5', ...
    '\epsilon = 1/2', ...
    '\epsilon = 1');

% Plot frequency response.
figure;
loglog(imag(jw), abs(Hjw), 'LineWidth', 2);
axis tight;
xlim([6 1200]);
xlabel('frequency (rad/sec)');
ylabel('magnitude');
legend( ...
    '\epsilon = 1/50', ...
    '\epsilon = 1/20', ...
    '\epsilon = 1/10', ...
    '\epsilon = 1/5', ...
    '\epsilon = 1/2', ...
    '\epsilon = 1');

or in Python:

import matplotlib.pyplot as plt

# Get residues of the system.
r = np.empty(N, dtype=complex)
r[::2] = c + 1j * d
r[1::2] = c - 1j * d

ep = [1/50, 1/20, 1/10, 1/5, 1/2, 1]  # Parameter epsilon.
jw = 1j * np.geomspace(6, 1.2e3, 5000)  # Frequency grid.

# Computations for all given parameter values.
p = np.zeros((len(ep), N), dtype=complex)
Hjw = np.zeros((len(ep), len(jw)), dtype=complex)
for k, epk in enumerate(ep):
    # Poles.
    p[k, :N//2] = epk * a + 1j * b
    p[k, N//2:] = epk * a - 1j * b
    # Frequency response.
    Hjw[k, :] = (r / (jw[:, np.newaxis] - p[k])).sum(axis=1)

# Plot poles.
fig, ax = plt.subplots()
for k, epk in enumerate(ep):
    ax.plot(p[k].real, p[k].imag, '.', label=fr'$\varepsilon$ = {epk}')
ax.autoscale(tight=True)
ax.set_xlabel('Re(p)')
ax.set_ylabel('Im(p)')
ax.legend()

# Plot frequency response.
fig, ax = plt.subplots()
for k, epk in enumerate(ep):
    ax.loglog(jw.imag, np.abs(Hjw[k]), label=fr'$\varepsilon$ = {epk}', linewidth=2)
ax.autoscale(tight=True)
ax.set_xlabel('frequency (rad/sec)')
ax.set_ylabel('magnitude')
ax.legend()

Dimensions

System structure:

$\begin{align} \dot{x}(t) &= (\varepsilon A_{\varepsilon} + A_{0})x(t) + Bu(t), \\ y(t) &= Cx(t) \end{align}$

System dimensions:

$A_{\varepsilon} \in \mathbb{R}^{N \times N}$, $A_{0} \in \mathbb{R}^{N \times N}$, $B \in \mathbb{R}^{N \times 1}$, $C \in \mathbb{R}^{1 \times N}$

System variants:

Synth_matrices: $N = 100$, arbitrary even order $N$ by using the script

Citation

To cite this benchmark and its data:

The MORwiki Community, Synthetic parametric model. hosted at MORwiki - Model Order Reduction Wiki, 2005. https://modelreduction.org/morwiki/Synthetic_parametric_model

@MISC{morwiki_synth_pmodel,
  author =       {{The MORwiki Community}},
  title =        {Synthetic parametric model},
  howpublished = {hosted at {MORwiki} -- Model Order Reduction Wiki},
  url =          {https://modelreduction.org/morwiki/Synthetic_parametric_model},
  year =         2005
}

Contact

Antonio Cosmin Ionita