AO short=AO, long=Adaptive Optics, GMT short=GMT, long=Giant Magellan
Telescope, ELT short=ELT, long=Extremely Large Telescope, ELTs short=ELTs,
long=Extremely Large Telescopes, VCM short=VCM, long=Voice Coil Motor, FE
short=FE, long=Finite Element, ODEs short=ODEs, long=Ordinary Differential
Equations, LTI short=LTI, long=Linear Time Invariant, SISO short=SISO,
long=Single Input Single Output, MIMO short=MIMO, long=Multiple Input Multiple
Output, SVD short=SVD, long=Singular Value Decomposition, IRKA short=IRKA,
long= Iterative Rational Krylov Algorithm, ITIA short=ITIA, long= Iterative
Tangential Interpolation Algorithm , ISTIA short=ISTIA, long= Iterative SVD
Tangential Interpolation Algorithm , FISTIA short=FISTIA, long= Frequency
Limited Iterative SVD Tangential Interpolation Algorithm , RB short=RB,
long=Reference Body, HO short=HO, long=Higher Order, LF short=LF, long=Loewner
Framework, BT short=BT, long=Balanced Truncation,

# High fidelity adaptive mirror simulations
with reduced order models

Bernadett Stadler RICAM, Altenbergerstraße 69, A-4040 Linz, Austria,
<EMAIL_ADDRESS>ronny.ramlau@ricam.oeaw.ac.at). Roberto
Biasi Microgate, Via Waltraud Gebert Deeg 3e, 39100, Bolzano, Italy,
<EMAIL_ADDRESS>mauro.manetti@microgate.it). Mauro
Manetti22footnotemark: 2 Andreas Obereder Industrial Mathematics Institute,
JKU Linz, Altenberger Str. 69, 4040 Linz, Austria,
<EMAIL_ADDRESS>ronny.ramlau@jku.at). Ronny
Ramlau11footnotemark: 1 33footnotemark: 3 Matteo Tintori A.D.S. International,
Via Pio Galli sindacalista 3, 23841, Annone di Brianza, Italy, (m.tintori@ads-
int.com).

# High fidelity adaptive mirror simulations
with reduced order models

Bernadett Stadler RICAM, Altenbergerstraße 69, A-4040 Linz, Austria,
<EMAIL_ADDRESS>ronny.ramlau@ricam.oeaw.ac.at). Roberto
Biasi Microgate, Via Waltraud Gebert Deeg 3e, 39100, Bolzano, Italy,
<EMAIL_ADDRESS>mauro.manetti@microgate.it). Mauro
Manetti22footnotemark: 2 Andreas Obereder Industrial Mathematics Institute,
JKU Linz, Altenberger Str. 69, 4040 Linz, Austria,
<EMAIL_ADDRESS>ronny.ramlau@jku.at). Ronny
Ramlau11footnotemark: 1 33footnotemark: 3 Matteo Tintori A.D.S. International,
Via Pio Galli sindacalista 3, 23841, Annone di Brianza, Italy, (m.tintori@ads-
int.com).

###### Abstract

In the design process of large adaptive mirrors numerical simulations
represent the first step to evaluate the system design compliance in terms of
performance, stability and robustness. For the next generation of Extremely
Large Telescopes increased system dimensions and bandwidths lead to the need
of modeling not only the deformable mirror alone, but also all the system
supporting structure or even the full telescope. The capability to perform the
simulations with an acceptable amount of time and computational resources is
highly dependent on finding appropriate methods to reduce the size of the
resulting dynamic models. In this paper we present a framework developed
together with the company Microgate to create a reduced order structural model
of a large adaptive mirror as a preprocessing step to the control system
simulations. The reduced dynamic model is then combined with the remaining
system components allowing to simulate the full adaptive mirror in a
computationally efficient way. We analyze the feasibility of our reduced
models for Microgate’s prototype of the adaptive mirror of the Giant Magellan
Telescope.

Keywords: model order reduction, modal truncation, balanced truncation, Krylov
subspace methods, moment matching, adaptive mirrors

AMS: 65-04, 74-10, 85-08, 85-10

##  1 Introduction

For ground-based telescopes, so-called AO systems are used to compensate the
image distortions in astronomical observations caused by atmospheric
turbulence, using the flexible shape of an optical surface. The control of an
AO system is complex and requires at least a wavefront sensor to get
information about the atmospheric turbulence causing the image distortions, a
known natural or laser guide star as reference source, a mirror shape command
generator, and a deformable mirror to compensate for the atmospheric
turbulence.

We consider the control of an adaptive mirror based on non-contacting voice-
coil actuators, which are co-located to capacitive position sensors [12]. This
deformable mirror technology was developed by the company
Microgate111https://engineering.microgate.it together with other partners (ADS
International222https://www.ads-int.com/, INAF-Osservatorio Astrofisico di
Arcetri333https://www.arcetri.inaf.it/en/, and the Aerospace Engineering
Department of Politecnico di Milano444https://www.aero.polimi.it/it/il-
dipartimento). It has been already deployed realizing the adaptive secondary
mirror of several large telescopes, including the Multiple Mirror Telescope,
the Large Binocular Telescope, and the Very Large Telescope. Microgate is
engaged in the final design and construction of the adaptive mirrors for the
next generation of ELTs. Sub-system and full-system multiphysics simulation
plays a key role in the design phase of such complex projects. Therefore,
numerical simulation has always been intensively pursued by the company.
Several design solutions need to be compared and the impact of different
modelling choices have to be verified.

The structural dynamic models used to describe the large adaptive mirror are
of high order, caused by the need to describe well high spatial order
deformations and by the inherent complexity of the system supporting
structure. Hence, a reduced order mirror model is required that guarantees
high fidelity results with a reasonable simulation time. The earliest methods
for model order reduction techniques go back to the 1960s in the field of
structural dynamics. These methods rely on the identification of
eigenfrequencies and are referred to as mode displacement methods [40, 22]. In
the 1980s the important reduced order method balanced truncation [35, 15] has
been developed in the system and control theory community. In the field of
numerical mathematics, approaches such as Pade-via-Lanczos and rational
interpolation methods came up in the 1990s and are still under research, e.g.
in [16, 27, 2]. Nowadays these methods are often used for the design and
analysis of large electronic circuits. In recent years data driven approaches
gained a lot of traction, see e.g. [8, 26, 24, 31, 4, 2, 25]. In the
literature various reviews on model order reduction techniques exist [13, 14,
28, 5, 17, 1, 3, 19, 11, 9, 10, 6].

In the framework of this paper a reduced order model is created by applying
different methods as a preprocessing step to the control system simulations.
In [33] the authors used balanced truncation for reducing the complexity of an
adaptive mirror model, which is a common approach for control theory
applications. However, in contrast to the model presented in this paper theirs
did not take into account the system supporting structure up to the full
telescope and was thus less complex. Moreover, here we also consider Krylov
subspace based methods and the Loewner framework. We combine the reduced order
structural model with the remaining system components and run full system
simulations using different numerical tools. Our developments are validated
via simulations of Microgate’s P72 [18], which is a $72$ actuation points
prototype of the GMT.

The outline of this paper is as follows: We start with a brief description of
the physical model of the adaptive mirror in Section 2. Section 3 is dedicated
to an overview of existing model order reduction methods from the field of
structural dynamics, control theory, and numerical mathematics. In Section 4
we describe our framework for performing high fidelity mirror simulations with
a reduced order structural model and we give some implementation details.
Numerical simulations for the GMT P72 adaptive mirror, including a performance
evaluation of different reduced order methods, are shown in Section 5. We end
with a summary and conclusions in Section 6.

##  2 Physical model

The physical modeling of the deformable mirrors treated in this paper usually
requires a multiphysics description. The level of modeling accuracy and
complexity can be set on the basis of the needs and usually involves the
description of (see Figure 1):

  * •

a thin deformable mirror, usually a Zerodur shell about 2 mm thick;

  * •

a reference structure called reference body, also made of Zerodur to grant the
required thermal stability;

  * •

the squeeze film action of the air trapped in the thin gap, about 100 microns,
between the mirror and the reference body, which significantly contributes to
the deformable mirror damping, affecting both the control loop stability and
the controlled system performance;

  * •

