Modeling and state-space identification of deformable mirrors

: To develop high-performance controllers for adaptive optics (AO) systems, it is essential to ﬁrst derive suﬃciently accurate state-space models of deformable mirrors (DMs). However, it is often challenging to develop realistic large-scale ﬁnite element (FE) state-space models that take into account the system damping, actuator dynamics, boundary conditions, and multi-physics phenomena aﬀecting the system dynamics. Furthermore, it is challenging to establish a modeling framework capable of the automated and quick derivation of state-space models for diﬀerent actuator conﬁgurations and system geometries. On the other hand, for accurate model-based control and system monitoring, it is often necessary to estimate state-space models from the experimental data. However, this is a challenging problem since the DM dynamics is inherently inﬁnite-dimensional and it is characterized by a large number of eigenmodes and eigenfrequencies. In this paper, we provide modeling and estimation frameworks that address these challenges. We develop an FE state-space model of a faceplate DM that incorporates damping and actuator dynamics. We investigate the frequency and time domain responses for diﬀerent model parameters. The state-space modeling process is completely automated using the LiveLink for MATLAB toolbox that is incorporated into the COMSOL Multiphysics software package. The developed state-space model is used to generate the estimation data. This data, together with a subspace identiﬁcation algorithm, is used to estimate reduced-order DM models. We address the model-order selection and model validation problems. The results of this paper provide essential modeling and estimation tools to broad AO and mechatronics scientiﬁc communities. The developed Python, MATLAB, and COMSOL Multiphysics codes are available online


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][5][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][10][11][12][13][14][15][16][17][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][19][20][21][22][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][27][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.symbol explanation symbol explanation A ∈ R s×s , B ∈ R s×m identified system matrices s identified system state order

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.

Faceplate model and analysis
We consider a circular mirror with a diameter of [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.

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]: where is the control input vector containing actuator forces (that will be defined in the sequel), z and 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][38][39][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 = αM 1 + βM 3 , where α, β ∈ R are the parameters that are determined by solving the following system of equations (for more details see Chapter 9 in [45]): where ω i , i = 1, 2, are the eigenfrequencies and ζ 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 where δ is the logarithmic decrement, and the variables η 1 , η 2 , η 3 , and η 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 η 2 = 12[µm] [17,46], we obtain ζ = 0.03.We solve the system in Eq. ( 2 .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. 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][47][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 γ m z a + γ d z a + γ s z a = F a , where γ m , γ d , and γ 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 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.
Next, we combine damped faceplate model with actuator models.We consider two actuator types, referred to as Actuator  [10,15], and this can be incorporated in our modeling approach.Since an estimate of the faceplate stiffness is 5 • 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) • 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 .

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 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 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: where y ∈ R r is the output vector that consists of the first r Zernike coefficients (that are functions of time), and the matrix C z ∈ R r×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]: where w ∈ R n , n = 2n 1 , and the matrices E ∈ R n×n , R ∈ R n×n , T ∈ R n×m , and Q ∈ R r×n are defined as follows: 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, . ..).Then, the discretized model has the following form: where w k ∈ R n is the approximation of the vector w, 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.

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: where x k ∈ R s is the state vector that is a (lower-dimensional and projected) estimate of the state w k of the data generating model in Eq. ( 6), e k ∈ R r is the innovation vector, A ∈ R s×s , B ∈ R s×m , and K ∈ R s×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 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 {y 0 , u 0 , y 1 , u 1 , . . ., y N , u N }.
We use a Subspace Identification (SI) algorithm presented in [58] and briefly summarized in [57] to estimate the model.

Subspace identification algorithm
To state the SI algorithm, we introduce the following notation.Let d k ∈ R g be an arbitrary g-dimensional vector.For two arbitrary positive integers c and q, c ≤ q , the notation d c,q defines the (q − c + 1)g dimensional vector The vector d c,q is called the lifted vector of the vector d k .Similarly, for positive integers c, q, l, c<q<l, the notation D (l) c,q denotes a ((q − c + 1)g) × (l + 1) matrix 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, . . ., q of the matrix D. Similarly, the notation D(:, c : q) denotes a matrix composed of the columns c, c + 1, . . ., q of the matrix D.
The model in Eq. ( 7) can be written as follows where , and v k ∈ 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.
where the matrices Y (l) p,p and V (l) 0,p−1 are constructed using the vectors y k and v k , respectively, and the notation in Eqs. ( 8) and ( 9).The solution is given by (2) Choose the parameter f = p and form the matrix ∆ ∈ R rf ×p(r+m) defined as follows where 0 i×j is an i × j matrix of zeros.Compute the Singular Value Decomposition (SVD) [26] of the matrix product ∆ • V (l) (3) Choose the state order s (see the comments below the algorithm description) and compute the estimated state sequence matrix X(l) p,p ∈ R s×(l+1) as follows X(l) p,p = Σ(1 : s, 1 : s) 1/2 V(1 : s, :).
(4) From the matrix X(l) p,p construct the matrices X(l−1) p+1,p+1 and X(l−1) p,p (to define these vectors we use the notation introduced in Eq. ( 9)).Using the vector v k construct the matrix V (l−1) p,p .Construct the matrix S 1 and compute the matrix S 2 The system matrices are estimated as follows where the notation [ •] 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 M in Eq. ( 12) a large number of times for increasing values of p = 1, 2, 3, . . ., p max .Due to the model large-scale nature, and a relatively high computational complexity of computing M, this approach is not feasible for large values of p max .Consequently, in this manuscript, we compute a sequence of the solutions, where p max is a relatively small number (in our case p 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, . . ., p max , and for every value of the state order s that is in the interval s min , s min + 1, . . ., s max , where s min and s 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 Σ 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 {u k } are white noise sequences.Similarly, the initial states are normal random vectors.In both cases, the simulated outputs y k are corrupted by white measurement noise.In this way, we generate two input-output sequences {u k , 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 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 ε = y 0,N − ŷ0,N 2 / y 0,N 2 , where y 0,N is the lifted output sequence of the validation data-set and ŷ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].

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 • 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 ∆ • V (l) 0,p−1 in Eq. ( 14). Figure 7(b) shows the validation error ε (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 Â 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.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.
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.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.

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.

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 .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 .

Fig. 3 .
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 .

Fig. 4 .
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].

Fig. 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).

Fig. 6 .
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•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).

( 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 M ∈ R r×(r+m)p of the Markov parameters[58] by solving the following least-squares problem min M Y(l)  p,p − MV(l)

Fig. 7 .
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.

Fig. 8 .
Fig.8.Eigenvalues of the matrices Â (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.

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

Fig. 10 .
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 ith and jth entry of the residual vector.