Anonymous
×
Create a new article
Write your page title here:
We currently have 106 articles on MOR Wiki. Type your article name above or click on one of the titles below and start writing!



Synthetic parametric model: Difference between revisions

Werner (talk | contribs)
Rewritten first sections, added planned sections.
Werner (talk | contribs)
Fully changed basic structure.
Line 12: Line 12:
<figure id="fig1">[[Image:synth_poles.png|frame|right|border|Poles of synthetic parametric system, for parameter values 1/50 (blue), 1/20 (green), 1/10 (red), 1/5 (teal), 1/2 (purple), 1 (yellow).]]</figure>
<figure id="fig1">[[Image:synth_poles.png|frame|right|border|Poles of synthetic parametric system, for parameter values 1/50 (blue), 1/20 (green), 1/10 (red), 1/5 (teal), 1/2 (purple), 1 (yellow).]]</figure>


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.
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 <xr id="fig1"/>).
Also, the decay of the Hankel singular values can be changed indirectly through the parameter.
Also, the decay of the Hankel singular values can be changed indirectly through the parameter.


Line 21: Line 21:
:<math>
:<math>
  \begin{align}
  \begin{align}
   H(s,\varepsilon) & = \sum_{k=1}^{n}\frac{r_{k}}{s-p_{k}}\\
   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})},
   & = \sum_{k=1}^{N}\frac{r_{k}}{s-(\varepsilon a_{k} + jb_{k})},
  \end{align}
  \end{align}
</math>
</math>
Line 32: Line 32:
:<math>
:<math>
  \begin{align}
  \begin{align}
   H(s,\varepsilon) = \widehat{C}(sI_{n} - (\varepsilon \widehat{A}_{\varepsilon} + \widehat{A}_{0}))^{-1}\widehat{B},
   H(s,\varepsilon) = \widehat{C}(sI_{N} - (\varepsilon \widehat{A}_{\varepsilon} + \widehat{A}_{0}))^{-1}\widehat{B},
  \end{align}
  \end{align}
</math>
</math>


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


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


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


:<math>
:<math>
Line 56: Line 56:
   p_{2} & = \varepsilon a_{1} - jb_{1},\\
   p_{2} & = \varepsilon a_{1} - jb_{1},\\
   & \ldots\\
   & \ldots\\
   p_{n-1} & = \varepsilon a_{m} + jb_{m},\\
   p_{N-1} & = \varepsilon a_{m} + jb_{m},\\
   p_{n} & = \varepsilon a_{m} - jb_{m}.
   p_{N} & = \varepsilon a_{m} - jb_{m}.
  \end{align}
  \end{align}
</math>
</math>
Line 68: Line 68:
   r_{2} & = c_{1} - jd_{1},\\
   r_{2} & = c_{1} - jd_{1},\\
   & \ldots\\
   & \ldots\\
   r_{n-1} & = c_{m} + jd_{m},\\
   r_{N-1} & = c_{m} + jd_{m},\\
   r_n & = c_{m} - jd_{m}.
   r_N & = c_{m} - jd_{m}.
\end{align}
\end{align}
</math>
</math>
Line 84: Line 84:
</math>
</math>


with <math>A_{\varepsilon, k} = \begin{bmatrix} a_{k} & 0  \\ 0 & a_{k} \end{bmatrix}</math>, <math>A_{0, k} = \begin{bmatrix} 0 & b_{k} \\ -b_{k} & 0 \end{bmatrix}</math>, <math>B_{k} = \begin{bmatrix} 2 \\ 0 \end{bmatrix}</math>, <math>C_{k} = \begin{bmatrix} c_k, & d_k \end{bmatrix}</math>.
with <math>A_{\varepsilon, k} = \begin{bmatrix} a_{k} & 0  \\ 0 & a_{k} \end{bmatrix}</math>, <math>A_{0, k} = \begin{bmatrix} 0 & b_{k} \\ -b_{k} & 0 \end{bmatrix}</math>, <math>B_{k} = \begin{bmatrix} 2 \\ 0 \end{bmatrix}</math>, <math>C_{k} = \begin{bmatrix} c_{k}, & d_{k} \end{bmatrix}</math>.
 
