Modeling and state-space identification of deformable mirrors

Aleksandar Haber; Michel Verhaegen

doi:10.1364/OE.382880

1. Introduction

Adaptive Optics (AO) is a widely used technology for improving the performance of optical systems, such as telescopes, microscopes, free-space communication links, optical lithography machines, high-power lasers, etc. [1]. The main components of AO control systems are WaveFront Sensors (WFSs), controllers and Deformable Mirrors (DMs). A WFS is used to measure the wavefront aberrations and on the basis of these measurements, a controller computes control signals for DM actuators. The actuators deform the reflective surface of a DM in such a way that the optical aberrations in the reflected wavefront are attenuated. A relatively recent overview of the control aspects of AO systems is presented in [2]. Due to the need for higher spatial and temporal correction performances of AO systems, the number and spatial resolution of DM actuators are constantly increasing [1,3]. Thus, DMs in AO systems for extremely large telescopes have hundreds or even thousands of spatially distributed actuators [3]. To maximize the correction bandwidth of AO systems it is essential to develop sufficiently accurate models of DMs, WFSs, and possibly of other processes that influence the AO dynamics, such as atmospheric turbulence in the case of telescopes [2,4–6].

There are various types of DM architectures and design concepts used in AO systems. Depending on how the reflective surface is designed, we can distinguish between segmented, faceplate, and membrane DMs [1,7,8]. Due to their widespread usage, in this manuscript, we consider DMs with continuous faceplates. The results can easily be extended to other mirror concepts. A typical continuous faceplate DMs is composed of a thin faceplate that is deformed by spatially distributed actuators placed on its backside.

Modeling and control problems for large-scale DMs have been considered in [5,9–18]. Distributed control and estimation problems or problems of efficient implementation of control and estimation algorithms for large-scale DMs and AO systems have been considered in [4,18–23]. Before these methods can be used, it is necessary to develop sufficiently accurate DM models. Furthermore, it is preferable that these models are in a state-space form that can be used for the design of advanced model-based control algorithms. Analytical mirror models that are based on the biharmonic plate equation and its modifications have been developed in [7,10,11,24,25]. However, such first-principle approaches assume that all the geometrical and physical parameters, as well as boundary conditions and multiphysics phenomena that describe the system dynamics, are known a priori. However, all critical system parameters might not be known accurately, or they vary with operational conditions. In these cases, the system model needs to be estimated (identified) from the experimental data.

The problem of estimating DM model parameters, such as damping, has been considered in [11]. However, this approach does not address the problem of selecting the model order. Subspace Identification (SI) algorithms are powerful estimation methods that can estimate model orders and state-space models by performing a few numerical steps, without the need to employ complex optimization methods [26–28]. SI algorithms have been used in [29,30] to estimate state-space models of thermally actuated DMs. An SI algorithm has also been used in [31] to estimate a state-space model of a bimorph mirror. In [32], an SI algorithm has been used to estimate a model of a DM with a small number of piezoelectric actuators. On the other hand, in [9] tensor-based identification methods have been used to estimate AO models. However, to the best of our knowledge, the potential of using SI algorithms for the estimation of large-scale faceplate DM models has not been systematically investigated. The problems of estimating model orders and state-space descriptions of faceplate DMs are challenging identification problems mainly because the dynamics of these mirror types are inherently infinite-dimensional and systems consist of large numbers of inputs and outputs. Furthermore, additional identification challenges are imposed by the fact that the DM dynamics is characterized by a large number of eigenmodes and eigenfrequencies (some of which are lightly damped).

Apart from the model estimation problems, when designing controllers and estimators we are often faced with another important modeling problem that is often overlooked in the literature. Namely, when all the system parameters are known, there is a need to automatically and quickly derive Finite Element (FE) state-space models [33] of DMs for arbitrary system geometries, actuator configurations, and boundary conditions.

In this manuscript, we develop a large-scale state-space FE model of a faceplate DM and its actuators. We investigate both the frequency and time domain system responses for different model parameters, such as damping and stiffness. We demonstrate that by using the LiveLink for MATLAB toolbox, that is incorporated into the COMSOL Multiphysics software package, it is possible to completely automatize the process of deriving large-scale FE state-space models. This modeling framework can easily incorporate arbitrary geometries, various actuator configurations, model damping, and actuator dynamics. Using the developed model as a data generating model, we obtain sets of input-output data. Using these data sets and the SI algorithm, we estimate a reduced-order DM model. We investigate how well the SI algorithm can be used to estimate the dominant eigenvalues of the DM model and the model order. Furthermore, we investigate and address a model overfitting problem, and perform the model validation on the basis of a correlation-based residual test. The results of this paper provide useful modeling and estimation tools to broad AO and mechatronics scientific communities. The codes used to generate the results in this paper are available online [34].

This paper is organized as follows. In Section 2 we present the DM and actuator models and perform frequency domain and time domain analyses. In Section 3, we define a measurement equation and a state-space model that is simulated to obtain the identification data. In Section 4, we summarize the used SI algorithm and present the identification results. Finally, in Section 5, we present conclusions. The summary of the symbols used in this paper is given in Table 1.

Table 1. Glossary of the main symbols used in the manuscript.

View Table

2. Modeling and analysis

In this section, we present the mirror model and perform eigenmode, eigenfrequency, and time-domain studies. Furthermore, to improve the model realisticity, we also model the system damping and actuator dynamics. The results presented in this section enable us to draw important conclusions on the system dynamics and to properly choose system identification parameters.

2.1 Faceplate model and analysis

We consider a circular mirror with a diameter of $1[m]$. The mirror thickness is $0.002[m]$. The mirror is made of ZERODUR, with the density, Young’s modulus, and Poisson’s ratio equal to $2530 [kg/m^{3}]$, $9.03 \cdot 10^{10} [Pa]$, and $0.24$, respectively. We assume a rectangular actuator configuration. The mirror deformation is modeled using the Kirchhoff-Love theory of plates [35]. There are three degrees of freedom: the out-of-plane displacement and two displacements of shell normals. We assume that the outer mirror edges are clamped. We assume that there are no initial deformations and that initial velocities of the material points are zero. A discretized mirror model is obtained using the Finite Element (FE) method implemented in the COMSOL Multiphysics software. Figure 1(a) shows the mirror geometry and the FE mesh. To quantify the faceplate dynamics, we first perform the frequency response analysis. We assume that a harmonic point load of the amplitude of $1[N]$ is acting at the point $P_{1}=(-0.15,-0.15)$. We calculate the magnitude responses at the points $P_{1}$ and $P_{2}=(0,0)$. The magnitude responses are shown in Figs. 1(b) and (c). While interpreting the results, it should be kept in mind that COMSOL Multiphysics artificially introduces structural damping in order to approximate the natural frequencies. To test the solution convergence and to obtain a more accurate results mesh refinement technique should be used. To further characterize the mirror model, we perform the modal analysis. The first few eigenfrequencies (natural frequencies) with the corresponding eigenmodes are shown in Fig. 2.

