General approach to precise deformable mirror control

We develop a simple and effective control method for accurate control of deformable mirrors (DMs). For a desired DM surface profile and using batches of observed surface profile data, the proposed method adaptively determines both a DM model (influence matrix) and control actions that produce the desired surface profile with good accuracy. In the first iteration, the developed method estimates a DM influence matrix by solving a multivariable least-squares problem. This matrix is then used to compute the control actions by solving a constrained least-squares problem. Then, the computed actions are randomly perturbed and applied to the DM to generate a new batch of surface profile data. The new data batch is used to estimate a new influence matrix that is then used to re-compute control actions. This procedure is repeated until convergence is achieved. The method is experimentally tested on a Boston Micromachines DM with 140 micro-electronic-mechanical-system actuators. Our experimental results show that the developed control approach can achieve accurate correction despite significant DM nonlinearities. Using only a few control iterations, the developed method is able to produce a surface profile root-mean-square error that varies from 5− 30 [nm] for most of the tested Zernike wave-front modes without using direct feedback control. These results can additionally be improved by using larger data batches and more iterations or by combining the developed approach with feedback control. Finally, as we experimentally demonstrate, the developed method can be used to estimate a DM model that can effectively be used for a single-step open-loop DM control. © 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

In this paper, we present a novel DM control method and we experimentally verify it on a MEMS DM with 140 actuators that is developed by Boston Micromachines [20,21]. MEMS DMs are relatively low-cost, have a fast response, and allow for a relatively large number of actuators within the DM active surface. One of the main characteristics of the used MEMS DM, as well as of DMs that are based on other actuation technologies and mechanisms, is that for a sufficiently large actuator operating range, the deformation response is a nonlinear function of the applied control signals. Intrinsic nonlinearities and nonlinear couplings between neighboring actuators become severe for large magnitudes and large variations of control signals that are necessary to produce certain wave-front shapes. If these nonlinearities are not properly modeled and if nonlinear effects are not properly addressed by control algorithms, then the overall wave-front correction performance can significantly be degraded.
Generally speaking, we can distinguish between closed-loop and open-loop approaches for DM control. Closed-loop approaches use DM-influenced measurements of system performance, such as wave-front observations or image quality measures, as feedback to the controller to facilitate iterative convergence to the desired control objective or wave-front shape. Open-loop approaches attempt to produce a desired wave-front shape in a single step without feedback or iteration. The performance of open-loop approaches can significantly be degraded in the presence of model uncertainties and external disturbances affecting the DM dynamics. On the other hand, closed-loop approaches can work relatively well even in the presence of the DM model uncertainties and disturbances, however, at the expense of larger numbers of wave-front observations and control iterations. In both approaches, improved model fidelity between control inputs and desired shapes will result in improved controller performance. These facts motivate us to work on the development of control approaches that are based on accurate DM models estimated from the observed wave-front data.
There are a large number of articles addressing control of DMs and related wave-front estimation problems. A comprehensive survey of all DM control approaches goes well beyond the length limits and the scope of this article. Consequently, we mention only the most relevant or recent approaches. In [22], multivariate adaptive regression splines are used to develop an open-loop approach for DM control. In [23,24], a least-squares method is used to estimate parameters of a discretized Partial Differential Equation (PDE) model of a DM. On the basis of such a model, in [24] an open-loop DM control approach is developed. In [25], a second-order DM model is used to develop a closed-loop control method. In [26], a modified iterative learning control method [27,28] is used to adjust DM control actions during closed-loop operation. The approach presented in [26] explicitly penalizes control actions, and consequently, it can be used to prevent actuator saturation. In [20], an open-loop DM control approach is developed that uses a first-principle PDE model to calculate DM control actions. Further refinements and developments of the open-loop control approach are presented in [29].
The above-summarized approaches can also be classified on the basis of how they take into account DM nonlinearities. The approach presented in [22] uses nonparametric estimation procedures and statistical analysis to estimate a DM model that can capture nonlinear system behavior. However, this approach requires a large number of wave-front samples in order to accurately fit the DM model. In [23,24], linear DM and actuator coupling models are assumed. In [26], a DM model is linearized around a working point. When used in a feedback control setting, the approaches presented in [23,24,26] can partially compensate for unmodelled nonlinearities. However, in the open-loop control framework, such approaches might produce sub-optimal control performance. In [20] and in the follow-up article [29], some of the DM model nonlinearities are taken into account by using a nonlinear PDE model of plate deformation [30].
A widely used DM control paradigm is to assume a constant DM model, without updating it on the basis of wave-front data collected during control iterations. For example, in a feedback control setting, used in [23,24,26], only control actions are being updated, and the DM model is kept constant. On the other hand, adaptive control techniques developed by control engineers [31], aim at adaptively updating both control actions and system models on the basis of the collected data. This control paradigm can effectively deal with model uncertainties and nonlinearities. However, despite several benefits of adaptive control approaches, to the best of the authors' knowledge, the potential of these approaches for DM control is not thoroughly investigated. There are only several approaches available in the literature that to some extent use and exploit the adaptive control idea. Thus, in [32] a DM influence function is iteratively calibrated during the closed-loop correction process. In [33], an approach based on a recursive least-squares method [34] is used to adaptively estimate a DM model. However, the approaches presented in [32,33] do not account for actuator saturation, which can significantly degrade the correction performance. Furthermore, since the approach developed in [33] works in closed-loop, control signals might become correlated with measurement noise and disturbances, and consistency of the adaptive estimation algorithm might not be guaranteed [35]. In addition, closed-loop data might be non-informative or the persistency of excitation condition might not be guaranteed [34]. If not properly addressed, these phenomena can significantly limit the DM correction performance.
In this paper, we develop a simple and easy-to-implement adaptive control approach for accurate DM control. For a desired DM surface profile (shape) and using batches of observed surface profile data, the proposed method adaptively determines both a DM model (influence matrix) and control actions that produce the desired surface profile with good accuracy. In the first iteration, the developed method estimates a DM influence matrix by solving a multivariable least-squares problem. This matrix is then used to compute the control actions by solving a constrained least-squares problem. Then, the computed actions are randomly perturbed and applied to the DM to generate a new batch of surface profile data. The new data batch is used to estimate a new influence matrix that is then used to re-compute control actions. This procedure is repeated until convergence is achieved. In our experimental setup, the DM deformation is observed using a partitioned aperture wavefront sensor [36][37][38]. Our experimental results show that the developed control approach can achieve accurate correction despite significant DM nonlinearities. Using only a few control iterations, the developed method is able to produce a surface profile root-mean-square error that varies from 5 − 30 [nm] for most of the tested Zernike wave-front modes without using direct feedback control. These results can additionally be improved by using larger data batches and more iterations or by combining the developed approach with feedback control. Finally, as we experimentally demonstrate, the developed method can be used to estimate a DM model that can effectively be used for a single-step open-loop DM control.
We briefly summarize the main contributions of the developed approach. Since the DM model is not estimated using closed-loop data, the developed control method does not suffer from issues related to the possible correlations between inputs and disturbances, and moreover, we ensure that the persistency of excitation condition is satisfied. This might not be the case for the closed-loop approach developed in [33]. In contrast to other adaptive DM control approaches, our approach explicitly takes into account actuator control limits when computing control actions, and thus can avoid performance loss due to actuators working in the saturation regime. On the other hand, most of the approaches proposed in the literature assume that the "true" DM model, represented by the influence matrix, is constant and global. That is, independent from the operating DM regime. We show that such an assumption can limit the performance of the control algorithms. On contrary, we relax this assumption by allowing DM influence matrices and control actions to depend on the desired surface profiles. In this way, we can improve the correction performance and obtain local models that are valid in certain operating regimes and that can better approximate the mirror nonlinear behavior. Finally, as we experimentally demonstrate, our approach can serve as the basis for the development of high-performance open-loop DM controllers. In the present form, the developed approach is ideal for correcting modest-speed wave-front aberrations, such as aberrations induced by thermal phenomena [39][40][41][42][43][44][45][46][47]. The developed methods can easily be adapted for applications where faster correction is required by implementing the main steps using recursive least squares techniques or using steepest descent optimization methods.
This manuscript is organized as follows. In Section 2, we briefly describe the used DM and experimental setup. In this section, we also quantify the DM nonlinearities and estimate an exponential model that relates control actions and deformations. In Section 3, we present the control methods. In Section 4 we present the experimental results, and in Section 5, we present conclusions and future work.