the cold plate, typically made by aluminum, where the VCM actuators are
mounted, which sustains the VCM reaction forces providing adequate mechanical
stability and taking care of the system cooling;

  * •

the system positioner, typically a hexapod, responsible for the system
alignment and offloading of the mirror low-order modes;

  * •

the mirror capacitive sensors and VCM actuators with their signal conditioning
and the digital feedforward-feedback control loop.

Figure 1: Deformable mirror and supporting structure.

The focus of this paper is on model order reduction methods for the structural
dynamics of the system components mentioned above, even if the effectiveness
of the reduced order models are verified also through the simulation of the
full system behavior. Details about the fluid dynamic model and the control
system are omitted and can be found in [34].

The deformable mirror is described using a linear structural model in the time
domain. One way to represent the structural dynamics is via its degrees of
freedoms using second order differential equations. Another way is given by
its system states with first order differential equations often referred to as
state space representation. This formulation is preferred by control engineers
as linear control system analyses and design methods are usually given in the
state space form. Besides the form of equations the model can be represented
via different coordinate systems. Commonly, nodal and modal coordinates are
employed. The nodal coordinates are expressed through displacements and
velocities of specific physical locations, also called nodes. In contrast,
modal coordinates are defined through structural eigenmodes [21].

###  2.1 Second order structural models

The structural model of the adaptive mirror shape can be represented by second
order linear differential equations using nodal coordinates, i.e., in terms of
displacement, velocity and acceleration. Typically, the FE method with
piecewise polynomial basis functions [44] is applied, where the structure is
described by the mass, stiffness and damping matrices. The discretized
displacement $u_{h}$, which depends on the spatial variable $z$ and the time
$t$, is represented via

$u_{h}(z,t)=\sum_{i=1}^{n_{g}}\psi_{i}(z)u_{i}(t),$ (1)

where $n_{g}$ is the number of degrees of freedom, $\psi_{i}(z)$ are the
linear independent FE basis functions, and
$(u_{i})_{i=1}^{n_{g}}\eqqcolon\vec{u}_{h}(t)\in\mathbb{R}^{n_{g}}$ are the
time dependent coefficients of the FE solution $u_{h}$. With Equation (1) we
can identify the FE function $u_{h}$ with its coefficient vector
$\vec{u}_{h}$. Assuming a linear structural model described by its degrees of
freedom the FE discretization leads to the following system of ODEs

$M_{g}\ddot{\vec{u}}_{h}+D_{g}\dot{\vec{u}}_{h}+K_{g}\vec{u}_{h}=B_{g}(f_{a}^{c}+f_{a}^{d}),$
(2)

where $M_{g}$ is the $n_{g}\times n_{g}$ mass matrix, $D_{g}$ is the
$n_{g}\times n_{g}$ damping matrix, $K_{g}$ is the $n_{g}\times n_{g}$
stiffness matrix and $B_{g}$ is the $n_{g}\times n_{a}$ force influence
matrix. The control and disturbance force vectors at the $n_{a}$ actuation
points are denoted by $f_{a}^{c}$ and $f_{a}^{d}$, respectively. For more
details about the mechanical model and the coupling with the fluid part we
refer to [33].

Equation (2) can be formulated in modal coordinates. When considering free
vibrations of a structure without damping, i.e., a structure without external
excitation and with a damping matrix $D=0$, the vibration modes are found
using a time harmonic representation of the displacement of the unforced
system in (2). This leads to the following generalized eigenvalue problem

$(K_{g}-\omega_{i}^{2}M_{g})\phi_{i}=0,$ (3)

where $\omega_{i}^{2}$ with $i\in\\{1,\dots,n_{m}\\}$ is an eigenfrequency and
$\phi_{i}$ the corresponding mode shape vector. The so-called modal
displacement can be written as

$q_{m}=\sum_{i=1}^{n_{m}}\phi_{i}q_{i},$ (4)

with the modal coefficients denoted by $q_{i}$. The displacement written in
matrix notation becomes

$q_{m}=\Phi_{g}q,\quad\Phi_{g}=[\phi_{1},\phi_{2},\dots,\phi_{n_{m}}]\in\mathbb{R}^{n_{g}\times
n_{m}},$ (5)

with the $n_{g}\times n_{m}$ modal shape matrix $\Phi_{g}$ normalized such
that

$\displaystyle M_{m}$ $\displaystyle\coloneqq\Phi_{g}^{T}M_{g}\Phi_{g}=I,$ (6)
$\displaystyle K_{m}$
$\displaystyle\coloneqq\Phi_{g}^{T}K_{g}\Phi_{g}=\operatorname{diag}(\omega_{i}^{2}).$
(7)

The adaptive mirrors considered in this paper experiences, at least most of
the time, only light damping. In this setting it is common, see e.g. [21], to
include a diagonal modal damping approximation with damping coefficients
$\xi_{i}$

$D_{m}\coloneqq\operatorname{diag}\\{2\xi_{i}\omega_{i}\\}\approx\Phi_{g}D_{g}\Phi_{g}^{T}.$
(8)

The force influence matrix is transformed similarly

$\displaystyle B_{m}\coloneqq\Phi_{g}^{T}B_{g}.$ (9)

In this way the equation of motion (2) in modal coordinates becomes fully
uncoupled

$\ddot{q_{i}}+2\xi_{i}\omega_{i}\dot{q}_{i}+\omega^{2}_{i}q_{i}=f_{m_{i}}^{c}+f_{m_{i}}^{d},\quad
i\in\\{1,\dots,n_{m}\\},$ (10)

where $f_{m_{i}}^{c}$ and $f_{m_{i}}^{d}$ denote the $i$-th component of the
modal force vectors

$\displaystyle f_{m}^{c}\coloneqq B_{m}f_{a}^{c},$ (11) $\displaystyle
f_{m}^{d}\coloneqq B_{m}f_{a}^{d}.$ (12)

Note that for the numerical simulations we focus on modal coordinates as the
FE model provided by the company is given in this representation.

###  2.2 State space representation

For control theory applications it is common to rewrite the second order
differential equations in (2) into a so-called state space representation,
which is a set of first order differential equations. We consider the LTI
state space representation of the form

$\displaystyle\dot{x}$ $\displaystyle=Ax+Bu$ (13) $\displaystyle y$
$\displaystyle=Cx,$ (14)

where $u\in\mathbb{R}^{m}$ is the input and $y\in\mathbb{R}^{p}$ the output of
the system. The state vector is denoted by $x\in\mathbb{R}^{n}$ and
$A\in\mathbb{R}^{n\times n}$ is commonly referred to as state matrix. The
input matrix is denoted by $B\in\mathbb{R}^{n\times m}$ and the output matrix
by $C\in\mathbb{R}^{p\times n}$.

By applying the Laplace transform the model can be represented via its
transfer function, which describes the input-output behavior of the system. In
this paper we assume to have a stable system, i.e., all eigenvalues of $A$
have a negative real part [21]. Moreover, we assume a minimal realization,
i.e., all components of the state vector contribute to the input-output
behavior. The transfer function of the system in (13)-(14) is given by

$H(s)=C(sI-A)^{-1}B,\quad s\in\mathbb{C}.$ (15)

Note that the transfer function is invariant under coordinate transformation.

The modal version of the state space representation (13)-(14) of flexible
structures has a special form, in which the state matrix $A_{m}$ in modal
coordinates is block diagonal

$\displaystyle A_{m}=\operatorname{diag}(A_{m_{i}}),$

with $A_{m_{i}}$ being $2\times 2$ blocks, see e.g. [21]. The input and output
matrices in modal coordinates are arranged accordingly

$B_{m}=[B_{m_{1}},\dots,B_{m_{n_{m}}}]^{T}\quad
C_{m}=[C_{m_{1}},\dots,C_{m_{n_{m}}}].$ (16)

The state vector is split into $n_{m}$ components that represent the state of
a certain mode

$x=[x_{1},\dots,x_{n_{m}}]^{T},\quad x_{i}=[x_{i_{1}},x_{i_{2}}].$ (17)

The $i$-th component, i.e., the $i$-th mode, has the state representation
$(A_{m_{i}},B_{m_{i}},C_{m_{i}})$.

We consider the following representation for the blocks of the state space
matrices