Fig. 1. (a) The FE mesh with the mesh quality. The magnitude response observed at (b) the point $P_{1}$ and (c) the point $P_{2}$. For the calculation of the magnitude responses, we assume that the harmonic force of the amplitude of $1$ [N] is acting at the point $P_{1}$.

Download Full Size | PDF

Fig. 2. The first few eigenfrequencies and the corresponding eigenmodes of the faceplate.

Download Full Size | PDF

2.2 Damping and actuator modeling

To obtain a more realistic model, we incorporate the system damping and actuator dynamics in the DM model. In addition, we analyze their effects on the system dynamics. The discretized FE model has the following form [33]:

(1)$$M_{1}\ddot{\mathbf{z}}+M_{2}\dot{\mathbf{z}}+M_{3}\mathbf{z}=B_{1}\mathbf{u},$$

where $M_{i}\in \mathbb {R}^{n_{1} \times n_{1}}$, $i=1,2,3$, are the FE matrices, $B_{1}\in \mathbb {R}^{n_{1}\times m}$ is the control input matrix, $\mathbf {z}=\mathbf {z}(t) \in \mathbb {R}^{n_{1}}$ is the discretized displacement vector, $t$ is time, $\mathbf {u}\in \mathbb {R}^{m}$ is the control input vector containing actuator forces (that will be defined in the sequel), $\dot {\mathbf {z}}$ and $\ddot {\mathbf {z}}$ denote the first and second time derivatives. The mass $M_{1}$, damping $M_{2}$, stiffness $M_{3}$, and $B_{1}$ model matrices are exported from the COMSOL Multiphysics environment to the MATLAB environment by using the LiveLink interface. Furthermore, the complete modeling procedure is parametrized using the LiveLink interface. An explanation of the LiveLink codes used to generate the models in this manuscript is presented in [36]. The codes used in this manuscript are available online [34].

The system damping depends on the DM concept and the actuator type and it positively influences the system stability. For DM designs with gaps between the faceplate and the back plate, the damping is provided by air trapped in gaps [37–40], and the surrounding air. We refer to this damping as the global damping. For plates oscillating in an air environment, it has been experimentally shown that damping depends on the amplitude and is a nonlinear function of the surface area [41]. Furthermore, the damping also originates from the faceplate internal structural damping and mechanical connections used as constraints [41]. In [41] it has been shown that for systems with a relatively large ratio of area to mass, the air damping is an order of magnitude greater than other damping sources. The local actuator damping also plays an important role in overall system damping. We refer to this damping as the local damping. The local actuator damping can be passive or active. For example, passive damping can be provided by fluid trapped between moving parts of voice-coil actuators [39,42]. The active damping can be achieved by using a local feedback update from a local position sensor [43,44].

Rayleigh damping method. The global actuator damping and in some cases the local actuator damping (for actuators with stiffness smaller than the faceplate stiffness) can be incorporated into the model using the Rayleigh method [45]. This method is based on modal system analysis and it selects the damping matrix as follows [45]: $M_{2}=\alpha M_{1}+\beta M_{3}$, where $\alpha , \beta \in \mathbb {R}$ are the parameters that are determined by solving the following system of equations (for more details see Chapter 9 in [45]):

(2)$$\alpha +\beta \omega_{1}^{2}=2\omega_{1}\zeta_{1}, \;\; \alpha +\beta \omega_{2}^{2}=2\omega_{2}\zeta_{2},$$

where $\omega _{i}$, $i=1,2$, are the eigenfrequencies and $\zeta _{i}$, $i=1,2$, are the damping ratios which are considered as constants. Since the air damping is more dominant than the structural damping, we estimate the air damping using the approach experimentally verified in [41]. The damping ratio is determined by $\zeta = 1/ \sqrt {1+(2\pi / \delta )^{2}}, \;\; \delta =22 \eta _{1} \eta _{2} \eta _{3}^{4/3}/ \eta _{4}$, where $\delta$ is the logarithmic decrement, and the variables $\eta _{1}$, $\eta _{2}$, $\eta _{3}$, and $\eta _{4}$ are air density, amplitude of oscillation, faceplate surface area, and faceplate mass, respectively, expressed in U.S. customary units. For the atmospheric pressure and for $\eta_{2}=12[\mu m]$ [17,46], we obtain $\zeta = 0.03$. We solve the system in Eq. (2) for $\zeta =\zeta _{1}=\zeta _{2}=0.03$, $\omega _{1}=190 [rad/s]$, and $\omega _{2}=314 [rad /s]$. The solution is $\alpha =7.1024$ and $\beta =1.19\cdot 10^{-4}$. Next, we compute a frequency response of the Rayleigh damped model. We insert a harmonic point load of $1[N]$ at the point $P_{1}$, and observe the responses at the points $P_{1}$ and $P_{2}$. Figure 3 shows the magnitude responses of the Rayleigh damped model. For comparison, we show the magnitude responses of the faceplate without damping. Due to the introduced damping, the eigenmodes and eigenfrequencies are complex and we omit them.

Fig. 3. The effect of introduced air damping by using the Rayleigh method. A harmonic force is applied at the point $P_{1}$. The response observed at (a) $P_{1}$ and at (b) $P_{2}$.

Download Full Size | PDF

Damping through actuator dynamics and actuator model. The local damping and actuators can be modeled as mass-spring-damper systems. This approach is widely used in literature for modeling linearized electromechanical and piezo-stack array actuators, see for example [15,24,42,46–48]. Mass-spring-damper systems are good first approximations of the actuator dynamics and can model various actuator configurations and structures. The LiveLink codes used in the paper can easily be modified to take into account other, possibly nonlinear actuator models. We incorporate models of electromechanical actuators (voice coils) in our model. The presented modeling approach can be applied to actuators with moving rods attached to the faceplate [48,49] or to actuators with permanent magnets attached to the faceplate [42,50]. The actuator stiffness can be tuned during the design phase, and is provided by an elastic suspension [49] or using some other mechanisms. We assume that the damping is provided by internal viscous damping (passive) [39,42]. Also, the approach can easily be modified to take into account the active damping established by feedback from local position sensors [43,48]. The actuator model is $\gamma _{m}\ddot {z}_{a}+\gamma _{d}\dot {z}_{a}+\gamma _{s}z_{a}=F_{a}$, where $\gamma _{m}$, $\gamma _{d}$, and $\gamma _{s}$, are the actuator mass, damping, and spring constant, $z_{a}$ is the displacement of the actuator, and $F_{a}$ is the electromagnetic force developed by the voice coil. We assume that $z_{a}$ is equal to the local out-of-plane faceplate deformation. In our simulations, the local actuator dynamics is coupled with the plate model with global damping introduced by the Rayleigh method. The actuator forces $F_{a}$ are entries of the vector $\mathbf {u}$ in Eq. (1).