DM nonlinearity analysis and preliminary estimation
In this section, we first briefly describe the used MEMS DM and analyze its nonlinear behavior and actuator couplings. We also estimate an analytical form of a function relating control actions and mirror deformations immediately under actuators. We show that an exponential function model provides a more optimal deformation fit compared to the quadratic function that is often postulated in the literature. The presented analysis enables us to postulate and estimate a DM model in the next section, as well as to develop the control algorithm.
Experiments in this paper are performed on a Boston Micromachines MEMS DM. The actuators are distributed on a 12 by 12 grid. Corner actuators are inactive, so in total, we can compute control actions for 140 actuators. The DM pitch is 400 [µm] and the stroke is 3.5 [µm]. The DM facesheet thickness is 3 [µm]. The mirror surface profile (shape) is observed over a circular mirror region (pupil) that partly covers the grid of 10 by 10 actuators (see the experimental results in Section 4). Some of the controlled actuators are outside of the defined pupil area. This is a typical practice in DM control. Namely, by locating some of the active actuators outside of the observed pupil we can improve the quality and fidelity of the shapes formed within the pupil, and additionally, in this way, we can decrease the negative influence of the edge effects. The behavior of this mirror was analyzed and characterized in a number of papers, see for example [20,29] and references therein. Consequently, we will analyze only nonlinear mirror behavior that is relevant to our work.
We develop a fairly generally applicable method for DM control. The only requirement for the implementation of the developed control method is that wave-front or surface shape observations are available in vector or matrix forms. The main ideas and numerical implementation of the developed method are independent of the used sensor type, camera, and configuration of the optical system, and consequently, it can easily be used in other optical systems. Consequently, for brevity and due to space limitations, we only provide a brief description of the optical system used to test the developed method. The mirror deformation is observed using a Partitioned Aperture Wave-front (PAW) sensor [36][37][38]. This sensor type has a relatively large dynamic range. It works with uncollimated light sources, it has a high resolution that is limited by camera pixels, and since it performs only a single observation of the surface profile, it is relatively fast. Additional advantages of the PAW sensor are that it is robust (no mechanical moving parts), speckle-free, and polarization-independent. A detailed description of the PAW sensor and its working principle are given in [37,38]. The DM and sensor are integrated in an optical setup with the source provided by a red LED (660 [nm], Thorlabs) followed by a system of optical components that direct and shape the beam, and illuminate the mirror surface. The PAW sensor uses a monochrome camera (Blackfly BFS-U3-123S6M-C). The sensor output is an image (represented by a 1001 by 999 matrix) corresponding to the DM surface deformation.
Our first task is to analyze and quantify actuator coupling nonlinearities. An ideal DM will have a linear response that is mathematically described in the sequel. Let w(x, y, u 1 , u 2 , . . . , u 140 ) ∈ R be a function describing the deformation of the DM surface. The scalar variables x, y ∈ R are spatial positions in the mirror plane and u 1 , u 2 , . . . , u 140 are the DM control inputs. Namely, in the ideal case of a linear DM, the superposition principle should hold true. The superposition principle for two actuators is mathematically described as where u i and u j are control inputs applied to arbitrary actuators i and j, respectively, w(x, y, u i , u j ) is the deformation when both actuators are active at the same time, w(x, y, u i ) is the deformation when only the actuator i is active, and similarly, w(x, y, u j ) is the deformation when only the actuator j is active. Any deviation from this equation implies that the response is nonlinear. Consequently, deviations from the Eq. (1) can be used to identify and quantify the coupling nonlinearities of two actuators. We define the superposition error as We test the superposition principle for several adjacent actuators. The superposition test for the actuators i and j consists of the following experiment. First, we apply a control action to actuator i, while keeping the actuator j inactive, and observe the response. Then we reverse the order. That is, we apply a control action to actuator j, while keeping the actuator i inactive, and observe the response. Finally, we apply the control actions to both actuators and we observe the response. This data enables us to compute the left-hand and right-hand sides of the Eq. (1) and the superposition error given by (2). Figure 1 shows the results of the superposition test for several adjacent actuators. Panels (a), (c), and (e) show the mirror deformation in the vicinity of the actuator 50 (located near the mirror center). Panels (b), (d), and (f) show the mathematically superimposed response and the observed mirror deformation response. The control actions applied to every actuator are 70% of the maximal control-action magnitude. We can observe that the actuator coupling nonlinearity is significant for the actuator in the immediate vicinity of the actuator 50 (spacing equal to 1), and the nonlinearity effect fades away as the actuator spacing increases. The same phenomenon is observed for other tested actuators. This implies that the actuator coupling is strongly nonlinear for neighboring actuators and the coupling becomes more and more linear as the distance between actuators increases.
Next, we investigate the superposition error (2) for two adjacent actuators when the controlaction magnitudes vary from 0% to 100% of the maximal control-action magnitude. The results are shown in Fig. 2(a). We can observe that the superposition error increases as the control-action magnitude increases. We can observe that the superposition error is larger than 150 [nm] for 100 % of control-action magnitude. Finally, we experimentally investigate the shapes of functions relating the mirror deformations immediately under the actuators and applied control inputs. To collect the data, we vary the control-action magnitudes of every actuator from 0% to 100%, and we observe the deformations under each of the 140 actuators. The results are shown in Fig. 2(b). This figure shows the deformations of mirror points that correspond to approximate centers of 140 actuators. The above-presented superposition and nonlinear analyses imply that the nonlinear effects will be more dominant for larger magnitudes of control actions. On the other hand, we know that such magnitudes are required to produce wave-front shapes with larger Peak-to-Valley (P-V) ratios. Also, nonlinear effects will become more dominant when neighboring actuators have to be excited at the same time (since this will increase the superposition error). From Fig. 2(b) it follows that the actual actuator deformation response strays from the paths of the fitted quadratic or exponential functions for higher values of applied control-input magnitudes. Since the exponential function will be included in the model presented in the next section, we can expect that model uncertainties will increase for actuator actions close to the maximal actuation limits.
Next, we present an effective procedure for fitting an exponential function modeling actuator deformation responses. This function will be integrated into the DM model developed in the next section. A large number of DM modeling approaches available in the literature assume a quadratic model describing the actuator deformation response, see for example [23]. However, our experimental results show that the quadratic model is not the most accurate one (at least for the DM used in our experiments). Instead, the experimental observations imply that an exponential function represents a more accurate model. The exponential function approximates the absolute deformation as follows where k 1 and k 2 are constants that we want to estimate, and w a,i is the absolute value of the deformation of a mirror point located close to the center of the actuator i. We estimate the model constants using the observed responses of all actuators. We proceed as follows. For a known value of u i , the dependence of the function w a,i on k 1 and k 2 is nonlinear. Consequently, it seems that to estimate these constants we would need to employ a nonlinear squares method. However, by taking the natural logarithms of both sides of (3), we can transform this function into a form that enables us to use a linear least-squares method To every actuator i we apply 10 values of the control-input magnitude varying from 10% to 100% of the maximal control-input magnitude with an increment of 10%. We scale-down these values to 0.1, 0.2, . . . , 1. Let w a,i,j , j = 1, 2, . . . , 10, be the observed deformation for the scaled applied control-input magnitude of 0.1j. Then, on the basis of (4), we can write where q i ∈ R 10 , D ∈ R 10×2 , and x ∈ R 2 is a vector that we want to estimate. For i = 1, 2, . . . , 140, from (5), we have We determine x by solving the following least-squares problem where ∥·∥ 2 is the vector 2-norm. The solution of the least-squares problem is given by Once the solution x is found, we can determine the constants k 1 and k 2 as follows: k 1 = e x 1 and k 2 = x 2 , where x 1 and x 2 are the first and second entries of x. In fact, to postulate the DM model in the next section, we only need the constant k 2 , since the constant k 1 is absorbed as a parameter in the model and estimated again. On the basis of our experimental data, we have obtained: k 1 = 2.2473 and k 2 = 1.7420. Figure 2(b) shows the results of fitting the exponential function. Since the quadratic model is often used in the literature to fit the actuator response, for comparisons, we also show the results of fitting the quadratic function (fitted using a similar procedure to the one explained above). We can observe that the exponential function more optimally fits the actuator response compared to the quadratic function. This improved fitting accuracy will directly improve the DM correction performance. In the next section, we postulate the DM model that partly eliminates the effect of the exponential nonlinearity of the actuator deformation response. This enables us to use a linear multivariable least-squares method to estimate the DM model.