$A_{m_{i}}=\begin{bmatrix}0&1\\\
\omega_{i}&-2\xi_{i}\omega_{i}\end{bmatrix},\quad
B_{m_{i}}=\begin{bmatrix}0\\\ b_{m_{i}}\end{bmatrix},\quad
C_{m_{i}}=\begin{bmatrix}0\\\ c_{m_{i}}\end{bmatrix}.$ (18)

The $i$-th state is given by

$x_{i}=[q_{m_{i}},\dot{q}_{m_{i}}]^{T},$ (19)

where $q_{m_{i}}$ is the $i$-th modal displacement and $\dot{q}_{m_{i}}$ the
modal velocity. Each component has a modal displacement and velocity, which
are related to the original one by (4). Note that this leads also to a special
form of the transfer function in modal coordinates which we denote by
$H_{m}(s)$. For more details we refer to [21].

##  3 Model order reduction

In practice, we deal with a very large number of degrees of freedom for the
structural model required to obtain a certain modeling accuracy. Moreover, for
the design and construction of the large adaptive mirrors a high number of
simulations has to be performed. This makes the direct use of a FE model
infeasible and a reduced order model, which represents well the dynamics of
the structure within the frequency band of interest, has to be used.

For the control of an adaptive mirror it is important that the reduced order
model preserve the input-output behavior. Thus, the quality of the reduced
model can be evaluated by comparing its output to the output of the high order
model. This is commonly done by considering the error in the transfer
function. The $\mathcal{H}_{2}$ error of a stable system is defined by

$\left\|H_{r}-H\right\|_{\mathcal{H}_{2}},\quad\left\|H\right\|_{\mathcal{H}_{2}}=\sqrt{\frac{1}{2\pi}\int_{-\infty}^{\infty}\text{tr}\\{H^{T}(-j\omega)H(j\omega)\\}d\omega},$
(20)

where we denote by $j$ the imaginary unit and $H$ and $H_{r}$ are the transfer
functions of the high and reduced order model, respectively. Attention is also
paid to the $\mathcal{H}_{\infty}$ error

$\left\|H_{r}-H\right\|_{\mathcal{H}_{\infty}},\quad\left\|H\right\|_{\mathcal{H}_{\infty}}=\text{sup}_{\omega\in\mathbb{R}}\sigma_{\text{max}}(H(j\omega)),$
(21)

where $\sigma_{\text{max}}(.)$ is the largest singular value. There exist
several algorithms that aim for minimizing the $\mathcal{H}_{2}$ error, i.e.,
the so-called first order optimality condition,

$\left\|H_{r}-H\right\|_{\mathcal{H}_{2}}\rightarrow\text{min}\text{, for
stable }H_{r}.$ (22)

In the following we summarize model order reduction approaches based on modal
truncation, moment matching and data-driven interpolation.

###  3.1 Modal truncation

The modal displacement method [22] computes a reduced model based on free
vibration modes of the structure. For an acceptable approximation of the
dynamics of the mirror a relative small number of modes may be sufficient. To
this end not all $n_{m}$ eigenvectors in (4) are used, but a limited number
$r<n_{m}$. Usually, one is interested in the response of the system for lower
frequencies, as most structural forcing terms operate at low frequencies.
Hence, the first $r$ eigenvectors are kept, which correspond to the $r$ lowest
eigenfrequencies, and the others are truncated. Note that to accurately
control complex mirror deformations a number of modes at least equal to the
number of actuators is required to correctly represent the number of
independent degrees of freedom of the control problem.

###  3.2 Balanced truncation

Balanced truncation [35, 15] is a very popular model order reduction method
which allows to reach the first order optimality condition in (22). The LTI
system is transformed into a so-called balanced realization, in which the
states are ordered according to their contribution to the input-output
behavior. To quantify the systems input-output behavior the notions of
controllability and observability are introduced. These system properties are
determined via the so-called controllability $W_{c}\in\mathbb{R}^{n\times n}$
and observability $W_{o}\in\mathbb{R}^{n\times n}$ Gramians [35]. We are
interested in the stationary or time invariant solution. In this case the
controllability and observability Gramians, see e.g. [43], can be obtained by
solving the Lyapunov equations

$\displaystyle AW_{c}+W_{c}A^{T}+BB^{T}=0,$ (23) $\displaystyle
A^{T}W_{o}+W_{o}A+C^{T}C=0,$ (24)

where $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$ and
$C\in\mathbb{R}^{p\times n}$. The Hankel singular values $\gamma_{i}$ with
$i=1,\dots,n$, used to determine the importance of a certain state for the
input-output behavior, are defined as the square roots of the eigenvalues of
the product $W_{c}W_{o}$. Note that in contrast to Gramians the Hankel
singular values are system invariants. Once the system is transformed into the
balanced representation the states with small influence on the system
dynamics, i.e., with small Hankel singular values, are discarded. The reduced
system is again in balanced form, therefore, it preserves the stability
property of the original system [37]. This is a very important point for
adaptive mirror control. Moreover, an upper bound on the error [15, 23] is
given by

$\left\|H_{r}(s)-H(s)\right\|_{\mathcal{H}_{\infty}}\leq
2\sum_{i=r+1}^{n}\gamma_{i},$ (25)

i.e., we obtain a good approximation as long as the truncated Hankel singular
values are small.

In balanced truncation the total frequency range is covered, however, for
adaptive mirror control there is a certain frequency range of interest. For
the frequency limited balanced truncation, the controllability and
observability Gramians are computed in finite frequency intervals. We refer to
[20] for more details on the computation. Similar extensions exist to, e.g.,
time-limited balanced truncation or frequency weighted balanced truncation,
e.g. in [28], but are not considered in this paper.

####  3.2.1 Modal approximation

A huge drawback of balanced truncation is that solving the Lyapunov equations
in (23) is computationally very demanding. Thus, applying balanced truncation
in this form would not be feasible for our large scale system. However, the
representation of the state space in modal form (18) leads to a special form
of the controllability and observability Gramians [21]. Assuming small damping
the Gramians can be approximated via

$\displaystyle W_{c}\approx\operatorname{diag}(w_{c_{i}}I),\quad
w_{c_{i}}=\frac{\left\|B_{m_{i}}\right\|_{2}^{2}}{4\xi_{i}\omega_{i}},$ (26)
$\displaystyle W_{o}\approx\operatorname{diag}(w_{o_{i}}I),\quad
w_{o_{i}}=\frac{\left\|C_{m_{i}}\right\|_{2}^{2}}{4\xi_{i}\omega_{i}},$ (27)

where $w_{c_{i}}$ and $w_{o_{i}}$ are called modal controllability and
observability factors. The approximate Hankel singular values in modal
coordinates are the geometric means of these two factors

$\gamma_{i}\approx\sqrt{w_{c_{ii}}w_{o_{ii}}}.$ (28)

Equations (26)-(27) can be evaluated very fast. Moreover, the system preserves
its original modal representation.

More recently also a data-driven approach of balanced truncation has been
developed in [24] in which the product of the Gramians is approximated
directly through frequency response data up to a desired accuracy.

###  3.3 Moment matching methods

Similar to balanced truncation, the goal of moment matching methods is to
reduce the system to one with fewer degrees of freedom, but with similar
input-output behavior, i.e., approximating the transfer function. Typically,
those methods are used for large electronic systems, but also for structural
vibrations, e.g., in [32]. The basic approach is to approximate the state
space $S$ of $x$ by a low dimensional subspace $\tilde{S}$ via projection.

For an arbitrary interpolation point $s_{0}\in\mathbb{C}$ the transfer
function can be written via its moment expansion

$\displaystyle H(s)=\sum_{i=0}^{\infty}(-1)^{i}M_{i}(s_{0})(s-s_{0})^{i},$
(29) $\displaystyle M_{i}(s_{0})\coloneqq C^{T}(s_{0}I-A)^{-i}(sI-A)^{-1}B.$
(30)

where the $M_{i}$ are the so-called moments. Note that for Multiple Input
Multiple Output (MIMO) systems the moments are $p\times m$ matrices. The idea
behind moment matching methods is to truncate the Taylor series in (29) and
thus find a reduced order transfer function $H_{r}$, which matches the first
$r$ moments of $H$ at the given expansion point $s_{0}$. The precision of the
moment matching methods depends not only on the number of moments matched but
also highly depends on the chosen expansion point $s_{0}$. The reason is that
the Taylor series in (29) is only a reasonable approximation within a certain
distance from $s_{0}$. To increase accuracy methods that use multiple
expansion points $s_{1}\dots,s_{\ell}\in\mathbb{C}$ are often used. In the
SISO case, the moment matching problem is given by finding an $H_{r}$ such
that