On the basis of the actuator stiffness, we can distinguish between the force and displacement actuators [15,24,51]. The stiffness of the force actuators is generally smaller than the stiffness of the plate. The stiffness of the displacement actuators (piezoelectric actuators for example) is equal to or larger than the faceplate stiffness [15,24]. Voice coil actuators can be designed such that they are force (more often) and displacement actuators (in some cases). On the other hand, piezoelectric actuators are often designed as displacement actuators. For clarity, we first investigate the influence of the actuator stiffness on the eigenfrequencies and eigenmodes. We assume an actuator spacing of $0.1[m]$ and actuator mass of $0.1[kg]$. This produces a total of $69$ actuators inside of the circular mirror aperture. Figure 4 shows the first and second eigenfrequencies and eigenmodes for different values of the actuator stiffness constant. As expected, the global frequencies increase with increased stiffness and the eigenmodes transition from the global modes to the modes with localized structures. This is in accordance with the results reported in [49]. Furthermore, as the stiffness increases, we observed that the modes tend to cluster around localized frequencies.

Fig. 4. The eigenfrequencies and eigenmodes of the DM for several values of the actuator stiffness. The upper (lower) row shows the first (second) eigenfrequencies and eigenmodes. The results are generated for the actuator spacing of $0.1[m]$ and actuator mass of $0.1[kg]$.

Download Full Size | PDF

Next, we combine damped faceplate model with actuator models. We consider two actuator types, referred to as Actuator 1 and Actuator 2 in the text. The parameters of Actuator 1 are $\gamma _{m}=0.1 [kg]$, $\gamma _{d}=40 [Ns/m]$, and $\gamma _{s}=40\cdot 10^3 [N/m]$. The damping ratio of this actuator is $\zeta _{A1}= 0.39$ and the bandwidth is $f_{A1}=138 [Hz]$. The parameters of Actuator 2 are $\gamma _{m}=0.1 [kg]$, $\gamma _{d}=50 [Ns/m]$, and $\gamma _{s}=100 \cdot 10^3 [N/m]$. The damping ratio of this actuator is $\zeta _{A2}= 0.25$ and the bandwidth is $f_{A2}=236.2670 [Hz]$. The actuator bandwidth can be increased by local feedback [10,15], and this can be incorporated in our modeling approach. Since an estimate of the faceplate stiffness is $5\cdot 10^{5} [N/m]$ [52], Actuator 1 (Actuator 2) is a force (displacement) actuator. Depending on the actuator structure and purpose, as well as on the desired control bandwidth, a variety of actuators with different dynamical parameters have been described in the literature [10,15,43,49]. The actuator dynamical parameters used in this manuscript are selected such that they are within the range of values reported in the literature for voice-coil actuators. The damping values are chosen such that they are within the range of damping values reported for $(10-100) \cdot 10^{-6} [m]$ gaps, see Fig. 8 in [42]. Actuator 1 (Actuator 2) stiffness is selected such that it is a force (displacement) actuator. The mass is chosen such that it corresponds to the mass of a cylindrical object with the diameter of $0.03[m]$, the height of $0.02 [m]$, that is made of iron. Furthermore, we have chosen the mass and stiffness parameters such that the actuator open-loop frequency response is similar to the one shown in Figs. 3 and 4 in [15] (where we have neglected the suction cup effect). The approach presented in this paper is general and applicable to arbitrary values of system parameters.

Figure 5 shows the global steady-state mirror deformation and step responses, respectively, at different points for the two actuators. In both cases, the force of $3[N]$ is acting at the point $P_{5}=(-0.2,-0.2)$. Figure 6 shows the magnitude responses for Actuator 1 model when a harmonic force is acting at the point $P_{5}$.

Fig. 5. (a,c) The global deformation and (b,d) step responses. The results shown in (a,b) are generated for Actuator 1 model and in (c,d) are generated for Actuator 2 model. In both cases, the force of $3[N]$ is acting at the point $P_{5}=(-0.2,-0.2)$.

Download Full Size | PDF

Fig. 6. The magnitude responses of the system models. The red line (”Rayleigh+stiffness”) corresponds to a system model obtained by coupling the damped faceplate model (Rayleigh damping) with actuators spaced at $0.1[m]$, actuator mass of $m=0.1[kg]$, stiffness of $4\cdot 10^{3}[N/m]$, and zero damping. The black line (”Rayleigh+full actuator model”) corresponds to a system model obtained by coupling the damped faceplate model (Rayleigh damping) with Actuator 1 model. The harmonic force of $1 [N]$ is acting at $P_{5}=(-0.2,-0.2)$ (forces of other actuators are zero). The responses are observed at (a) $P_{5}$ and (b) $P_{2}=(0,0)$.

Download Full Size | PDF

3. Measurement equation and data generating model

In this section, we introduce a measurement equation and a data-generating model. By simulating the data-generating model, we obtain the estimation data set that is used in Section 4 to identify the model.

The displacement vector $\mathbf {z}$ in Eq. (1) is usually not completely observable. This important modeling fact is often overlooked in the literature. First of all, among all the components of the vector $\mathbf {z}$, only the out-of-plane components can be directly observed. Secondly, only a deformation of a portion of the mirror top surface can be directly observed. That is, the measurement radius is smaller than the total DM radius. In our case, the measurement radius is $0.4 [m]$ (the total DM radius is $0.5 [m]$). This deformation can be directly or indirectly observed using various devices and techniques, such as for example, using WFSs or interferometers [4]. We express these displacements (belonging to the measurement area) in the Zernike basis [53]. Our approach can easily be modified such that the deformation is expressed using some other basis functions. Consequently, our measurement equation has the following form:

(3)$$\mathbf{y}=C_{z}\mathbf{z},$$

where $\mathbf {y}\in \mathbb {R}^{r}$ is the output vector that consists of the first $r$ Zernike coefficients (that are functions of time), and the matrix $C_{z}\in \mathbb {R}^{r\times n_{1}}$ is a matrix that maps the surface deformations into the Zernike coefficients. The model in Eq. (1) and the measurement equation in Eq. (3) can be written in the descriptor state-space form [54]:

(4)$$E\dot{\mathbf{w}}=R\mathbf{w}+T\mathbf{u}, \;\; \mathbf{y} =Q\mathbf{w},$$

where $\mathbf {w} \in \mathbb {R}^{n}$, $n=2n_{1}$, and the matrices $E\in \mathbb {R}^{n\times n }$, $R\in \mathbb {R}^{n\times n}$, $T\in \mathbb {R}^{n\times m}$, and $Q\in \mathbb {R}^{r \times n}$ are defined as follows:

(5)$$\begin{aligned}\mathbf{w}=\begin{bmatrix}\mathbf{z} \\ \dot{\mathbf{z}} \end{bmatrix},\;\; E=\begin{bmatrix} I & 0 \\ 0 & M_{1} \end{bmatrix}, \;\; R=\begin{bmatrix} 0 & I \\ -M_{3} & -M_{2} \end{bmatrix}, \;\; T=\begin{bmatrix} 0 \\ B_{1} \end{bmatrix}, \;\; Q=\begin{bmatrix}C_{z} & 0 \end{bmatrix}. \end{aligned}$$