Model formulation and control method
In this section, we first formulate the DM model. Then, on the basis of the formulated model, we present the adaptive control method that for the desired DM surface profile shape iteratively updates both the DM model and control inputs. Finally, at the end of the section, we explain how the developed approach can be used to estimate a DM model that can be used for a single-step open-loop DM control.

Motivation, model formulation, and simple adaptive control
Let w(i∆x, j∆y) denote the observed mirror deformation (here, for notation simplicity, we omit the dependence on control inputs) at the spatial position i∆x, j∆y, where ∆x and ∆y represent the step sizes in the x and y directions of the masked observation grid (pupil). Next, let w ∈ R n be a vector that collects all the observed deformations w(i∆x, j∆y), where n is a total number of observation points. As experimentally verified in the previous section, the mirror behavior is nonlinear and for special cases of boundary conditions and applied loads, it is possible to derive the deformation as a function of applied forces. However, this function involves series of trigonometric functions [30], and as such, it is not suitable for the development of computationally efficient control methods. In contrast, our approach is to postulate a simplified model that can be adaptively estimated using a multivariable least-squares method.
In order to motivate the development of our method, in the sequel, we point out several drawbacks of a mirror modeling approach widely used in the literature. In the literature, the following mirror model is widely used and where g ∈ R 140 , Q const ∈ R n×140 is a constant influence matrix, and u i , i = 1, 2, . . . , 140 are control inputs applied to actuators. We assume that the control inputs are scaled down such that they belong to the interval [0, 1], where 0 represents no control action, and 1 represents 100% of control action. The matrix Q const is usually estimated from the wave-front measurements by formulating a multivariable least-squares problem. For example, in [26] we have used a multivariable least-squares method to estimate a constant influence matrix using batches of the input-output data. The main idea of the approach presented in [26] is to generate random input data and to excite random wave-front shapes (DM surface profiles). However, this and similar approaches have the following limitations. First of all, our experimental results show that randomly generated control inputs (for example, generated from a random uniform distribution) and the observed mirror deformations, produce estimates of Q const that in most cases result in significantly sub-optimal wave-front correction performance. The reason is that completely random data sets result in the influence matrix model that is only valid for a limited DM operating range that is determined by random fluctuations of the applied control inputs. We can actually say that completely random inputs might produce data sets that are not sufficiently information-rich to estimate a global influence function model (unless the shapes that we want to produce are actually in the operating range determined by random inputs). In fact, due to the nonlinear nature of the system, the model of the form (10) might not be the most optimal representation of the global model, or it might be difficult to postulate a global model. For example, in the previous section, we have seen that coupling nonlinearities become dominant when adjacent actuators are activated at the same time and when the actuation signals take larger values. Consequently, global linear models of the form (10) might be significantly suboptimal.
Another limitation of the approaches that rely upon the model (10), is that they assume a quadratic function of the applied control actions. This assumption might be valid for certain mirror types, however, it is more scientifically and practically justified to directly estimate the function form from the experimental observations. For example, our results, shown in Fig. 2(b) demonstrate that a more appropriate model is the exponential model.
To overcome these difficulties, we postulate a different mirror model. Let w d ∈ R n be the desired surface shape (surface profile) that we want to produce. We assume that around every desired shape that we want to produce, exists a model of the following form whereQ ∈ R n×140 is the influence matrix that to some extent depends on the shape we want to produce,ẑ ∈ R 140 is a vector consisting of exponents of optimal inputsû i , i = 1, 2, . . . , 140, that produce the desired shape, and k 2 = 1.7420 is the constant estimated in the previous section. It should be kept in mind that this constant value is valid only for the DM used in our experimental study and if other DMs are used, this constant needs to be re-estimated. It should be emphasized that we do not explicitly assume that the influence matrixQ is constant. We allow this matrix to depend to some degree on the desired shape that we want to produce. When the control inputs are not optimal for a certain shape they are simply denoted by u i , and the corresponding vector z is defined in the same spirit to (11) A few comments are in order. By defining the vector z (and its optimal equivalent) that consists of the terms u k 2 i , we are able to some degree to eliminate the input nonlinearities shown in Fig. 2(b). That is, if the influence matrixQ is known, then the model is linear with respect to the vectorẑ. However, the overall model (11) is not linear since both the influence matrix and the vectorẑ are parameters that need to be determined for the desired surface shape. Next, if we would know the influence matrixQ, then we would be able to invert the model (11) to estimate the vectorẑ, and consequently, to estimate the control inputsû 1 ,û 2 , . . . ,û 140 that produce the desired shape. However, the difficulty is that we do not know this matrix a priori. Our idea is to adaptively estimate the most appropriate form ofQ andẑ that will produce the desired shape. Before we provide the details of the proposed method, we introduce the following notation. Let the vector u s ∈ R 140 be an arbitrary control input vector applied to the DM at the discrete-time instant s, and let the vector z s ∈ R 140 be constructed from the entries of u s . That is, these vectors are defined by where u i,s is the control input applied to the actuator i at the discrete-time instant s. We have that u i,s = z 1/k 2 i,s . Let w s ∈ R n be a vector of mirror deformations produced by the control vector u s at the discrete-time instant s. Next, for s = 1, 2, . . . , p, we define matrices of batch data where U p ∈ R 140×p , Z p ∈ R 140×p , W p ∈ R n×p , and p is the number of collected data samples. For a chosen desired shape w d , we first perform the following preliminary estimation and control steps. First, we choose the number of data samples p. The minimum number of data samples is 140, in order to ensure the invertibility of matrices that are defined in the sequel. At every time sample s, we select the entries of the vector u s , s = 1, 2, . . . , p from the uniform distribution on the interval [0, 1] (note here that this selection is only used in the preliminary step). This selection ensures that the applied control actions are within the actuation limits, otherwise, the actuators will saturate and this will introduce significant control and estimation errors. Also, as explained later, this selection ensures that data matrices are invertible. Then, we apply the control inputs, collect the data, and we form the matrices Z p and W p . Then, we estimate a preliminary influence matrix Q (0) , by solving the following multivariable least-squares problem min where ∥·∥ F is the Frobenius norm and the superscript "(0)" in the notation Q (0) denotes the fact that this is a preliminary estimation step. The solution is given bŷ The fact that the control inputs are random ensures that the matrix Z p Z T p is invertible, and consequently, the solutionQ (0) is well-defined. Then, we compute the preliminary control actions u (0) by solving the following optimization problem where the relation symbol ≤ is applied element-wise, and 0 ∈ R 140 and 1 ∈ R 140 are vectors of zeros and ones, respectively. The constraint (18) ensures that the computed control actions are within the actuator physical limits. We solve the problem (17)-(18) using the interior-point method implemented in the MATLAB function lsqlin(·). Let the solution of (17)-(18) be denoted byẑ (0) . Then, the control inputs for the preliminary iteration are computed bŷ whereẑ (0) i is the ith entry ofẑ (0) . On the basis of computedû (0) i we form the control input vector u (0) . This control vector is then used to initialize the steps 1 and 2 that are described below. These steps are iteratively performed for l = 1, 2, 3, . . ., until convergence or until satisfactory control accuracy is achieved.
where α (l) s ∈ R 140 is a random vector whose entries are sampled from the standard normal distribution for every s, and ε (l) is a sufficiently small number that is gradually decreased for increasing l (for example the starting value for l = 1 can be 0.01 or smaller, for additional recommendations see the comments given at the end of this section). In order to ensure that the entries of the vectorũ (l) s do not exceed the saturation limits, we need to constrain them to the interval [0, 1]. Thus, if an entry of the vectorũ (l) s is larger than 1, we set that entry to 1. Similarly, if the entry is smaller than 0, set that entry to 0. After constructing the set of inputs (20), we apply them to the DM, and we observe the corresponding response s . Then, we can construct new batch matrices for s = 1, 2, . . . , h, as follows where the entries ofz (l) i , i = 1, 2, . . . , h, are computed by raising the entries ofũ (l) s to the power of k 2 . Next, we compute a new influence matrix by solving The solution is given by:Q The fact that the control inputs are randomly selected in (20) ensures that the matrix in the above expression is invertible and that the solution exists.