$H_{r}(s_{i})=H(s_{i})\quad\text{ for }i=1,\dots,N,$ (31)

which is also known as rational interpolation. In the MIMO case the moments
are matrices and one commonly interpolates along certain directions, which
leads to the so-called tangential interpolation problem, see e.g. [42]. The
aim is to find a rational matrix $H_{r}$ for given interpolation points
$\lambda_{i}\in\mathbb{C}$ and $\mu_{i}\in\mathbb{C}$ and the corresponding
right and left tangential directions $r_{i}\in\mathbb{C}^{m}$ and
$\ell_{i}\in\mathbb{C}^{p}$ such that

$\displaystyle H_{r}(\lambda_{i})r_{i}=H(\lambda_{i})r_{i}\quad\text{ for
}i=1,\dots,k,$ (32)
$\displaystyle\ell_{j}^{T}H_{r}(\mu_{j})=\ell_{j}^{T}H(\mu_{j})\quad\text{ for
}j=1,\dots,q.$ (33)

The interpolation points and tangential directions have to be selected in
advance to realize certain model reduction goals.

####  3.3.1 Rational Krylov subspace methods

A common way to deal with the moment matching problem is rational
interpolation by projection. Here the moments are used to construct the
biorthogonal projection matrices $V,W\in\mathbb{R}^{n\times r}$, with which
the reduced order model is obtained via Petrov-Galerkin projection

$\displaystyle A_{r}=W^{T}AV,\quad B_{r}$ $\displaystyle=W^{T}B,\quad
C_{r}=C^{T}V,$ (34) $\displaystyle x$ $\displaystyle\approx Vx_{r},$ (35)

where $A_{r}\in\mathbb{R}^{r\times r}$, $B_{r}\in\mathbb{R}^{r\times m}$,
$C_{r}\in\mathbb{R}^{p\times r}$ and $x_{r}\in\mathbb{R}^{r}$. The columns of
$V$ form an orthonormal basis of a subspace $\tilde{S}\subset S$ of the space
of the state vector. Assuming the SISO case with a single expansion point
$s_{0}\in\mathbb{C}$, $\tilde{S}$ is spanned by

$\displaystyle\text{range}\\{V\\}=\text{span}\\{(s_{0}I-A)^{-1}B,(s_{0}I-A)^{-1}(s_{0}I-A)^{-1}B,\dots,$
(36) $\displaystyle((s_{0}I-A)^{-1})^{r-1}(s_{0}I-A)^{-1}B\\}.$

For the right projection matrix $W$ a similar property holds

$\displaystyle\text{range}\\{W\\}=\text{span}\\{(s_{0}I-A)^{-T}C^{T},(s_{0}I-A)^{-T}(s_{0}I-A)^{-T}C^{T},\dots,$
(37) $\displaystyle((s_{0}I-A)^{-T})^{r-1}(s_{0}I-A)^{-1}C^{T}\\}.$

The range of $V$ and the range of $W$ span so-called Krylov subspaces. Note
that when multiple expansion points $s_{1}\dots,s_{\ell}\in\mathbb{C}$ are
used, the Krylov subspaces get extended. In the SISO case $B$ and $C$ are
vectors, thus, for the construction of $V$ and $W$ rational Krylov subspace
algorithms, such as rational Lanczos or rational Arnoldi, can be used [27].
These methods are very efficient as only matrix factorization, forward and
backward substitution as well as matrix-vector multiplications are involved.
Their complexity is of $\mathcal{O}(nr^{2})$. For adaptive mirror simulations
we are dealing with the MIMO case, where $B$ and $C$ are matrices and the
tangential directions need to be included into Krylov subspaces. Here, block
rational Arnoldi or Lanczos methods can be employed.

Recently methods that start from an initial group of interpolation points and
iteratively update them have been developed. One of those methods is the IRKA.
Initially this algorithm was developed for SISO systems in [29]. In [42] it
was extended to the MIMO case called the ITIA. This method creates a reduced
order model that fulfills the first order $\mathcal{H}_{2}$ optimality
condition in (22). The updated expansion points are selected as the negative
poles of the transfer function of the reduced model $H_{r}$ or equivalently as
the negative eigenvalues of the state matrix[42].

In contrast to balanced truncation, moment matching methods, in general, do
not preserve the stability of the system. For systems with special structures
there exist approaches where the reduced order model is guaranteed to be
stable, see [36]. In [39] the authors proposed the ISTIA, which combines
balanced truncation and Krylov subspace based methods. For one of the
projection matrices the observability or controllability Gramian is used,
whereas the other projector is computed using Krylov subspace based methods.
The authors show that due to the usage of a Gramian preserves stability. When
frequency limited controllability or observability Gramians are used, the
method is commonly referred to as FISTIA [39].

An advantage of the Krylov subspace based methods it that their efficiency
does not rely on a modal representation, but only depends on the system size.
The more moment vectors are included into the subspace, the more accurate the
approximation is, however, the reduced order model is enlarged. In practice, a
trade-off between accuracy and the reduction size has to be made. If both
projection matrices are constructed via Krylov subspace methods it is commonly
referred to as two sided, whereas if only one of the projection matrices is
constructed in this way it is called one sided. Using the two sided method
results in $2r$ matched moments, whereas the one sided method matches $r$
moments [7].

####  3.3.2 Loewner framework

The Loewner framework is a data-driven interpolation method for which the
frequency response data, i.e., samples of the transfer function $H$ at certain
frequency points $s$, is used. In the MIMO case, where the moments are
$p\times m$ matrices, the tangential directions are included into the data set
and the data is split into a right part

$\displaystyle
M=\operatorname{diag}(\mu_{1},\dots,\mu_{q})\in\mathbb{C}^{q\times q},\quad
L=[\ell_{1}^{T},\dots,\ell_{q}^{T}]^{T}\in\mathbb{C}^{q\times p},$ (38)
$\displaystyle\quad
V=[v_{1}^{T},\dots,v_{q}^{T}]^{T}=[\ell_{1}^{T}H(\mu_{1}),\dots,\ell_{q}^{T}H(\mu_{q})]^{T}\in\mathbb{C}^{q\times
m},$

and a left part

$\displaystyle\Lambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{k})\in\mathbb{C}^{k\times
k},\quad R=[r_{1}^{T},\dots,r_{k}^{T}]^{T}\in\mathbb{C}^{m\times k},$ (39)
$\displaystyle\quad
W=[w_{1}^{T},\dots,w_{k}^{T}]^{T}=[H(\lambda_{1})r_{1},\dots,H(\lambda_{q})r_{k}]^{T}\in\mathbb{C}^{p\times
k}.$

The aim is to find a rational function $H_{r}$ that fulfills the conditions in
(32). The data matrices defined in Equation (38) and Equation (39), as well as
the Loewner matrix $\mathbb{L}$ and its shifted version $\mathbb{L}_{s}$
defined via

$\mathbb{L}:=\left(\frac{v_{i}^{T}r_{i}-\ell_{j}w_{j}}{\mu_{i}-\lambda_{j}}\right)_{i=1,\dots,q}^{j=1,\dots,k},\quad\mathbb{L}_{s}:=\left(\frac{\mu_{i}v_{i}^{T}r_{i}-\lambda_{j}\ell_{j}^{T}w_{j}}{\mu_{i}-\lambda_{j}}\right)_{i=1,\dots,q}^{j=1,\dots,k},$
(40)

are used to construct the unprocessed Loewner model
$\\{W,\mathbb{L},\mathbb{L}_{s},V\\}$. If the pencil
$(\mathbb{L},\mathbb{L}_{s})$ is regular,
$H_{r}(s):=W(\mathbb{L}_{s}-s\mathbb{L})^{-1}V$ satisfies the interpolation
condition in (32). In many applications the pencil
$(\mathbb{L},\mathbb{L}_{s})$ is singular and a post-processing step is
required. The dominant features of the data are extracted and redundancies
removed. Commonly this is done by applying the singular value decomposition
(SVD) to the augmented Loewner matrix