In practice, the measurement data is obtained at the equidistant discrete-time instants. Consequently, it is more convenient to express Eq. (4) in the discrete-time domain. Due to its stability and robustness, we use the backward Euler discretization method to obtain a discrete-time model [55,56]. Let $h$ and $k$ denote a discretization constant (period) and a discrete-time instant, respectively ($t=hk$, where $k=0,1,2,\ldots$). Then, the discretized model has the following form:

(6)$$\mathbf{w}_{k}=A_{1} \mathbf{w}_{k-1} +B_{1} \mathbf{u}_{k}, \;\; \mathbf{y}_{k}= Q\mathbf{w}_{k},$$

where $\mathbf {w}_{k} \in \mathbb {R}^{n}$ is the approximation of the vector $\mathbf {w}$, $A_{1}=\big (E-hR \big )^{-1}E$, and $B_{1}=h\big (E-hR \big )^{-1}T$.

The model in Eq. (6) represents the data generating model. Since we do not have access to an experimental setup, by simulating the data generating model for chosen initial conditions and input sequences, we obtain the identification input-output data. Then, from this input-output data, we want to estimate a state-space model that reproduces as accurately as possible the input-output behavior of the data generating model. In Section 4, we solve this problem and present the identification results.

4. Subspace identification

In this section, we first formulate the target model that we want to estimate and then summarize the identification algorithm. At the end of the section, we present the identification results.

Our main goal is to estimate a Kalman innovation state-space model [26,57] that approximates the input-output behavior of Eq. (6). The Kalman innovation state-space model has the following form:

(7)$$\mathbf{x}_{k+1} =A\mathbf{x}_{k}+B\mathbf{u}_{k}+K\mathbf{e}_{k}, \;\; \mathbf{y}_{k} =C\mathbf{x}_{k}+\mathbf{e}_{k},$$

where $\mathbf {x}_{k}\in \mathbb {R}^{s}$ is the state vector that is a (lower-dimensional and projected) estimate of the state $\mathbf {w}_{k}$ of the data generating model in Eq. (6), $\mathbf {e}_{k}\in \mathbb {R}^{r}$ is the innovation vector, $A\in \mathbb {R}^{s\times s}$, $B \in \mathbb {R}^{s\times m}$, and $K \in \mathbb {R}^{s \times r}$ are the system matrices. The parameter $s$ is an estimate of the state order $n$ of the data generating system in Eq. (6). Equation (7) is generic in the sense that through the innovation vector $\mathbf {e}_{k}$ and the matrix $K$ it implicitly models the measurement noise and the process disturbances [26]. In this way, we can incorporate the measurement noise and process disturbances (or stochastic unmodeled dynamics) that always affect the DM dynamics.

Formally speaking, the identification problem is to estimate the model order $s$, and the system matrices $A$, $B$, $K$, and $C$ of the state-space model in Eq. (7), from the set of the observed input-output data $\{ \mathbf {y}_{0},\mathbf {u}_{0},\mathbf {y}_{1},\mathbf {u}_{1},\ldots , \mathbf {y}_{N},\mathbf {u}_{N} \}$.

We use a Subspace Identification (SI) algorithm presented in [58] and briefly summarized in [57] to estimate the model.

4.1 Subspace identification algorithm

To state the SI algorithm, we introduce the following notation. Let $\mathbf {d}_{k}\in \mathbb {R}^{g}$ be an arbitrary $g$-dimensional vector. For two arbitrary positive integers $c$ and $q$, $c \le q$ , the notation $\mathbf {d}_{c,q}$ defines the $(q-c+1)g$ dimensional vector

(8)$$\mathbf{d}_{c,q}=\begin{bmatrix}\mathbf{d}_{c}^{T} & \mathbf{d}_{c+1}^{T} & \ldots & \mathbf{d}_{q}^{T} \end{bmatrix}^{T}.$$

The vector $\mathbf {d}_{c,q}$ is called the lifted vector of the vector $\mathbf {d_{k}}$. Similarly, for positive integers $c,q,l$, $c<q<l$, the notation $D_{c,q}^{(l)}$ denotes a $((q-c+1)g)\times (l+1)$ matrix

(9)$$D_{c,q}^{(l)}=\begin{bmatrix}\mathbf{d}_{c,q} & \mathbf{d}_{c+1,q+1} & \ldots & \mathbf{d}_{c+l,q+l} \end{bmatrix}.$$

Let $D$ be an arbitrary matrix, and let $c$ and $q$, $c<q$, be two integers. The notation $D(c:q,:)$ denotes a new matrix composed of the rows ${c,c+1,\ldots , q}$ of the matrix $D$. Similarly, the notation $D(:,c:q)$ denotes a matrix composed of the columns ${c,c+1,\ldots , q}$ of the matrix $D$.

The model in Eq. (7) can be written as follows

(10)$$\mathbf{x}_{k+1} =\tilde{A}\mathbf{x}_{k}+\tilde{B}\mathbf{v}_{k}, \;\; \mathbf{y}_{k} =C\mathbf{x}_{k}+\mathbf{e}_{k},$$

where $\tilde {A}=A-KC$, $\tilde {B}=\begin {bmatrix}B & K \end {bmatrix}$, $\mathbf {v}_{k}=\begin {bmatrix} \mathbf {u}_{k}^{T} & \mathbf {y}_{k}^{T} \end {bmatrix}^{T}$, and $\mathbf {v}_{k}\in \mathbb {R}^{r+m}$.

To formulate the SI algorithm, we need to introduce two parameters that are selected by users. They are the past window $p$ and the future window $f$. The identification algorithm consists of the following steps.