2.
Step 2: Compute control actions. On the basis ofQ (l) , form the following constrained least-squares problem, and compute the control actionsû (l) subject to: 0 ≤ z (l) ≤ 1.
Then similarly to (19), compute the control actionsû (l) i from the entries of the vectorẑ (l) that solves (24)- (25). Apply the computed control inputs and compute the surface profile error. If the results are not satisfactory or if the convergence is not achieved, return to step 1.
A few comments and recommendations for parameter selection of the proposed method are in order. First of all, the number of data samples for preliminary estimation and for forming the batch matrices, that is, the numbers p and h, should be larger or equal to the number of actuators, that is in our case 140. This is a necessary condition to ensure the invertibility of the corresponding data matrices for the estimation of the influence function. Then, the parameter ε (l) in (20) should be a relatively small number to ensure that the computed control actions do not exceed the saturation limits. Also, this number should gradually be decreased as the iteration l progresses, since the control actions will gradually converge, and we do not want them to be significantly perturbed such that the convergence is violated and the correction performance is decreased. We vary this parameter from 0.01 in initial iterations to 0.0001 in final iterations, with a proper increment. This parameter is a tuning parameter of the control approach. Our experimental results (see Section 4) show that the convergence is achieved after only several iterations. Here it should be emphasized that the initial estimation step, given by the Eq. (16), is practically independent of the desired wave-front shape. That is, to save data collection time, this step can be performed only once, and the preliminary estimates of the control actions should be then used as a seed for any desired wave-front shape.
Lastly, the above-presented approach is based on batches of input and wave-front data. However, the approach can be generalized such that only a single wave-front measurement and a single control input are used in every iteration l to update both the influence function and control actions. The main idea is to use either a recursive least squares algorithm or a steepest descent method to iteratively solve a least-squares problem similar to the problem given in (22). Then after every iteration, we can use the estimate of the influence matrix to solve (24) and to compute the control actions. However, special care should be taken to ensure the convergence of this approach. Finally, the influence function and control actions computed by the developed approach can be used to initialize a feedback control method to additionally improve the wave-front correction performance. Experimental verification of these ideas is a future research topic.