$\begin{bmatrix}\mathbb{L}&\mathbb{L}_{s}\end{bmatrix}=Y\tilde{\Sigma}\tilde{X}^{H},\quad\begin{bmatrix}\mathbb{L}\\\
\mathbb{L}_{s}\end{bmatrix}=\tilde{Y}\Sigma X^{H},$ (41)

with $\Sigma,\tilde{\Sigma}\in\mathbb{C}^{r\times r}$,
$Y\in\mathbb{C}^{q\times r}$, $\tilde{Y}\in\mathbb{C}^{2q\times r}$,
$X\in\mathbb{C}^{k\times r}$ and $\tilde{X}\in\mathbb{C}^{r\times 2k}$. The
truncation index $r$ is chosen depending on the application and data size. The
projected Loewner model is then given by

$\quad A_{r}=-Y^{T}\mathbb{L}_{s}X,\quad B_{r}=Y^{T}V,\quad C_{r}=WX.$ (42)

The SVD provides optimal low-rank solutions, however, the full SVD has cubic
complexity. An alternative way was proposed in [30] and uses the CUR
decomposition instead, which is less accurate, but has a lower asymptotic time
complexity. Moreover, it has the benefit that the rows and columns of the
decomposed matrices $\mathcal{C}$ and $\mathcal{R}$ are interpretable. In our
case, this means that the dominant interpolation points $\lambda$ and $\mu$
are given directly by the columns and rows of $\mathcal{C}$ and $\mathcal{R}$.

There exists an iterative extension of the Loewner framework called AAA
algorithm [4]. Within this method the fitted rational approximants are
expressed in a numerically stable way using a baricentric representation. In
every iteration the next interpolation points are selected via a greedy
method, i.e. interpolation is enforced at data points where the error was
maximal. Recently a version of the algorithm enforcing real-valued and
strictly proper rational approximants was proposed in [25].

##  4 Reduced order high fidelity mirror simulations

In the process of designing a specific adaptive mirror there is the need to
perform simulations that accurately represent the structural dynamics, i.e.,
realistic operative conditions. The simulations have to guarantee high
fidelity results, but with an affordable computational load. As the size of
the models is very large, model order reduction methods are required.

Figure 2 illustrates the process of performing high fidelity adaptive mirrors
simulations with reduced order models. In a first step the accurate and
complex FE models describing the deformable mirror, all the system supporting
structure and in some cases even the full telescope are created. These models
have been provided by the company A.D.S. International555https://www.ads-
int.com/. The FE models are then reduced by modal truncation to a predefined
frequency range of interest for the certain adaptive mirror. As a second step
of model order reduction, methods described in Section 3 are used to reduce
the system size further. The resulting models are then combined with the rest
of the system, i.e., fluid dynamics and control system, and the dynamic
analysis of the mirror is performed. Especially, the second step of model
order reduction makes it possible for Microgate to run their many simulations
required in the adaptive mirror’s design process in a feasible time frame.

Figure 2: Graphical illustration of performing high fidelity adaptive mirror
simulations with reduced order models.

In the upcoming sections the model order reduction and system simulations are
described in more detail.

###  4.1 Creating the reduced order mirror model

The FE model is given in modal form (4). The first step of model reduction is
applied by truncation to a given frequency range of the mirror as described in
Section 3.1. Based on the resulting reduced system the LTI state space
representation in modal form is set up as shown in Equation (18). A second
step of model order reduction is applied using either balanced truncation,
rational Krylov subspace methods or the Loewner framework. More details on the
implementation of the methods is given in the upcoming sections. Our code is
mainly implemented in MATLAB. Computationally demanding parts are programmed
in C. Parts of the C implementation use the Matrix Equation Sparse Solver
(M.E.S.S) library [41] to efficiently deal with matrices, factorization and
solving systems of equations.

####  4.1.1 Balanced truncation

Algorithm 1 shows the steps to perform the model order reduction with balanced
truncation. In order to make the method feasible for the given system size,
the modal approximation as shown in Equation (26) and Equation (27) is used.
The modal controllability factor $w_{c_{i}}$ and observability factor
$w_{o_{i}}$ and the Hankel singular values $\gamma_{i}$ are computed per mode
$i=1,\dots,n$. The Hankel singular values are then ordered and all values
below a certain threshold removed. As there is a one-to-one mapping between
the Hankel singular values and the modes, the corresponding modes are
discarded and in this way the system dimension reduced.

Algorithm 1 Balanced truncation in modal form

1:function ModalBT($B_{m},C_{m},\xi,\omega,t$)

2: for each mode $i=1,\dots,n$ do

3: Compute the modal controllability factors $w_{c_{i}}$ via Equation (26).

4: Compute the modal observability factors $w_{o_{i}}$ via Equation (27).

5:

6: Compute the Hankel singular values $\gamma_{i}$ using Equation (28).

7:

8: Order the Hankel singular values
$\gamma_{1}\leq\gamma_{2}\dots\leq\gamma_{n}$.

9: Discard all Hankel singular values below the given threshold $t$

10: $\gamma_{1}\leq\gamma_{2}\dots\leq\gamma_{j-1}$ where $\gamma_{j-1}\leq
t<\gamma_{j}$.

11:

12: Compute $B_{m}^{r},C_{m}^{r},\xi^{r}$ and $\omega^{r}$ from the higher
order model by discarding

13: all modes $i$ corresponding to the Hankel singular values
$\gamma_{1},\dots,\gamma_{j-1}$.

####  4.1.2 Rational Krylov subspace methods

Algorithm 2 illustrates the model order reduction procedure with Krylov
subspace based methods. Our studies showed that the quality of these methods
depends highly on the chosen interpolation points and that methods which use
multiple expansion points and update them iteratively perform best. Hence, we
consider here ITIA,ISTIA andFISTIA, see Section 3 for more details. In a first
step of the algorithm the initial interpolation points are chosen. The choice
of the points is important as the methods provide only local convergence. The
eigenvalues of the state matrix would be a good initial guess. However,
computing them for a large scale system is not feasible. Because our input FE
model is given in modal form, we use their modal approximation, see e.g. [21].
In general, one would choose the initial interpolation directions $b$ and $c$
as the corresponding left and right eigenvectors. Due to computational reasons
we omit the computation and set $b$ and $c$ to vectors of ones. The projection
matrices $V$ and $W$ are computed and used to obtain the reduced state space
system $(A_{r},B_{r},C_{r})$. For details on the computation of the projection
matrices we refer to [42, 39]. The interpolation points for the next iteration
are chosen as the eigenvalues of the state matrix $\lambda(A_{r})$ and the
interpolation directions as the corresponding left and right eigenvectors. In
this way the interpolation points and directions are updated in an iterative
way until a maximum number of iterations is reached. Although the input to the
reduced order methods is in modal coordinates, the Krylov subspace based
methods produce a general state space representation. In order to study the
system performance in an efficient way we transform these state space matrices
into modal form.

Algorithm 2 Rational Krylov subspace method

1:function RationalKrylov($A_{m},B_{m},C_{m},maxIter$)

2:

3: for each mode $i=1,\dots,n$ do

4: Compute the interpolation points $s_{i}=\xi_{i}\omega_{i}\pm
j\omega_{i}\sqrt{1-\xi^{2}_{i}}$.

5: Set up the interpolation directions $b_{i}=[1,\dots,1],c_{i}=[1,\dots,1]$.

6:

7: for $i=1,\dots,maxIter$ do

8: Compute the projection matrices $V$ and $W$.

9: Compute the reduced system via Equation (34).

10: Compute eigenvalues $\lambda$ and left and right eigenvectors
$U_{\ell},U_{r}$ of $A_{r}$.

11: Extract the new interpolation points $s=-\lambda$.

12: Extract the new interpolation directions $b=U_{\ell}$, $c=U_{r}$.

13:

14: Verify that the reduced system $(A_{r},B_{r},C_{r})$ is stable.

15: Transform $(A_{r},B_{r},C_{r})$ into modal coordinates
$(A_{m}^{r},B_{m}^{r},C_{m}^{r})$.

####  4.1.3 Loewner framework