(1) Choose the parameter $p$ (see the comments below the algorithm description). Let the parameter $l$ be selected as follows $l= N-p-1$. Estimate the matrix $\mathcal {M}\in \mathbb {R}^{r \times (r+m)p}$ of the Markov parameters [58] by solving the following least-squares problem $(11)$$\min_{\mathcal{M}} \left\|Y_{p,p}^{(l)} -\mathcal{M} V_{0,p-1}^{(l)} \right\|_{F}^{2},$$$ where the matrices $Y_{p,p}^{(l)}$ and $V_{0,p-1}^{(l)}$ are constructed using the vectors $\mathbf {y}_{k}$ and $\mathbf {v}_{k}$, respectively, and the notation in Eqs. (8) and (9). The solution is given by $(12)$$\hat{\mathcal{M}}=Y_{p,p}^{(l)}\big(V_{0,p-1}^{(l)} \big)^{T}\big( V_{0,p-1}^{(l)} \big( V_{0,p-1}^{(l)} \big)^{T} \big)^{-1}.$$$
(2) Choose the parameter $f=p$ and form the matrix $\hat {\Delta } \in \mathbb {R}^{rf \times p(r+m)}$ defined as follows $(13)$$\hat{\Delta}=\begin{bmatrix} \hat{\mathcal{M}} \\ \begin{bmatrix} 0_{r\times (r+m)} & \hat{\mathcal{M}}\big(:,1:(p-1)(r+m)\big) \end{bmatrix} \\ \begin{bmatrix} 0_{r\times 2(r+m)} & \hat{\mathcal{M}}\big(:,1:(p-2)(r+m) \big) \end{bmatrix} \\ \vdots \\ \begin{bmatrix} 0_{r\times (f-1)(r+m)} & \hat{\mathcal{M}}\big(:,1:(p-f+1)(r+m) \big) \end{bmatrix} \end{bmatrix}.$$$ where $0_{i\times j}$ is an $i\times j$ matrix of zeros. Compute the Singular Value Decomposition (SVD) [26] of the matrix product $\hat {\Delta } \cdot V_{0,p-1}^{(l)}$ $(14)$$\hat{\Delta} \cdot V_{0,p-1}^{(l)}= \mathcal{U}\Sigma \mathcal{V}^{T},$$$
(3) Choose the state order $s$ (see the comments below the algorithm description) and compute the estimated state sequence matrix $\hat {X}_{p,p}^{(l)}\in \mathbb {R}^{s \times (l+1)}$ as follows $\hat {X}_{p,p}^{(l)}=\Sigma (1:s,1:s)^{1/2} \mathcal {V}(1:s,:)$.
(4) From the matrix $\hat {X}_{p,p}^{(l)}$ construct the matrices $\hat {X}_{p+1,p+1}^{(l-1)}$ and $\hat {X}_{p,p}^{(l-1)}$ (to define these vectors we use the notation introduced in Eq. (9)). Using the vector $\mathbf {v}_{k}$ construct the matrix $V_{p,p}^{(l-1)}$. Construct the matrix $S_{1}$ and compute the matrix $S_{2}$ $(15)$$\begin{aligned} S_{2}=\hat{X}_{p+1,p+1}^{(l-1)}S_{1}^{T}\big(S_{1}S^{T}_{1} \big)^{-1}, \;\; S_{1}=\begin{bmatrix} \hat{X}_{p,p}^{(l-1)} \\ V_{p,p}^{(l-1)} \end{bmatrix}. \end{aligned}$$$ The system matrices are estimated as follows $(16)$$\begin{aligned} & \hat{\tilde{A}}=S_{2}(:,1:s),\;\; \hat{B}=S_{2}(:,s+1,s+m), \hat{C}= Y_{p,p}^{(l)}\big(\hat{X}_{p,p}^{(l)} \big)^{T}\big( \hat{X}_{p,p}^{(l)} \big( \hat{X}_{p,p}^{(l)} \big)^{T} \big)^{-1},\\ &\hat{K}=S_{2}(:,s+m+1:s+m+r),\;\; \hat{A}=\hat{\tilde{A}}+\hat{K}\hat{C}. \end{aligned}$$$ where the notation $\hat {[\cdot ]}$ denotes estimates of the matrices in Eqs. (7) and (10).

Choice of the parameters $p$ and $s$. The parameter $p$ is the number of past measurements that need to be taken such that the effect of the system unknown initial state on the current system state can be neglected. The more stable the system is, the smaller is the value of $p$. However, we do not know this parameter a priori. Usually, a relatively good estimate of $p$ is obtained by minimizing a measure based on the Akaike Information Criterion (AIC) [57], that roughly speaking, provides a trade-off between the number of model parameters and the goodness of fit. However, to estimate the past window using this method, we need to compute $\hat {\mathcal {M}}$ in Eq. (12) a large number of times for increasing values of $p=1,2,3,\ldots , p_{\textrm {max}}$. Due to the model large-scale nature, and a relatively high computational complexity of computing $\mathcal {M}$, this approach is not feasible for large values of $p_{\textrm {max}}$. Consequently, in this manuscript, we compute a sequence of the solutions, where $p_{\textrm {max}}$ is a relatively small number (in our case $p_{\textrm {max}}=10$). One of the approaches to reduce the computational burden of the SI method is to use the recursive implementation [59], and it will be explored in future publications. For every value of the past window $p=1,2,\ldots , p_{\textrm {max}}$, and for every value of the state order $s$ that is in the interval $s_{\textrm {min}},s_{\textrm {min}}+1,\ldots ,s_{\textrm {max}}$, where $s_{\textrm {min}}$ and $s_{\textrm {max}}$ are minimal and maximal state orders selected by the user, we estimate a model. The best model order $s$, and consequently, the best model among all the models, is the model that produces the smallest validation error or the Variance Accounted For (VAF) that are explained in the sequel. In some cases, the model order $s$ can be selected by visually inspecting the singular values in the matrix $\Sigma$ in Eq. (14). A wide gap in singular values usually indicates a model order [26].

Data Generation and Input Sequences. By simulating the data generating model in Eq. (6) for two independent input sequences and two random initial conditions, we generated two data sets. For both sets, the input sequences $\{ \mathbf {u}_{k}\}$ are white noise sequences. Similarly, the initial states are normal random vectors. In both cases, the simulated outputs ${\mathbf {y}_{k}}$ are corrupted by white measurement noise. In this way, we generate two input-output sequences $\{\mathbf {u}_{k},\mathbf {y}_{k} \}$. The first sequence, referred to as the identification data set, is used to identify the model, and the second sequence, referred to as the validation data set, is used to validate the model using the procedure explained in the sequel.

Model validation. Using the input sequence that is used to generate the validation data set, we simulate the identified model (determined by the parameters $p$ and $s$), to obtain the simulated output sequence. The identified model is simulated in open-loop, that is by setting $\mathbf {e}_{k}=0$ in Eq. (7). To perform the simulations, we need to estimate an initial state of the identified model. This has been done using the validation data set and by using a simple least-squares approach explained in [57]. The model error is defined by $\varepsilon =\left \| \mathbf {y}_{0,N}-\hat {\mathbf {y}}_{0,N} \right \|_{2} / \left \| \mathbf {y}_{0,N} \right \|_{2}$, where $\mathbf {y}_{0,N}$ is the lifted output sequence of the validation data-set and $\hat {\mathbf {y}}_{0,N}$ is the vector composed of the simulated outputs of the identified system (we are using the notation introduced in Eq. (8)). Using the validation data we also compute a VAF measure that additionally quantifies the final model estimation quality. For the definition of VAF see [26].

Once the final model has been estimated, we further analyze its quality using a residual correlation analysis [57]. On the basis of the error between the output sequence of the identified model and the output sequence in the validation data set, we estimate cross-correlation matrices and perform a white-noise hypothesis test. In the ideal case, the error should have the properties of a white-noise sequence. Consequently, we test the error sequence under the white-noise hypothesis. A detailed procedure is explained in [57].

4.2 Identification results

