Description
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
with the poles of the system, the imaginary unit, and the residues. The parameter 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
with the corresponding system matrices , , , and given by
One notices that the system matrices and have complex entries. For rewriting the system with real matrices, we assume that is even, , and that all system poles are complex and ordered in complex conjugate pairs, i.e.,
Corresponding to the system poles, also the residues are written in complex conjugate pairs
Using this, the realization of the dynamical system can be written with matrices having real entries by
with , , , .
Numerical Values
We construct a system of order . The numerical values for the different variables are
- equally spaced in the interval ,
- equally spaced in the interval ,
- ,
- ,
- .
The frequency response of the transfer function is plotted for parameter values in Fig. 2.
Other interesting plots result for small values of the parameter . For example, for or , the peaks in the frequency response become more pronounced, since the poles move closer to the imaginary axis.
For , 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 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:
System dimensions:
, , ,
System variants:
Synth_matrices: , arbitrary even order 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. http://modelreduction.org/index.php/Synthetic_parametric_model
@MISC{morwiki_synth_pmodel, author = {{The MORwiki Community}}, title = {Synthetic parametric model}, howpublished = {hosted at {MORwiki} -- Model Order Reduction Wiki}, url = {http://modelreduction.org/index.php/Synthetic_parametric_model}, year = 2005 }