Algorithm 3 shows how to create a reduced order model using the Loewner
framework. In a first step $N$ logarithmically spaced interpolation points are
created over the frequency range of interest. The reason for choosing
logarithmically distributed points is that in adaptive mirror control it is
more important to match the lower frequencies. Hence, having a dense sampling
grid in the beginning and a more sparse distribution towards the end of the
frequency range is more appropriate for our application. The initial
interpolation directions $b$ and $c$ are chosen as random values. Using the
input state space system in modal coordinates $(A_{m},B_{m},C_{m})$, the
transfer function $H_{m}(s)$ is constructed and the right and left data set is
computed as shown in Equation (38) and Equation (39). Note that as the state
space system is given in modal coordinates the evaluation of the transfer
function and thus the computation of the input-output pairs for the data set
is much faster. For more details we refer to [20]. The data set is then used
to compute the Loewner matrix $\mathbb{L}$ and its shifted version
$\mathbb{L}_{s}$ using Equation (40). The constructed complex representation
$(W,\mathbb{L},\mathbb{L}_{s},V)$ of an underlying dynamical system is
transformed into a real model
$(W^{r},\mathbb{L}^{r},\mathbb{L}_{s}^{r},V^{r})$, see e.g. [30]. To identify
dominant subsets of the data on which interpolation is enforced we apply the
SVD to the Loewner matrices using Equation (41). Note that the CUR
decomposition would provide a faster way of achieving this, however, it is
only an approximation and performed worse in our experiments. The reduced
order model is then obtained through projection using Equation (42). Similar
as for the Krylov subspace based method we transform the system into modal
coordinates as a last step.

Algorithm 3 Loewner framework with SVD

1:function LoewnerSVD($A_{m},B_{m},C_{m},N$)

2: Compute the interpolation points $s=\text{logspace}(-1,3,N)$.

3: Set up interpolation directions $b=\text{rand}(n)$ and $c=\text{rand}(n)$.

4:

5: Compute the transfer function $H_{m}$.

6: Construct the right data set $M,L$ and $V$ using Equation (38).

7: Construct the left data set $\Lambda,R,W$ via Equation (39).

8:

9: Compute the Loewner pencil $(\mathbb{L},\mathbb{L}_{s})$ via Equation (40).

10: Transform the complex system into a real system
$(W^{r},\mathbb{L}^{r},\mathbb{L}_{s}^{r},V^{r})$.

11:

12: Compute the SVD of the Loewner matrix using Equation (41).

13: Compute the reduced model $(A_{r},B_{r},C_{r})$ via Equation (42).

14:

15: Verify that the reduced system $(A_{r},B_{r},C_{r})$ is stable.

16: Transform $(A_{r},B_{r},C_{r})$ into modal coordinates
$(A_{m}^{r},B_{m}^{r},C_{m}^{r})$.

In our experiments we also considered the AAA algorithm [4], which is an
adaptive and iterative extension of the Loewner framework that chooses
automatically the reduction size for a given accuracy. However, for our
simulation setting the results have been worse compared to the Loewner
framework and thus we omit more details on this method.

###  4.2 System simulation

The reduced order model of the system structural dynamics are combined with
the remaining system and exploited to run dynamic analysis. Different
numerical tools have been developed for this purpose:

  * •

dynamic simulation to evaluate the system behavior in time, this simulation is
mainly focused on the analysis of the system performance and stability;

  * •

frequency response to have a characterization in the frequency domain, very
useful to understand possible critical frequencies in the dynamic response of
the system, but also to provide a measurement of system robustness;

  * •

root locus analysis to obtain a description of the system in the Laplace
domain, having a characterization of the system stability and stability
sensitivity, i.e. robustness, with respect to some specific design parameters.

The three numerical tools above mentioned can be used all together to provide
a comprehensive overview of the system behavior in terms of both stability and
performance. The accuracy of the system description provided by these
different tools is usually not equivalent and the fidelity level of each one
can be tuned to find the best compromise between accuracy of the results and
computational time and resources necessary to perform the analysis.

In the framework of this paper Microgate’s mirror time simulator has been
upgraded to a new version, based on a multi-threaded parallel C++
implementation. This new code can provide the most reliable system
multiphysics simulation currently available and contains all the features of
the previous implementation written in C language [33, 34]. The new simulator
has the additional relevant capability to take into account not only the
structural dynamics of the deformable shell, but of the whole system
supporting structure, potentially up to the telescope one. The order of the
structural dynamic description must be reduced as much as possible to allow
fully exploiting this new feature of the simulator, in order to prevent that
the simulation time and/or the amount of computational resources become a show
stopper for the analysis. The codes computational efficiency has been further
improved by using the blaze666https://github.com/parsa/blaze library, which is
an open-source, high-performance, C++ math library for dense and sparse
arithmetic, and OpenMP777https://www.openmp.org/ for parallelization.

##  5 Numerical analysis

To validate our developments we perform numerical analysis for the model of
the P72 prototype of the GMT adaptive secondary mirror system [18]. The GMT is
currently under construction in the Atacama desert in Chile and will become
one of the new ELTs [18]. The prototype has $72$ actuators and a diameter of
$354$ mm, featuring the four innermost rings of actuators of the on-axis
adaptive secondary mirror segment of the GMT, see Figure 3. The shell is $2$
mm thick and made of Zerodur, the lateral flexures connecting the shell to the
RB are the same of the final unit and also the VCM actuators are exactly the
final ones. The RB is made in Zerodur, while the cold plate, where the
electronics is mounted, is in aluminum. The cooling system is based on the
direct expansion gas cooling concept recently introduced in this application
field[38]. All the simulators mentioned in Section 4.2 exploit the very same
FE model, shown in Figure 3, to capture the system structural behavior,
through the use of eigenmodes.

Figure 3: Picture (left) and FE model (right) of the GMT P72 [18].

###  5.1 Dynamic model order reduction

The FE model of the GMT P72 prototype is given in modal form and we truncate
it to the pre-known frequency range $0-6$kHz. Based on the resulting reduced
system we set up the LTI state space model and compare BT, ITIA,ISTIA and the
LF. The HO model after modal truncation has a state matrix with dimension
$1672\times 1672$. For ISTIA we use the controllability Gramian as one of the
projection matrices. This choice provided slightly better results than using
the observability Gramian.

Figure 4 shows a logarithmic plot of the relative $\mathcal{H}_{\infty}$ error
of the transfer functions between the HO and the reduced order models. In plot
(a) we show the model reduction error for a reduced state matrix of dimension
$144\times 144$ and in (b) for a size of $330\times 330$. The reason for
choosing $144$ as lower limit is that the P72 consists of $72$ degrees of
freedom resulting in a state matrix of $144\times 144$. Hence, this is the
lowest order allowed to be able to statically set the mirror deformation over
$72$ actuation points. During our studies we obtained the best trade-off
between quality and speed with a state matrix of size $330\times 330$. We
compare here BT marked in blue, ITIA in orange, ISTIA in yellow and LF in
green. Note that for this setting ISTIA and FISTIA performed very similar,
hence, we omit FISTIA. The reason for this might be that the model was already
truncated to the desired frequency range by modal truncation. In Table 1 we
list the mean of the relative error over the whole frequency range. For both
simulated cases BT and the Krylov subspace methods performed very similar.
With LF we were not able to produce an appropriate reduced order model with a
size of $144\times 144$, hence, we do not show it in Figure 4(a). For a
reduced system size of $330\times 330$ LF yields the highest
$\mathcal{H}_{\infty}$ error. This might be caused by not optimally chosen
interpolation points. For the other methods an increased reduction size
provides a smaller mean relative error.

(a)

(b)

Figure 4: Logarithmic plot of the relative $\mathcal{H}_{\infty}$ error in dB between the HO and the reduced order transfer functions with reduced state matrices of dimension $144\times 144$ (left) and $330\times 330$ (right) for BT (blue), ITIA (orange), ISTIA (yellow) and LF (green). | Mean relative $\mathcal{H}_{\infty}$ error
---|---
Reduction size |  BT |  ITIA |  ISTIA |  LF
$144\times 144$ |  $0.627$ |  $0.768$ |  $0.691$ |  $-$
$330\times 330$ |  $0.158$ |  $0.153$ |  $0.368$ |  $0.556$
Table 1: Mean of the relative $\mathcal{H}_{\infty}$ error as shown in Figure
4. The mean is calculated over the desired frequency interval $0-6$kHz for BT
and the Krylovs subspace based methods ITIA and ISTIA and LF.