Open-loop control
Here, we explain how the developed approach can be used to generate open-loop control inputs. Namely, the basic idea is to use the above-presented approach to compute control inputs corresponding to mirror surface shapes belonging to a certain basis. For example, we can compute control actions that produce a set of Zernike basis shapes, or we can use some other 2D basis functions, such as 2D Walsh basis functions [48,49], depending on the application. We refer to these control actions as basis control actions. Then, we use these basis control actions to generate a wave-front data set that is used to estimate a global influence matrix Q ol . The basis control actions produce observed wave-front data sets that contain more useful information and variability compared to observed wave-front data sets produced by purely random control actions. Consequently, the data generated by basis control actions will contain more valuable information for estimating the influence matrix. This estimation can be done in the calibration phase. In the operating phase (open-loop control phase), for an arbitrary wave-front shape, that is not being used to generate basis control actions, we can compute open-loop control actions by solving a problem similar to (24)- (25), in which the matrixQ (l) is substituted by the estimated global matrix Q ol .
Let w b,1 , w b,2 , . . . , w b,v be the selected desired surface basis shapes, and letû b,1 ,û b,2 , . . . ,û b,v be the corresponding basis control actions that approximately produce these desired basis shapes and that are computed using the developed method. Here, v is the number of desired basis surface shapes. Then, for every computed basis control actionû b,i , we can generate batch matricesZ h,b,i andW h,b,i using the identical procedure to the procedure explained in the step 1 of the developed approach, see the Eqs. (20) and (21). Then, we can estimate a global influence matrix, by solving the following optimization problem The solution is given by: Once we have estimated this global matrix, we have completed the calibration step. Now, let w d,ol represent an arbitrary desired shape that we want to produce and that is generally accurately representable using the v wave-front basis shapes used in the calibration step. Then, we can compute the control actions that produce this desired shape by solving subject to: 0 ≤ z ol ≤ 1.
It should be emphasized that in (28)- (29), we are not using wave-front observations to produce the desired shape. All the model information is encapsulated in the estimated matrixQ ol . The solution computed by this open-loop approach can also be used to initialize feedback control algorithms.
In the next section, we experimentally demonstrate that the open-loop approach produces relatively good results.

