View source for Moment-matching PMOR method

Parametric model order reduction (PMOR) methods are designed for model order reduction of parametrized
systems, where the parameters of the system play an important role
in practical applications such as Integrated Circuit (IC) design,
MEMS design, Chemical engineering etc.. The parameters could be the variables describing
geometrical measurement, material property, damping of the
systems or component flow-rate. The reduced models are constructed such that all the
parameters can be preserved with acceptable accuracy. 
Usually the time of simulating the reduced models is much shorter
than directly simulating the original large system. However, the
time of constructing the reduced model increases with the dimension
of the original system. If the original system is very large, the
process of obtaining the reduced model could become extremely slow.
The recycling algorithm considered in this paper tries to accelerate
the above process and reduce the time of deriving the reduced model
to a reasonable range.

The method introduced here is from [1][2], and applies to a linear parametrized system,
which has the following form in the frequency domain:

<math>
(E_0+s_1E_1+s_2E_2+\ldots +s_pE_p)x=Bu(s_p),                      
y=L^{\mathrm{T}}x,    \quad \quad \quad \quad (1)           
</math>                                              

where <math>s_1, s_2, \ldots, s_{p}</math> are the parameters of the system. They can be any scalar functions of some source parameters, like <math>s_1=e^t</math>, where <math>t</math> is time, or combination of several physical parameters like <math>s_1=\rho v</math>, where <math>\rho</math> and <math>v</math> are two physical parameters.   

<math>x(t)\in \mathbb{R}^n</math> is the state vector, <math>u \in \mathbb{R}^{d_I}</math> and <math>y \in
\mathbb{R}^{d_O}</math> are, respectively, the inputs and outputs of the
system. To obtain the reduced model in (1), a
projection matrix <math>V</math> which is independent of all the parameters has
to be computed.

<math>V^T(E_0+s_1E_1+s_2E_2+\ldots +s_pE_p)Vx=V^TBu(s_p),
y=L^{\mathrm{T}}Vx.
</math>

The matrix $V$ is derived by orthogonalizing a number of moment
matrices of the system in (1)[1][2].

By defining <math>B_M=\tilde{E}^{-1}B, M_i=-\tilde{E}^{-1}E_i,i=1,2,\ldots,p</math> and

<math>
\tilde{E}=E_0+s_1^0E_1+s_2^0E_2+\cdots+s_p^0E_p,
</math>

we can expand <math>x</math> in (1) at <math>s_1, s_2, \ldots, s_p</math> around a set of
expansion points <math>p_0=[s_1^0,s_2^0,\cdots,s_p^0]</math> as below,

<math>
 x=[I-(\sigma_1M_1+\ldots +\sigma_pM_p]^{-1}B_Mu(s_p)
 =\sum\limits_{i=0}^{\infty}(\sigma_1M_1+\ldots+\sigma_pM_p)^iB_Mu(s_p).
</math>

Here <math>\sigma_i=s_i-s_i^0, i=1,2,\ldots,p</math>. We call the coefficients
in the above series expansion moment matrices of the parametrized
system, i.e. <math>B_M, M_1B_M, \ldots, M_pB_M, M_1^2B_M, (M_1M_2+M_2M_1)B_M, \ldots, (M_1M_p+M_pM_1)B_M, M_p^2B_M, M_1^3B_M, \ldots</math>. The corresponding moments are those moment
matrices multiplied by <math>L^{\mathrm{T}}</math> from the left. The matrix <math>V</math> can be
generated by first explicitly computing some of the moment matrices
and then orthogonalizing them as is suggested in~\cite{Daniel04}.
The resulting <math>V</math> is desired to expand the subspace:

<math>
\mathop{\mathrm{range}}\{V\}=\mathop{\mathrm{span}}\{B_M, \ M_1B_M,\ldots, M_pB_M,\ M_1^2B_M, </math>
 
<math>(M_1M_2+M_2M_1)B_M, \ldots, (M_1M_p+M_pM_1)B_M, 
 M_p^2B_M, M_1^3B_M,\ldots, M_1^rB_M, \ldots,M_p^rB_M \}.        \quad \quad \quad   \quad      (2)
</math>

However, <math>V</math> does not really span the whole subspace, because the
latterly computed vectors in the subspace become linearly dependent
due to numerical instability. Therefore, with this matrix <math>V</math> one
cannot get an accurate reduced model which matches all the moments
included in the subspace.

Instead of directly computing the moment matrices in (2)[1], a
numerically robust method is proposed in [1] (the
detailed algorithm is described in [3]), which combines
the recursions in (4) with the modified Gram-Schmidt
process to implicitly compute the moment matrices. The computed <math>V</math>
is actually an orthonormal basis of the subspace as below,

<math>
\mathop{\mathrm{range}}\{V\}=\mathop{\mathrm{span}}\{R_0, R_1,\ldots, R_r \}.  \quad \quad \quad  \quad (3)
</math>


<math> R_0 </math>=[ <math>B_M  </math> ],

<math>R_1=[M_1R_0,\ldots, M_pR_0], </math>

<math>R_2=[M_1R_1,\ldots, M_pR_1], \quad \quad \quad \quad (4) </math>                      

<math>\vdots,</math> 

<math>R_r=[M_1R_{r-1},\ldots, M_pR_{r-1}]</math>

<math>\vdots.</math>

It can be proved that the subspace in(2) is included in the
subspace in(3). Due to the numerical stability properties of
the repeated modified Gram-Schmidt process employed in
[2][3], the reduced model derived from <math>V</math>
in(3) is computed in a numerically stable and accurate way.