Throughout the simulations we observed that the results of the Krylov subspace
based methods and LF highly depend on the choice of the initial interpolation
points and directions. Note that for ITIA and ISTIA we use the modal
approximation of the eigenvalues of the HO state matrix. This offers an easy
way of running the algorithms without tuning it by hand. For LF we use $6000$
logarithmically spaced points.

In terms of run-time BT clearly outperforms the other methods, because
computing the modal approximation of the Hankel singular values is very fast.
For the Krylov subspace based methods the run-time depends mainly on the
number of iterations and the reduction size, since in every iteration the
eigenvalues and eigenvectors of the reduced order model have to be computed.
In our simulations we obtained convergence within less than $50$ iterations.
For the data-driven LF we observed that a fairly high amount of data points is
required to obtain an appropriate reduced order model, which increases the
computational time as well as the memory requirements significantly.

Further analysis showed that all reduced order methods are not able to
represent $72$ independent degrees of freedom with a state matrix of
$144\times 144$. Hence, this reduction size is too low to represent properly
the system behavior and we continue our analysis with a reduced order model of
size $330\times 330$. The performance of ISTIA and ITIA is very similar for
the P72 mirror, hence, we focus for the upcoming analysis on ITIA.

###  5.2 Performance of the system

In this section the reduced order models representing the system structural
dynamics are exploited to simulate the full system behavior and a comparison
of the results with the system using the original full set of eigenmodes is
done. The whole system modeling can be summarized by the following list:

  * •

The structural dynamics of the system, including: the deformable mirror, the
reference body, the cold plate and the P72 supporting structure.

  * •

The fluid dynamic modeling of the air trapped between the deformable mirror
and the reference body.

  * •

The inner control loop, controlling the mirror shape, including: the
deformable mirror control law, the voice coil motors dynamics, the capacitive
sensors dynamics, the digital and analog signals modeling.

The optical loop, which is responsible to generate the mirror commands, is not
taken into account within the present analysis, i.e, the mirror is simulated
as a stand-alone component of the AO system.

####  5.2.1 Stability

To analyze the stability of the control system we use root locus plots, i.e.,
we plot the poles of the transfer function in the Laplace domain. These plots
allow to study the sensitivity of the system against certain parameters. Here
we consider the feedback gain with a varying scaling factor. We use a circle
to indicate a stable system, whereas crosses for unstable systems. In terms of
model order reduction, it is important that the reduced model does not
indicate a stable system when the original model was not stable or vice versa.
Figure 5 shows the root locus plots of the higher order model (black) and of
reduced order models of size $330\times 330$, using BT (blue), ITIA (orange)
and LF (green). Note that only the important section around the $0$ real axis
is plotted, just for one half of the complex conjugate poles characterized by
a positive imaginary part. We observe that all models indicate a stable system
with gain scaling $1$ and an unstable system for a scaling of $2$. The
critical poles, i.e. the ones having positive real part, slightly change for
different methods. Our simulations showed that for smaller reduction sizes the
reduced model may lead to stable systems where the original model was
unstable.

Figure 5: Root locus plot of the HO model in black and the reduced order model
with size $330\times 330$ with BT in blue, ITIA in orange and LF in green. The
gain scaling is varied between $1$ and $2$. We show only the important section
around the real $0$ axis (dashed red line). Moreover, we show only the poles
with positive imaginary part and not their conjugate complex counterpart.
Circles indicate stable systems, whereas crosses indicate unstable systems.

####  5.2.2 Input-output behavior

To assess the difference in the input-output behavior of the full system, we
study the minimum and maximum singular values of the complementary sensitivity
transfer function matrix, defined to have the deformable mirror commands as
inputs and the mirror deformation at the actuation points as output. The
frequency response of the complementary sensitivity is commonly used to asses
the tracking capabilities of the closed-loop system as function of the command
frequency content. Figure 6 shows a logarithmic plot of the minimum (dashed)
and maximum (solid) singular values in dB for the HO system in black and the
reduced systems with BT in blue,ITIA in orange and LF in green. We observe
that the reduced order models obtained with BT and ITIA yield very similar
maximum and minimum singular values compared to the original system. In
contrast, the singular values produced by LF are significantly different.

Figure 6: Logarithmic plot of the minimum (dashed) and maximum (solid)
singular values in dB of the high and reduced order transfer functions in the
frequency interval $0-6$ kHz. We compare the HO model in black and the reduced
order models with BT in blue, ITIA in orange and LF in green.

###  5.3 Mirror simulator step response

The last system analysis uses the reduced order models as input to the C++
mirror simulator and studies the simulation over time of the control system
response to a certain mirror deformation command. In Figure 7 we show the step
response plots to a tilt (a) and trefoil (b) command. We omit the reduced
model produced by LF as it lead to too large forces for the voice coil motor
and thus the simulation was stopped. The plots on the left represent the shell
displacement, whereas the ones on the right show the modal control force. We
provide a zoom around the y-axis in order to be able to see the performance of
the reduced order methods BT in blue and ITIA in orange compared to the HO
model in black. We observe that for both mirror commands all methods provide
good results. However, the step response obtained with BT is more in line with
that of the original model, especially for the trefoil command case.

(a) Step response for a tilt command. The left plot shows the shell
displacement over time, whereas on the right the modal control force is shown.

(b) Step response for a trefoil command. The left plot shows the shell
displacement over time, whereas on the right the modal control force is shown.

Figure 7: Plot of the step response for a tilt command (a) and a trefoil
command (b) for the HO model in black and the reduced models with BT (BT) in
blue and ITIA in orange. In addition, a zoom around the y-axis is shown.

##  6 Conclusion

In the design phase of an adaptive mirror, numerical simulations are crucial
to evaluate the system design compliance. The FE models for the adaptive
mirrors of ELTs have a high complexity, hence, model order reduction
techniques that facilitate the computationally efficient analysis are
required. In this paper we present a framework for performing high fidelity
control system simulations within a reasonable time frame. In a preprocessing
step a reduced order model is created, which is then used to simulate the full
adaptive mirror. We perform a feasibility study of different reduced order
methods exploiting the numerical model of the P72 prototype of the GMT
secondary mirrors. Starting with the FE model in modal coordinates several
model order reduction algorithms are applied. The quality is analyzed by means
of the $\mathcal{H}_{\infty}$ error of the transfer functions. The reduced
order model of the structural dynamics is combined with the remaining system
modeling and its performance is analyzed regarding accuracy, stability and
robustness. BT and ITIA yield both very similar stable and accurate models
when reducing to a state matrix of size $330\times 330$ or larger. In terms of
computational load, BT with modal approximation is significantly faster than
the others. The data-driven LF was only able to produce an appropriate reduced
model when considering a larger reduction size. The implemented framework for
performing highly accurate simulation with reduced order models together with
the developed numerical analysis tools allow to improve the design process of
adaptive mirrors and the relevant research and development activity at
Microgate. Although the framework was tested here on a relatively small
prototype the developments can be directly applied to larger systems and will
allow simulations within a reasonable computational time.

## Acknowledgement

This research was funded in part by the Austrian Science Fund (FWF) SFB
10.55776/F68 “Tomography Across the Scales”, project F6805-N36 (Tomography in
Astronomy). For open access purposes, the authors have applied a CC BY public
copyright license to any author-accepted manuscript version arising from this
submission. Furthermore, the authors were supported by the Austrian Research
Promotion Agency (FFG) project number FO999888133, as well as the NVIDIA
Corporation Academic Hardware Grant Program. LW is partially supported by the
State of Upper Austria.