The codes used in this paper are available online [34]. We obtain the identification results assuming Actuator 1 model with actuator spacing of $0.1[m]$ and with the parameters defined in Section 2.2. The model has $m=69$ actuators and $r=32$ Zernike coefficient outputs. The data-generating model order is $n=4878$. The discretization constant is $h=1\cdot 10^{-4}[s]$. This value of the discretization constant is chosen by analyzing the step responses presented in Fig. 5. We want to ensure that we have 15-20 samples during the rise time, but at the same time we do not want to select too many samples that can negatively affect the model identifiability, see [26,57] for more details on the selection of the discretization (sampling time) constant. The identification and validation input sequences are equal to normal white noise sequences with zero mean and unit variance. The initial states for the identification and validation are normal random vectors with zero means and variances of $0.1$. The proposed approach can also be applied to mirrors with larger actuator numbers. We have created such models and they are available online [34].

First, we test the SI method performance when the outputs of the data generating model are not corrupted by the measurement noise. We set $p=8$ and $N=1200$. Figure 7(a) shows the singular values of the matrix $\hat {\Delta }\cdot V_{0,p-1}^{(l)}$ in Eq. (14). Figure 7(b) shows the validation error $\varepsilon$ (introduced in the previous subsection) and VAF for different model orders $s$ . To test how well the SI method can estimate the eigenvalues of the matrix $A_{1}$ of the data-generating model in Eq. (6), we compute the eigenvalues of the matrix $\hat {A}$ of the identified model and compare them with the eigenvalues of the matrix $A_{1}$. The results are shown in Fig. 8. The results presented in Fig. 8 demonstrate that the SI method is able to relatively accurately estimate the dominant eigenvalues.

Fig. 7. (a) Singular values. (b) VAF and errors for different values of the state order - s. The outputs are not corrupted by measurement noise.

Download Full Size | PDF

Fig. 8. Eigenvalues of the matrices $\hat {A}$ (identified model) and $A_{1}$ (real system-data generating model in Eq. (6)) for estimated state orders of $s=10$, $s=50$, $s=150$, and $s=200$, from left to right, respectively. The results are generated for $N=1200$, past and future windows of $8$. The outputs are not corrupted by measurement noise.

Download Full Size | PDF

Next, we test the performance of the SI algorithm when the measurements are corrupted by noise. The approximate value of the Signal to Noise Ratio (SNR) is $20$. The identification and validation data sets consist of $N= 2500$ discrete-time samples. Figure 9 shows the error and VAF for several model orders and past windows. A relatively long computational time of the SI method prevents us from testing the method for larger values of $p$ and for more realizations of the measurement noise. A remedy for this is to use the recursive SI algorithms. We see that by increasing the past window $p$ we can improve the identification results. The overfitting phenomenon probably occurs for large values of $p$.

Fig. 9. (a) The identification error and (b) VAF in the case of the measurement noise. The results are generated for SNR of $20$.

Download Full Size | PDF

Figures 10(a)–(c) show the outputs (Zernike coefficients) of the data generating model and the identified model. The results are generated for $p=10$ and $s=300$. Figures 10(d)–(f) show the correlations between the entries of the residual vectors. The red dashed lines denote $5 \%$ confidence intervals for the white-noise hypothesis test. The percentages of entries in plots (d), (e), and (f) that exceed this limit are $5.8\%$, $6\%$, and $8 \%$. This indicates that there is a small correlation between the residuals, and this can be completely eliminated by increasing $p$, however, at the expense of the increased computational complexity.

Fig. 10. (a-c) The output of the identified (predicted output) and the data generating model (true output), in the case of measurement noise. (d)–(f) Correlations of the entries of the residual vectors (sequence). The notation $(i,j)$ is the correlation value between the $i$th and $j$th entry of the residual vector.

Download Full Size | PDF

Figure 10 further confirms the good performance of the SI algorithm. We can see that the data-generating model can be relatively accurately estimated even in the presence of the measurement noise. The accuracy can be further improved by increasing the past window $p$.

5. Conclusion

In this paper, we have presented state-space modeling and estimation approaches for faceplate Deformable Mirrors (DMs) used in Adaptive Optics (AO) systems. We have shown that by using the Finite Element (FE) method and the LiveLink for COMOSL Multiphysics package, we are able to automatize the derivation of DM state-space models. Furthermore, we have shown how to incorporate damping and actuator dynamics in such a model. This modeling framework can easily handle complex geometries and different actuator configurations. Using this model and the subspace identification methods, we have shown that it is possible to accurately estimate reduced-order state-space models of DMs. The future research directions will be oriented toward decreasing the computational complexity of the subspace identification algorithm and toward experimental verification of the used methods.

Funding

H2020 European Research Council (ERC 339681); City University of New York (PSC-CUNY Award A (61303-00 49), PSC-CUNY Award A (62267-00 50)).

Disclosures

The authors declare no conflicts of interest.

References

1. R. Tyson, Principles of Adaptive Optics (CRC Press, 2015).

2. C. Kulcsár, H.-F. Raynaud, C. Petit, and J.-M. Conan, “Minimum variance prediction and control for adaptive optics,” Automatica 48(9), 1939–1954 (2012). [CrossRef]

3. S. Hippler, “Adaptive optics for extremely large telescopes,” arXiv preprint arXiv:1808.02693 (2018).

4. A. Haber and M. Verhaegen, “Framework to trade optimality for local processing in large-scale wavefront reconstruction problems,” Opt. Lett. 41(22), 5162–5165 (2016). [CrossRef]

5. C. Yu and M. Verhaegen, “Structured modeling and control of adaptive optics systems,” IEEE Trans. Contr. Syst. Technol. 26(2), 664–674 (2018). [CrossRef]

6. M. Manetti, M. Morandini, P. Mantegazza, R. Biasi, D. Gallieni, and A. Riccardi, “Experimental validation of massively actuated deformable adaptive mirror numerical models,” Control. Eng. Pract. 20(8), 783–791 (2012). [CrossRef]

7. C. R. Vogel and Q. Yang, “Modeling, simulation, and open-loop control of a continuous facesheet MEMS deformable mirror,” J. Opt. Soc. Am. A 23(5), 1074–1081 (2006). [CrossRef]

8. S. Kuiper, N. Doelman, J. Human, R. Saathof, W. Klop, and M. Maniscalco, “Advances of TNO’s electromagnetic deformable mirror development,” Proc. SPIE 10706, 1070619 (2018). [CrossRef]

9. B. Sinquin and M. Verhaegen, “Tensor-based predictive control for extremely large-scale single conjugate adaptive optics,” J. Opt. Soc. Am. A 35(9), 1612–1626 (2018). [CrossRef]

10. R. Heimsten, M. Owner-Petersen, T. Ruppel, D. G. MacMynowski, and T. Andersen, “Suppressing low-order eigenmodes with local control for deformable mirrors,” Opt. Eng. 51(2), 026601 (2012). [CrossRef]

11. T. Ruppel, S. Dong, F. Rooms, W. Osten, and O. Sawodny, “Feedforward control of deformable membrane mirrors for adaptive optics,” IEEE Trans. Contr. Syst. Technol. 21(3), 579–589 (2013). [CrossRef]