<figure id="fig2">[[Image:synth_freq_resp.png|frame|border|right|Frequency response of synthetic parametric system, for parameter values 1/50 (blue), 1/20 (green), 1/10 (red), 1/5 (teal), 1/2 (purple), 1 (yellow).]]</figure>


===Numerical Values===
===Numerical Values===


We construct a system of order <math>n = 100</math>. The numerical values for the different variables are
<figure id="fig3">[[Image:synth_hsv.png|frame|border|right|Hankel singular values of synthetic parametric system, for parameter values 1/50 (blue), 1/20 (green), 1/10 (red), 1/5 (teal), 1/2 (purple), 1 (yellow).]]</figure>


* <math>a_i </math> equally spaced in <math> [-10^3, -10]</math>,
We construct a system of order <math>N = 100</math>.
The numerical values for the different variables are


* <math>b_i </math> equally spaced in <math>[10, 10^3]</math>,
* <math>a_{k}</math> equally spaced in the interval <math>[-10^3, -10]</math>,


* <math> c_i = 1</math>,
* <math>b_{k}</math> equally spaced in the interval <math>[10, 10^3]</math>,


* <math> d_i = 0</math>,
* <math>c_{k} = 1</math>,


* <math>\varepsilon</math><math> \in [1/50,1]</math>.
* <math>d_{k} = 0</math>,


* <math>\varepsilon \in \left[\frac{1}{50}, 1\right]</math>.


In MATLAB, the system matrices are easily formed as follows:


<source lang="matlab">
The frequency response of the transfer function <math>H(s,\varepsilon) = C(sI_{N}-(\varepsilon A_{\varepsilon} + A_{0}))^{-1}B</math> is plotted for parameter values <math>\varepsilon \in \left[\frac{1}{50}, \frac{1}{20}, \frac{1}{10}, \frac{1}{5}, \frac{1}{2}, 1\right]</math> in <xr id="fig2"/>.
n = 100;
a = -linspace(1e1,1e3,n/2).';  b = linspace(1e1,1e3,n/2).';
c = ones(n/2,1);                d = zeros(n/2,1);
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;
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);
</source>


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