Experimental tests
Here, we experimentally test the performance of the developed approach. First, we demonstrate the performance of adaptively estimating selected Zernike basis modes and we investigate the convergence of the developed method. Then, we experimentally investigate the performance of the open-loop approach. Finally, we experimentally verify our modeling assumption that the influence matrix models are inherently localized and dependent on the desired mirror surface shape. That is, we demonstrate that estimated influence matrices are significantly different for different desired surface shapes.

Simple adaptive control
We use p = 500 data samples for the preliminary estimation of the influence matrix, and in steps 1 and 2 of the developed method, we use h = 200 data samples. We vary the parameter ε (l) in (20) from 0.01 for the initial iterations to 0.0001 for the final iterations. Figures 3-6 show the control performance for the Zernike modes Z 0 2 , Z 1 3 , Z 3 5 , and Z 5 7 . Panels (a) and (b) in these figures show the desired and produced DM surface shapes, respectively. Panels (c) show the surface profile errors obtained as a difference between desired and produced DM surface shapes. Other details, such as P-V ratios and RMS errors are given in the captions of these figures. Panels (d) show the middle cross-sections of the desired and produced wave-front shapes. We can see that the correction performance varies from 10.2 to 44.1 [nm] RMS error for the full mirror surface and from 5.8 to 15.6 [nm] RMS error for the central part of the mirror surface. We can also observe that the errors are larger at the edges, and they become smaller and more uniform as we approach the central DM part. Next, we show the convergence of the developed approach. Figure 7 shows the convergence of the RMS surface error for generating the lower-order Zernike modes Z 0 2 , Z 1 3 , Z 0 4 , and Z 2 4 , for   the full (panel (a)) and for the central part (panel (b)) of the mirror surface. On the other hand, Fig. 8 shows the convergence of the RMS surface error for the higher-order Zernike modes Z 3 5 , Z 5 5 , Z 2 6 , and Z 5 7 , for the full (panel (a)) and for the central part (panel (b)) of the mirror surface. We can observe that the control performance slightly deteriorates for the higher-order Zernike modes. The experimentally observed convergence results confirm that completely random inputs (such as the ones generated from the uniform distribution) produce DM models that result in suboptimal correction performance. Namely, the iteration 1 in Figs. 7 and 8 corresponds to the preliminary iteration in which random control signals from the uniform distribution are used for estimating the DM influence matrix. We can observe that the RMS error value in iteration 1 is much larger than the RMS error values in the subsequent iterations. In the subsequent iterations, the control inputs and their random perturbations are becoming more optimized, and consequently, the estimated influence matrices are becoming more optimal for certain desired shapes. This experimentally confirms our statement that the global models of the form (10) estimated using completely random control inputs are generally suboptimal.