12. M. Manetti, M. Morandini, P. Mantegazza, R. Biasi, M. Andrighettoni, and D. Gallieni, “VLT DSM, the control system of the largest deformable secondary mirror ever manufactured,” Proc. SPIE 9148, 91484G (2014). [CrossRef]

13. G. Agapito, S. Baldi, G. Battistelli, D. Mari, E. Mosca, and A. Riccardi, “Automatic tuning of the internal position control of an adaptive secondary mirror,” Eur. J. Control 17(3), 273–289 (2011). [CrossRef]

14. M. Manetti, M. Morandini, and P. Mantegazza, “Completely, partially centralized and fully decentralized control schemes of large adaptive mirrors,” J. Vib. Control 22(5), 1288–1305 (2016). [CrossRef]

15. T. Andersen, O. Garpinger, M. Owner-Petersen, F. Bjoorn, R. Svahn, and A. Ardeberg, “Novel concept for large deformable mirrors,” Opt. Eng. 45(7), 073001 (2006). [CrossRef]

16. M. Manetti, M. Morandini, and P. Mantegazza, “Self-tuning shape control of massively actuated adaptive mirrors,” IEEE Trans. Contr. Syst. Technol. 22(3), 838–852 (2014). [CrossRef]

17. R. Hamelinck, R. Ellenbroek, N. Rosielle, M. Steinbuch, M. Verhaegen, and N. Doelman, “Validation of a new adaptive deformable mirror concept,” Proc. SPIE 7015, 70150Q (2008). [CrossRef]

18. D. G. MacMartin, “Control challenges for extremely large telescopes,” Proc. SPIE 5054, 275–286 (2003). [CrossRef]

19. D. G. MacMynowski, R. Heimsten, and T. Andersen, “Distributed force control of deformable mirrors,” Eur. J. Control 17(3), 249–260 (2011). [CrossRef]

20. R. Ellenbroek, M. Verhaegen, N. Doelman, R. Hamelinck, N. Rosielle, and M. Steinbuch, “Distributed control in adaptive optics: deformable mirror and turbulence modeling,” Proc. SPIE 6272, 62723K (2006). [CrossRef]

21. P. Massioni, C. Kulcsár, H.-F. Raynaud, and J.-M. Conan, “Fast computation of an optimal controller for large-scale adaptive optics,” J. Opt. Soc. Am. A 28(11), 2298–2309 (2011). [CrossRef]

22. R. Fraanje, P. Massioni, and M. Verhaegen, “A decomposition approach to distributed control of dynamic deformable mirrors,” Int. J. Optomechatronics 4(3), 269–284 (2010). [CrossRef]

23. D. R. Jenkins, A. Basden, and R. M. Myers, “ELT-scale adaptive optics real-time control with the Intel Xeon Phi many integrated core architecture,” Mon. Not. R. Astron. Soc. 478(3), 3149–3158 (2018). [CrossRef]

24. S. K. Ravensbergen, R. F. H. M. Hamelinck, P. C. J. N. Rosielle, and M. Steinbuch, “Deformable mirrors: design fundamentals for force actuation of continuous facesheets,” Proc. SPIE 7466, 74660G (2009). [CrossRef]

25. M. Böhm and O. Sawodny, “Optimal sensor placement for modal based estimation of deformable mirror shape,” in 2015 IEEE Conference on Control Applications (CCA), (IEEE, 2015), pp. 418–423.

26. M. Verhaegen and V. Verdult, Filtering and System Identification: a Least Squares Approach (Cambridge University, 2007).

27. A. Haber, F. Pecora, M. Chowdhury, and M. Summerville, “Identification of temperature dynamics using subspace and machine learning techniques,” in Proceedings of the ASME 2019 Dynamic Systems and Control Conference, (2019, in press).

28. A. Haber and M. Verhaegen, “Subspace identification of large-scale interconnected systems,” IEEE Trans. Autom. Control 59(10), 2754–2759 (2014). [CrossRef]

29. A. Haber, A. Polo, I. Maj, S. F. Pereira, H. P. Urbach, and M. Verhaegen, “Predictive control of thermally induced wavefront aberrations,” Opt. Express 21(18), 21530–21541 (2013). [CrossRef]

30. A. Haber, A. Polo, S. Ravensbergen, H. P. Urbach, and M. Verhaegen, “Identification of a dynamical model of a thermally actuated deformable mirror,” Opt. Lett. 38(16), 3061–3064 (2013). [CrossRef]

31. A. Chiuso, R. Muradore, and E. Marchetti, “Dynamic calibration of adaptive optics systems: a system identification approach,” IEEE Trans. Contr. Syst. Technol. 18(3), 705–713 (2010). [CrossRef]

32. H. Song, R. Fraanje, G. Schitter, G. Vdovin, and M. Verhaegen, “Controller design for a high-sampling-rate closed-loop adaptive optics system with piezo-driven deformable mirror,” Eur. J. Control 17(3), 290–301 (2011). [CrossRef]

33. O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, vol. 1 (Butterworth-Heinemann, 2005).

34. A. Haber, “Machine learning and system identification for adaptive optics,” (2019) [retrieved 06 November 2019], https://github.com/AleksandarHaber/Machine-Learning-and-System-Identification-for-Adaptive-Optics, .

35. S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells (McGraw-Hill, 1959).

36. A. Haber, F. Pecora, M. Chowdhury, and M. Verhaegen, “Estimation and machine learning of large-scale elastic structures using LiveLink for MATLAB and Python,” in 2019 COMSOL Multiphysics Conference, Boston, (2019, in press).

37. G. Brusa-Zappellini, A. Riccardi, V. Biliotti, C. Del Vecchio, P. Salinari, P. Stefanini, P. Mantegazza, R. Biasi, M. Andrighettoni, C. Franchini, and D. Gallieni, “Adaptive secondary mirror for the 6.5-m conversion of the multiple mirror telescope: first laboratory testing results,” Proc. SPIE 3762, 38–49 (1999). [CrossRef]

38. G. Brusa and C. Del Vecchio, “Design of an adaptive secondary mirror: a global approach,” Appl. Opt. 37(21), 4656–4662 (1998). [CrossRef]

39. G. Brusa, A. Riccardi, M. Accardo, V. Biliotti, M. Carbillet, C. Del Vecchio, S. Esposito, B. Femenia, O. Feeney, L. Fini, and S. Gennari, “From adaptive secondary mirrors to extra-thin extra-large adaptive primary mirrors,” in European Southern Observatory Conference and Workshop Proceedings, vol. 57 (2000), p. 181.

40. M. Manetti, M. Morandini, and P. Mantegazza, “High precision massive shape control of magnetically levitated adaptive mirrors,” Control. Eng. Pract. 18(12), 1386–1398 (2010). [CrossRef]

41. D. G. Stephens and M. A. Scavullo, “Investigation of air damping of cicular and rectangular plates, a cylinder, and a sphere,” Tech. rep., NASA Langley Research Center (1965).