The above system matrices <math>A_\varepsilon, A_0, B, C</math> are also available in [http://math.nist.gov/MatrixMarket/formats.html MatrixMarket] format  [[Media:Synth_matrices.tar.gz|Synth_matrices.tar.gz]].
For <math>\varepsilon \in \left[\frac{1}{50}, \frac{1}{20}, \frac{1}{10}, \frac{1}{5}, \frac{1}{2}, 1\right]</math>, we also plotted the decay of the Hankel singular values in <xr id="fig3"/>.
Notice that for small values of the parameter, the decay of the Hankel singular values is very slow.


== Plots ==
==Data and Scripts==


We plot the frequency response <math>H(s,\varepsilon) = C\big(sI-\varepsilon A_\varepsilon - A_0\big)^{-1}B</math> and poles for parameter values <math>\varepsilon \in [1/50, 1/20, 1/10, 1/5, 1/2, 1] </math>.
This benchmark includes one data set. The matrices can be downloaded in the [http://math.nist.gov/MatrixMarket/formats.html MatrixMarket] format:
* [[Media:Synth_matrices.tar.gz|Synth_matrices.tar.gz]] (1.28 kB)
The matrix name is used as an extension of the matrix file.


System data of arbitrary order <math>N</math> can be generated in MATLAB or Octave by the following script:


:[[Image:synth_freq_resp.png|frame|border|left|Frequency response of synthetic parametric system, for parameter values 1/50 (blue), 1/20 (green), 1/10 (red), 1/5 (teal), 1/2 (purple), 1 (yellow).]]
<div class="thumbinner" style="width:25%;text-align:left;"><!--[[Media:ilp.m|-->
<source lang="matlab">
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);


In MATLAB, the plots are generated using the following commands:
% 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);
</source>
<!--]]--></div>
Beside that, the plots in <xr id="fig1"/> and <xr id="fig2"/> can be generated in MATLAB and Octave using the following script:
<div class="thumbinner" style="width:25%;text-align:left;"><!--[[Media:ilp.m|-->
<source lang="matlab">
<source lang="matlab">
r(1:2:n-1,1) = c+1j*d;     r(2:2:n,1) = c-1j*d;
% Get residues of the system.
ep = [1/50; 1/20; 1/10; 1/5; 1/2; 1];                       % parameter epsilon
r(1:2:n-1, 1) = c + 1j * d;
jw = 1j*linspace(0,1.2e3,5000).';                           % frequency grid
r(2:2:n, 1)   = c - 1j * d;
for j = 1:length(ep)
 
  p(:,j) = [ep(j)*a+1j*b; ep(j)*a-1j*b];                   % poles
ep = [1/50; 1/20; 1/10; 1/5; 1/2; 1]; % Parameter epsilon.
  [jww,pp] = meshgrid(jw,p(:,j));
jw = 1j * linspace(0, 1.2e3, 5000).'; % Frequency grid.
  Hjw(j,:) = (r.')*(1./(jww-pp));                           % freq. resp.
 
end
% Computations for all given parameter values.
figureloglog(imag(jw),abs(Hjw),'LineWidth',2)
p  = zeros(2 * length(a), length(ep));
          axis tight,    xlim([6 1200])
Hjw = zeros(length(ep), 5000);
          xlabel('frequency (rad/sec)')
for k = 1:length(ep)
          ylabel('magnitude')
    p(:, k)   = [ep(k) * a + 1j * b; ep(k) * a - 1j * b]; % Poles.
          title('Frequency response for different \epsilon')
    [jww, pp] = meshgrid(jw, p(:, k));
figureplot(real(p),imag(p),'.')
    Hjw(k, :) = (r.') * (1 ./ (jww - pp)); % Frequency response.
          title('Poles for different \epsilon')
end
 
% Plot poles.
figure;
loglog(imag(jw), abs(Hjw), 'LineWidth', 2);
axis tight;
xlim([6 1200]);
xlabel('frequency (rad/sec)');
ylabel('magnitude');
title('Frequency response for different \epsilon');
 
% Plot frequency response.
figure;
plot(real(p), imag(p), '.');
title('Poles for different \epsilon');
</source>
</source>
<!--]]--></div>


==Dimensions==
System structure:
:<math>
\begin{align}
\dot{x}(t) &= (\varepsilon A_{\varepsilon} + A_{0})x(t) + Bu(t), \\
y(t) &= Cx(t)
\end{align}
</math>


Other interesting plots result for small values of the parameter. For example, for <math>\varepsilon = 1/100, 1/1000 </math>, the peaks in the frequency response become more pronounced, since the poles move closer to the imaginary axis.


System dimensions:


Next, for <math>\varepsilon \in [1/50, 1/20, 1/10, 1/5, 1/2, 1] </math>, we also plot the decay of the Hankel singular values. Notice that for small values of the parameter, the decay of the Hankel singular values is very slow.
<math>A_{\varepsilon} \in \mathbb{R}^{N \times N}</math>,
<math>A_{0} \in \mathbb{R}^{N \times N}</math>,
<math>B \in \mathbb{R}^{N \times 1}</math>,
<math>C \in \mathbb{R}^{1 \times N}</math>


[[Image:synth_hsv.png|frame|border|center|Hankel singular values of synthetic parametric system, for parameter values 1/50 (blue), 1/20 (green), 1/10 (red), 1/5 (teal), 1/2 (purple), 1 (yellow).]]


==Data==
System variants:


==Dimensions==
<tt>Synth_matrices</tt>: <math>N = 100</math>


==Citation==
==Citation==
To cite this benchmark and its data:
:: Oberwolfach Benchmark Collection '''Synthetic parametric model'''. hosted at MORwiki - Model Order Reduction Wiki, 2005. http://modelreduction.org/index.php/Synthetic_parametric_model
@MISC{morwiki_synth_pmodel,
  author =      {Oberwolfach Benchmark Collection},
  title =        {Synthetic parametric model},
  howpublished = {hosted at {MORwiki} -- Model Order Reduction Wiki},
  url =          {<nowiki>http://modelreduction.org/index.php/Synthetic_parametric_model</nowiki>},
  year =        2005
}


==Contact==
==Contact==


''[[User:Ionita]]''
''[[User:Ionita]]''

Revision as of 21:49, 5 April 2018


Description

Poles of synthetic parametric system, for parameter values 1/50 (blue), 1/20 (green), 1/10 (red), 1/5 (teal), 1/2 (purple), 1 (yellow).

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 xx--CrossReference--dft--fig1--xx). 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

H(s,ε)=k=1Nrkspk=k=1Nrks(εak+jbk),

with pk=εak+jbk the poles of the system, j the imaginary unit, and rk 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

H(s,ε)=C^(sIN(εA^ε+A^0))1B^,

with the corresponding system matrices A^εN×N, A^0N×N, B^N, and C^TN given by

εA^ε+A^0=ε[a1aN]+[jb1jbN],B^=[1,,1]T,C^=[r1,,rn].

One notices that the system matrices A^0 and 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.,

p1=εa1+jb1,p2=εa1jb1,pN1=εam+jbm,pN=εamjbm.

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

r1=c1+jd1,r2=c1jd1,rN1=cm+jdm,rN=cmjdm.

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

Aε=[Aε,1Aε,m],A0=[A0,1A0,m],B=[B1Bm],C=[C1,,Cm],

with Aε,k=[ak00ak], A0,k=[0bkbk0], Bk=[20], Ck=[ck,dk].

Frequency response of synthetic parametric system, for parameter values 1/50 (blue), 1/20 (green), 1/10 (red), 1/5 (teal), 1/2 (purple), 1 (yellow).

Numerical Values

Hankel singular values of synthetic parametric system, for parameter values 1/50 (blue), 1/20 (green), 1/10 (red), 1/5 (teal), 1/2 (purple), 1 (yellow).

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

  • ak equally spaced in the interval [103,10],
  • bk equally spaced in the interval [10,103],
  • ck=1,
  • dk=0,
  • ε[150,1].


The frequency response of the transfer function H(s,ε)=C(sIN(εAε+A0))1B is plotted for parameter values ε[150,120,110,15,12,1] in xx--CrossReference--dft--fig2--xx.

Other interesting plots result for small values of the parameter ε. For example, for ε=1100 or 11000, the peaks in the frequency response become more pronounced, since the poles move closer to the imaginary axis.

For ε[150,120,110,15,12,1], we also plotted the decay of the Hankel singular values in xx--CrossReference--dft--fig3--xx. 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:

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

System data of arbitrary 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);

Beside that, the plots in xx--CrossReference--dft--fig1--xx and xx--CrossReference--dft--fig2--xx 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;
loglog(imag(jw), abs(Hjw), 'LineWidth', 2);
axis tight;
xlim([6 1200]);
xlabel('frequency (rad/sec)');
ylabel('magnitude');
title('Frequency response for different \epsilon');

% Plot frequency response.
figure;
plot(real(p), imag(p), '.');
title('Poles for different \epsilon');

Dimensions

System structure:

x˙(t)=(εAε+A0)x(t)+Bu(t),y(t)=Cx(t)


System dimensions:

AεN×N, A0N×N, BN×1, C1×N


System variants:

Synth_matrices: N=100

Citation

To cite this benchmark and its data:

Oberwolfach Benchmark Collection Synthetic parametric model. hosted at MORwiki - Model Order Reduction Wiki, 2005. http://modelreduction.org/index.php/Synthetic_parametric_model
@MISC{morwiki_synth_pmodel,
  author =       {Oberwolfach Benchmark Collection},
  title =        {Synthetic parametric model},
  howpublished = {hosted at {MORwiki} -- Model Order Reduction Wiki},
  url =          {http://modelreduction.org/index.php/Synthetic_parametric_model},
  year =         2005
}

Contact

User:Ionita