Open-loop control and influence matrix locality
We test the open-loop control approach by generating a desired surface shape that is a linear combination of the following Zernike basis modes where Z j i are the Zernike basis modes and c ij are the coefficients. We first use the developed adaptive control method to compute basis control actions (explained in Section 3.2) that produce the Zernike basis modes Z j i in (30). We group the coefficients in the set C = {c 20 , c 22   Next, we experimentally demonstrate the fact that the estimated influence functionsQ are dependent on the desired surface shapes. That is, the estimated influence functions are localized and together with control inputs, they define local models. Figure 11(a) and (b) show the 5000th row and 5th column, respectively, of the estimated influence matrices for the desired shapes  We can clearly observe that the entries of the estimated influence matrices have significantly different values. This shows that the estimated influence matrices are significantly different, implying that the DM models are localized.

Conclusions and future work
In this manuscript, we presented a novel data-driven method for model estimation and control of Deformable Mirrors (DMs). Using batches of control and observed DM surface data, and for the desired surface shape, the method iteratively estimates both the DM influence matrix and control actions that produce the desired shape. Our experimental results show that the developed control approach can achieve accurate correction despite significant DM nonlinearities. Using only a few control iterations, the developed method is able to produce a mirror surface root-mean-square error that varies from 5 − 30 [nm] for most of the tested Zernike modes without using direct feedback control. These results can additionally be improved by using larger data batches and more iterations or by combining the developed approach with feedback control. Finally, as we experimentally demonstrate, the developed method can be used to estimate a DM model that can effectively be used for a single-step open-loop DM control. In future work, we will use recursive estimation techniques and feedback control to increase the correction speed and performance of the developed method. Data availability. Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.