42. G. Brusa-Zappellini, A. Riccardi, S. D. Ragland, S. Esposito, C. Del Vecchio, L. Fini, P. Stefanini, V. Biliotti, P. Ranfagni, P. Salinari, D. Gallieni, R. Biasi, P. Mantegazza, G. Sciocco, G. Noviello, and S. Invernizzi, “Adaptive secondary P30 prototype: laboratory results,” Proc. SPIE 3353, 764–775 (1998). [CrossRef]

43. S. Kuiper, N. Doelman, E. Nieuwkoop, T. Overtoom, T. Russchenberg, M. van Riel, J. Wildschut, M. Baeten, J. Human, H. Spruit, S. Brinkers, and M. Maniscalco, “Electromagnetic deformable mirror development at TNO,” Proc. SPIE 9912, 991204 (2016). [CrossRef]

44. B. Sedghi, M. Dimmler, and M. Müller, “Active damping strategies for control of the E-ELT field stabilization mirror,” Proc. SPIE 8444, 84441Z (2012). [CrossRef]

45. K.-J. Bathe, Finite Element Procedures (Prentice Hall, 1996).

46. J.-C. Sinquin, J.-M. Lurçon, and C. Guillemard, “Deformable mirror technologies for astronomy at CILAS,” Proc. SPIE 7015, 70150O (2008). [CrossRef]

47. A. Preumont, Vibration Control of Active Structures: An Introduction (Springer, 2018).

48. R. Heimsten, D. G. MacMynowski, T. Andersen, and M. Owner-Petersen, “Concept, modeling, and performance prediction of a low-cost, large deformable mirror,” Appl. Opt. 51(5), 515–524 (2012). [CrossRef]

49. R. F. M. M. Hamelinck, “Adaptive Deformable Mirror: Based on Electromagnetic Actuators,” Ph.D. thesis, Eindhoven University of Technology (2010).

50. M. Cayrel, “E-ELT optomechanics: overview,” Proc. SPIE 8444, 84441X (2012). [CrossRef]

51. S. Ravensbergen, “Adaptive optics for extreme ultraviolet lithography : actuator design and validation for deformable mirror concepts,” Ph.D. thesis, Eindhoven University of Technology, (2012).

52. S. K. Ravensbergen, P. C. J. N. Rosielle, and M. Steinbuch, “Deformable mirrors with thermo-mechanical actuators for extreme ultraviolet lithography: Design, realization and validation,” Precis. Eng. 37(2), 353–363 (2013). [CrossRef]

53. R. Tyson, Adaptive Optics Engineering Handbook (CRC Press, 1999).

54. A. Haber and M. Verhaegen, “Sparsity preserving optimal control of discretized PDE systems,” Comput. Method Appl. M. 335, 610–630 (2018). [CrossRef]

55. A. Iserles, A First Course in the Numerical Analysis of Differential Equations (Cambridge University, 2009).

56. A. Haber, F. Molnar, and A. E. Motter, “State observation and sensor selection for nonlinear networks,” IEEE Trans. Control Netw. Syst. 5(2), 694–708 (2018). [CrossRef]

57. A. Haber, “Subspace identification of temperature dynamics,” arXiv preprint arXiv:1908.02379 (2019).

58. I. Houtzager, J.-W. van Wingerden, and M. Verhaegen, “Varmax-based closed-loop subspace model identification,” in Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference, (IEEE, 2009), pp. 3370–3375.

59. I. Houtzager, J.-W. van Wingerden, and M. Verhaegen, “Fast-array recursive closed-loop subspace model identification,” in 15th IFAC Symposium on System Identification, vol. 42 (Elsevier, 2009), pp. 96–101.

symbol	explanation	symbol	explanation
$P_{i}$ , $i = 1, \dots, 5$	faceplate domain points	$z \in R^{n_{1}}$	FE displacement vector
$M_{1}, M_{2}, M_{3} \in R^{n_{1} \times n_{1}}$	mass, damping, stiffness	$u \in R^{m}$	vector of actuator forces
$B_{1} \in R^{n_{1} \times m}$	FE control input matrix	$α, β$	Rayleigh damping parameters
$δ$	logarithmic decrement	$ω, ζ$	frequency and damping ratio
$η_{1}$	air density	$η_{2}$	amplitude of faceplate oscillation
$η_{3}$	faceplate surface area	$η_{4}$	faceplate mass
$z_{a}$	actuator displacement	$F_{a}$	actuator force
$γ_{m}, γ_{d}$	actuator mass, damping	$γ_{s}, f_{A}$	actuator stiffness, bandwidth
$C_{z} \in R^{r \times n_{1}}$	matrix mapping $z$ to $y$	$y \in R^{r}$	Zernike coefficient output vector
$w = [z^{T} {\dot{z}}^{T}]^{T} \in R^{n}$	global state vector	$w_{k}$	discrete-time version of $w$
$E, R \in R^{n \times n}$	state-space matrices ( $w)$	$T \in R^{n \times m}$	input matrix ( $w$ )
$Q \in R^{r \times n}$	matrix mapping $w$ to $y$	$y_{k} \in R^{r}$	discrete-time version of $y$
$A_{1} \in R^{n \times n}$	state-space matrix ( $w_{k}$ )	$B_{1} \in R^{n \times m}$	input matrix ( $w_{k}$ )
$A \in R^{s \times s}, B \in R^{s \times m}$	identified system matrices	$s$	identified system state order
$x_{k} \in R^{s}$	identified state	$e_{k} \in R^{r}$	innovation vector
$K \in R^{n \times r}$	Kalman gain	$p$	past window
$\tilde{A} \in R^{n \times n}$	$\tilde{A} = A - K C$	$N$	number of data samples
$\tilde{B} \in R^{n \times m}$	$\tilde{B} = [B K]$	$v_{k} \in R^{r + m}$	$v_{k} = [u_{k}^{T} y_{k}^{T}]^{T}$
$f$	future window	$l$	identification parameter
$Y_{p, p}^{(l)}, V_{0, p - 1}^{(l)}$	data matrices	$M$	Markov parameter matrix
$\hat{}$	”hat” notation for estimates	$U, Σ, V$	SVD matrices
$\hat{Δ}, S_{1}, S_{2}$	matrices used in SI	$s_{m i n}, s_{m a x}$	minimal and maximal state orders
$ε$	model error	$y_{0, N}, {\hat{y}}_{0, N}$	lifted output, lifted estimated output
$p_{m a x}$	maximal past window	${\hat{X}}_{p, p}^{(l)}$	estimated state matrix
$h$	discretization constant

Modeling and state-space identification of deformable mirrors

Abstract

1. Introduction

2. Modeling and analysis

2.1 Faceplate model and analysis

2.2 Damping and actuator modeling

3. Measurement equation and data generating model

4. Subspace identification

4.1 Subspace identification algorithm

4.2 Identification results

5. Conclusion

Funding

Disclosures

References

Cited By

Figures (10)

Tables (1)

Equations (16)

Optics Express