## References

  * [1] A. C. Antoulas.  Approximation of Large-Scale Dynamical Systems.  Society for Industrial and Applied Mathematics, 2005.
  * [2] A. C. Antoulas, C. A. Beattie, and S. Güğercin.  Interpolatory Methods for Model Reduction.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2020\.
  * [3] A.C. Antoulas.  An overview of approximation methods for large-scale dynamical systems.  Annual Reviews in Control, 29(2):181–190, 2005.
  * [4] Athanasios C. Antoulas and Brian. D. O. Anderson.  On the scalar rational interpolation problem.  Ima Journal of Mathematical Control and Information, 3:61–88, 1986\.
  * [5] Z. Bai.  Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems.  Applied Numerical Mathematics, 43(1):9–44, 2002.  19th Dundee Biennial Conference on Numerical Analysis.
  * [6] Christopher Beattie, Peter Benner, Mark Embree, Serkan Gugercin, and Sanda Lefteriu, editors.  Realization and Model Reduction of Dynamical Systems.  Springer Cham, 2022.
  * [7] P. Benner and L. Feng.  3 Model order reduction based on moment-matching, pages 57–96.  De Gruyter, Berlin, Boston, 2021.
  * [8] Peter Benner, Stefano Grivet-Talocia, Alfio Quarteroni, Gianluigi Rozza, Wil Schilders, and Luís Miguel Silveira, editors.  Volume 1 System- and Data-Driven Methods and Algorithms.  De Gruyter, Berlin, Boston, 2021.
  * [9] Peter Benner, Serkan Gugercin, and Karen Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM Review, 57(4):483–531, 2015.
  * [10] Peter Benner, Mario Ohlberger, Albert Cohen, and Karen Willcox.  Model Reduction and Approximation.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017\.
  * [11] B. Besselink, U. Tabak, A. Lutowska, N. van de Wouw, H. Nijmeijer, D.J. Rixen, M.E. Hochstenbach, and W.H.A. Schilders.  A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control.  Journal of Sound and Vibration, 332(19):4403–4422, 2013.
  * [12] Roberto Biasi, Daniele Gallieni, Piero Salinari, Armando Riccardi, and Paolo Mantegazza.  Contactless thin adaptive mirror technology: past, present, and future.  In Brent L. Ellerbroek, Michael Hart, Norbert Hubin, and Peter L. Wizinowich, editors, Adaptive Optics Systems II, volume 7736, page 77362B. International Society for Optics and Photonics, SPIE, 2010.
  * [13] R. R. Craig and M. C. C. Bampton.  Coupling of substructures for dynamic analyses.  AIAA Journal, 6(7):1313–1319, 1968.
  * [14] D. de Klerk, D. J. Rixen, and S. N. Voormeeren.  General framework for dynamic substructuring: History, review and classification of techniques.  AIAA Journal, 46(5):1169–1181, 2008.
  * [15] Dale F. Enns.  Model reduction with balanced realizations: An error bound and a frequency weighted generalization.  In The 23rd IEEE Conference on Decision and Control, pages 127–132, 1984.
  * [16] P. Feldmann and R.W. Freund.  Efficient linear circuit analysis by Pade approximation via the Lanczos process.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 14(5):639–649, 1995.
  * [17] Roland W. Freund.  Model reduction methods based on Krylov subspaces.  Acta Numerica, 12:267–319, 2003.
  * [18] Daniele Gallieni, Matteo Tintori, Pierluigi Fumi, Lorenzo Crimella, Roberto Biasi, Gerald Angerer, Mauro Manetti, Antonin Bouchez, Glenn Brossus, Frank Groark, and Peter M. Thompson.  GMT adaptive secondary mirrors subsystem final design.  In Laura Schreiber, Dirk Schmidt, and Elise Vernet, editors, Adaptive Optics Systems VII, volume 11448, page 114485K. SPIE, 2020.
  * [19] K. A. Gallivan, E. Grimme, and P. M. Van Dooren.  Model Reduction of Large-Scale Systems Rational Krylov Versus Balancing Techniques, pages 177–190.  Springer Netherlands, 1999.
  * [20] W. Gawronski and J. Juang.  Model reduction in limited time and frequency intervals.  International Journal of Systems Science, 21(2):349–376, 1990.
  * [21] W. K. Gawronski.  Advanced Structural Dynamics and Active Control of Structures.  Mechanical Engineering Series. Springer New York, NY, 2004.
  * [22] M. Geradin and D. J. Rixen.  Mechanical Vibrations: Theory and Application to Structural Dynamics.  John Wiley and Sons, 1997.
  * [23] K. Glover.  All optimal Hankel-norm approximations of linear multivariable systems and their $L_{\infty}$-error bounds.  International Journal of Control, 39(6):1115–1193, 1984.
  * [24] Ion Victor Gosea, Serkan Gugercin, and Christopher Beattie.  Data-driven balancing of linear dynamical systems.  SIAM Journal on Scientific Computing, 44(1):A554–A582, 2022.
  * [25] Ion Victor Gosea, Charles Poussot-Vassal, and Athanasios C. Antoulas.  On loewner data-driven control for infinite-dimensional systems.  In 2021 European Control Conference (ECC), pages 93–99, 2021.
  * [26] Ion Victor Gosea, Charles Poussot-Vassal, and Athanasios C. Antoulas.  Chapter 15 - data-driven modeling and control of large-scale dynamical systems in the loewner framework: Methodology and applications.  In Emmanuel Trélat and Enrique Zuazua, editors, Numerical Control: Part A, volume 23 of Handbook of Numerical Analysis, pages 499–530. Elsevier, 2022.
  * [27] E. J. Grimme.  Krylov projection methods for model reduction.  Dissertation, 1997.
  * [28] S. Gugercin and A. C. Antoulas.  A Survey of Model Reduction by Balanced Truncation and Some New Results.  International Journal of Control, 77(8):748–766, 2004.
  * [29] S. Gugercin, A. C. Antoulas, and C. Beattie.  $\mathcal{H}_{2}$ Model Reduction for Large-Scale Linear Dynamical Systems.  SIAM Journal on Matrix Analysis and Applications, 30(2):609–638, 2008.
  * [30] Dimitrios S. Karachalios, Ion Victor Gosea, and Athanasios C. Antoulas.  6 The Loewner framework for system identification and reduction, pages 181–228.  De Gruyter, Berlin, Boston, 2021.
  * [31] J. Nathan Kutz, Steven L. Brunton, Bingni W. Brunton, and Joshua L. Proctor.  Dynamic Mode Decomposition.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2016\.
  * [32] H. Liew and P. M. Pinsky.  Matrix-padé via lanczos solutions for vibrations of fluid–structure interaction.  International Journal for Numerical Methods in Engineering, 84(10):1183–1204, 2010.
  * [33] M. Manetti, M. Morandini, and P. Mantegazza.  High precision massive shape control of magnetically levitated adaptive mirrors.  Control Engineering Practice, 18(12):1386–1398, 2010.
  * [34] M. Manetti, M. Morandini, and P. Mantegazza.  Servo-fluid-elastic modeling of contactless levitated adaptive secondary mirrors.  Computational Mechanics, 50:85–98, 2012.
  * [35] B. Moore.  Principal component analysis in linear systems: Controllability, observability, and model reduction.  IEEE Transactions on Automatic Control, 26(1):17–32, 1981.
  * [36] A. Odabasioglu, M. Celik, and L.T. Pileggi.  PRIMA: passive reduced-order interconnect macromodeling algorithm.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 17(8):645–654, 1998.
  * [37] L. Pernebo and L. Silverman.  Model reduction via balanced state space representations.  IEEE Transactions on Automatic Control, 27(2):382–387, 1982.
  * [38] Dietrich Pescoller, Roberto Biasi, Matteo Tintori, Mauro Manetti, Marco Mantegazza, and Alberto Merler.  Direct expansion gas cooling technology applied to the E-ELT M4 adaptive mirror.  In AO4ELT5, 2017.
  * [39] C. Poussot-Vassal.  An iterative SVD-tangential interpolation method for medium-scale MIMO systems approximation with application on flexible aircraft.  In 2011 50th IEEE Conference on Decision and Control and European Control Conference, pages 7117–7122, 2011.
  * [40] J.W.S. Rayleigh.  The Theory of Sound, volume 1.  Dover Publications, 1945.
  * [41] J. Saak, M. Köhler, and P. Benner.  M-M.E.S.S.-2.1 – the matrix equations sparse solvers library, February 2022.  see also:https://www.mpi-magdeburg.mpg.de/projects/mess.
  * [42] P. Van Dooren, K.A. Gallivan, and P.-A. Absil.  H2-optimal model reduction of MIMO systems.  Applied Mathematics Letters, 21(12):1267–1273, 2008.
  * [43] K. Zhou, J. C. Doyle, and K. Glover.  Robust and Optimal Control.  Prentice-Hall, Inc., USA, 1996.
  * [44] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu.  The Finite Element Method: Its Basis and Fundamentals, Sixth Edition.  Butterworth-Heinemann, 6 edition, May 2005.