Dual-Update Data-Driven Control of Deformable Mirrors Using Walsh Basis Functions

Deformable Mirrors (DMs) are one of the main components of Adaptive Optics (AO) systems [1, 2]. A typical DM consist of a reflective optical surface that is deformed by a set of actuators. By precisely shaping the surface of a DM, we can compensate for wave-front aberrations in AO systems [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16].

In this paper, we consider the problem of developing control algorithms for DMs. There are a large number of approaches for DM control. A complete survey of all the methods goes well beyond the length limits of this manuscript. Consequently, we briefly mention only the most relevant or recent approaches. Most of the DM control approaches are based on the following control paradigm [17, 18, 19, 20, 21]. First, a model of the DM (influence function matrix) is estimated before the correction process. Then, such a model is either directly inverted, or it is used together with an iterative feedback control algorithm to calculate the DM control actions. That is, in most of the feedback control approaches, only control actions (control voltages) for the DM are updated on the basis of the observed wave-front aberrations, and the mirror model is kept constant. However, this approach that is based on linear and time-invariant DM models, might not be able to produce satisfactory wave-front correction performance in at least two scenarios that are encountered in practice. The first scenario is when the DM behavior changes over time. For example, an environment in which a DM operates might induce temperature fluctuations that can create thermo-mechanical deformations of mirror components and other associated phenomena that can significantly alter the DM dynamics [22, 23, 24, 25, 26, 27, 28, 12, 13, 29, 10]. This is especially the case for optical systems operating in space and for optical systems operating with high-power laser sources, where the absorbed heat increases the temperature of optical components. The second scenario is when the DM control actions significantly deviate from an initial linearization point (determined for example by bias control actions) that is used to obtain the linear model. Namely, most of the developed DM devices and prototypes exhibit a nonlinear behavior for sufficiently large magnitudes of control actions that are necessary to produce wave-front shapes with larger Peak-to-Valleys (PVs). That is, to correct for certain wave-front shapes, control actions need to significantly deviate from the nominal control actions that are used to obtain linearized DM models. In both of these scenarios, the performance of the DM control algorithms developed on the basis of time-invariant and linear models estimated for nominal operating conditions will gradually degrade.

On the other hand, it is a very well-known fact in control theory literature that control methods that dynamically update systems’ models on the basis of available observations are able to cope with time-varying dynamics, model uncertainties, and system nonlinearities. Furthermore, even in the case of linear systems, adaptive control methods are more accurate, faster, and in many cases, more optimal than classical feedback algorithms. All these facts motivate the development of a dual-update control algorithm in this paper. Under the term “dual-update” control algorithm, we understand a control algorithm that is able to update both control actions and a model (influence function) during DM operation. Here is it should be emphasized that these types of algorithms are often referred to as adaptive algorithms in control theory literature. However, since the word “adaptive” has different meanings in optics and control theory communities, we use the word “dual-update”, not to confuse an interested reader.

To the best of our knowledge, adaptive control approaches for DM control have received less attention and interest compared to the classical optimal control approaches [30, 16]. Thus, the approach proposed in [31] iteratively calibrates a DM influence function during the correction process. The approach developed in [32] uses a recursive least-squares method to dynamically estimate DM influence function during DM operation. The main limitation of the approaches developed in [32, 31] is that they do not explicitly take into account actuator saturation. A reliable control method has to take into account actuator physical constraints and has to be able to deal with actuator saturation. Recently, in [33], we have developed a DM control approach that updates the DM model on the basis of batches of observed wave-front data. This method takes into account actuator saturation. However, once the DM model is updated, the control is performed in an open-loop setting. Consequently, this method is not suitable for online control, where adaptive DM model updating is coupled with feedback control. Apart from these methods, adaptive filters have been used in [34] to predict and control wave-front disturbances.

A widely used approach for DM control is to express the observed wave-front using a Zernike modal basis. In this way, we implicitly perform the spatial model-order reduction of the DM control problem (the control problem is transformed from the spatial domain to the modal basis domain). On the other hand, for certain mirror types and for desired wave-front shapes consisting of steep peaks and valleys, Walsh basis functions might be a more suitable option than the Zernike basis functions [35, 36, 37]. Apart from this, in the general case, since Walsh basis functions are orthogonal, they can be used instead of Zernike polynomials in classical adaptive optics applications.

In this paper, we develop a novel dual-update DM control approach. The developed approach explicitly takes into account actuator constraints and it couples a feedback control algorithm with a recursive estimation of the DM influence function model. In this way, we are able to dynamically update the DM model and at the same time compute control actions that produce the desired shape. Furthermore, by expressing the desired and observed mirror surface shapes as sums of Walsh pattern matrices, we formulate the control problem in the 2D Walsh basis domain. We experimentally verify the developed approach on a 140 actuator MEMS DM, developed by Boston Micromachines.

The main contributions of this paper are summarized in the sequel. We propose a novel dual-update control approach for accurate DM control. Furthermore, in contrast to other dual-update or adaptive DM control approaches proposed in the literature, our approach explicitly takes into account actuator saturation. In this way, we can avoid the loss of performance that comes from actuator saturation. In addition, in this paper, we thoroughly experimentally investigate the performance of Walsh basis functions for DM control. This is important since only a handful of articles have explored the possibility of using Walsh basis functions for wave-front reconstruction and control [35, 36, 37], and the true potential of using Walsh basis functions for DM control is largely unexplored. It should be emphasized here that although we have performed experiments on a continuous face sheet DM, the developed control approach is applicable to other DM types. The developed approach is especially suitable for segmented DMs. Finally, we can easily modify the develop control algorithm to use Zernike basis functions instead of the Walsh basis functions.

This paper is organized as follows. In Section II, we present the procedure for approximating the mirror surface shape as a sum of Walsh pattern matrices. In Section III, we present the control method. In Section IV, we present the experimental results. In Section V, we present conclusions and briefly discuss future work.

In this section, we present a simple numerical procedure for approximating the mirror surface shape as a sum of Walsh pattern matrices.

Mirror surface deformation is usually represented by a matrix. That is, every entry of this matrix is a mirror deformation at a fixed spatial location. We refer to this matrix as the mirror deformation matrix. Let $W\in\mathbb{R}^{n\times n}$ be the mirror deformation matrix, where $n$ is the total size (measured in pixels) along $x$ and $y$ dimensions of the observed mirror surface. We decompose this matrix as follows

where $a_{p,q}\in\mathbb{R}$ are coefficients, and $Z_{p,q}\in\mathbb{R}^{n\times n}$ are Walsh pattern matrices with the entries that can either be $-1$ or $1$. The number $n$ should be selected such that $n=2^{V}$, where $V$ is a user-selected positive integer. The total number of the Walsh pattern matrices in (1) is equal to $M^{2}$, where $M\leq n$. We form the Walsh pattern matrices using Walsh basis functions. Here it should be emphasized that since this paper presents a proof of concept and due to mathematical simplicity and brevity, we use Walsh basis functions defined over a square domain. The mirror surface shape can also be represented using Walsh basis functions defined over a circular domain (polar Walsh basis functions) [38]. Consequently, the developed method can easily be used in the case of circular correction domains. However, even without using polar Walsh basis functions, with some modifications, the used approach that is based on square-domain Walsh basis functions, can be used for wave-front correction over circular correction domains. On the other hand, the control approach that is developed in Section III is practically independent of the type of basis functions for expanding the mirror surface shape. Thus, instead of using Walsh basis functions, we can also use Zernike basis functions in the developed control approach.

In the sequel, we first introduce a procedure for constructing the Walsh pattern matrices $Z_{p,q}$. Then, we introduce a procedure for computing the coefficients $a_{p,q}$. First, we choose the constant $V$. We select the constant $V$ such that a deformation submatrix with the dimensions of $2^{V}$ by $2^{V}$ pixels is within the limits of the maximal active mirror surface area that is observable by the used sensor (more details about the sensor used in our experiments can be found in Section IV). In our case, we use $V=8$ or $V=9$, giving us the submatrices with the dimensions of $256$ by $256$, and $512$ by $512$, respectively (see Section IV for more details).

Let the entries of the vector $\boldsymbol{\gamma}_{i}^{(V)}\in\mathbb{R}^{n}$, $n=2^{V}$, represent the values of the Walsh function $i$ of the order $V$. Here the index $i$ takes the values from $1$ to $2^{V}$. For example for $V=2$, the vectors $\boldsymbol{\gamma}_{i}^{(2)}$, $i=1,2,3,4$, representing the Walsh basis functions take the following form

The vectors $\boldsymbol{\gamma}_{i}^{(V)}$ can easily be constructed by using the MATLAB function $\text{hadamard}(\cdot)$. The rows of a matrix returned by this function represent Walsh basis functions. However, the Walsh basis functions that are represented by the rows of this matrix are not arranged in increasing order. Consequently, the rows of the matrix returned by the function $\text{hadamard}(\cdot)$ need to be permuted. For more details, see the MATLAB tutorial page [39] explaining the construction process of the Walsh basis functions.

Once the vectors $\boldsymbol{\gamma}_{i}^{(V)}$ are constructed, for selected $V$, we use the following procedure to construct the Walsh pattern matrices. First, we need to select the constant $M$. While selecting the constant $M$ we have to keep in mind several competing factors. Generally speaking, we have to make a trade-off between representation accuracy that is increased by increasing $M$, and dimensions of the matrices of the control algorithm that increase with the factor of $M^{2}$. Larger matrix dimensions increase the computational and memory complexities of the decomposition process, as well as the computational and memory complexities of the control and estimation algorithms that are introduced in Section III. We have tested the control algorithm for $M$ up to $120$ and this value is sufficient for desired shapes used in our experiments. Our lab computer has $64$ GB RAM, and for $M\geq 120$, MATLAB programming language that is used to control the DM, runs out of memory. Larger values of $M$ are possible if the control algorithm is implemented in a more memory-efficient way. One of the possible pathways to decrease the computational and memory complexities is to exploit the structure of the control matrices using approaches similar to the approaches developed in [40, 41, 11, 42, 43, 44, 45, 46, 47, 48]. However, this requires additional research, implementation, and testing efforts that are left for future research.

Once we have selected $V$ and $M$, we can proceed with the construction of the pattern matrices. We perform the following steps:

1. 1.

   Step 1: In this step, we construct the matrices $Z_{1,q}\in\mathbb{R}^{n}$, where $q=1,2,\ldots,M$, and $n=2^{V}$. The matrix $Z_{1,q}$ is constructed by transposing the vector $\boldsymbol{\gamma}_{q}^{(V)}$ and by stacking the newly formed row vectors on top of each other $n$ times:

   $\displaystyle Z_{1,q}=\begin{bmatrix}\Big{(}\boldsymbol{\gamma}_{q}^{(V)}\Big{)}^{T}\\ \Big{(}\boldsymbol{\gamma}_{q}^{(V)}\Big{)}^{T}\\ \vdots\\ \Big{(}\boldsymbol{\gamma}_{q}^{(V)}\Big{)}^{T}\end{bmatrix}.$

   (3)

2. 2.

   Step 2: In this step, we construct the matrices $Z_{q,1}$, where $q=1,2,\ldots,M$. The matrices $Z_{q,1}$ are constructed by transposing the matrices $Z_{1,q}$ that are formed in step 1, that is, $Z_{q,1}=Z_{1,q}^{T}$.

3. 3.

   Step 3: In this step, we construct the matrices $Z_{p,q}$, for the indices $p\geq 2$ and $q\geq 2$. The matrix $Z_{p,q}$ is calculated as follows

   $\displaystyle Z_{p,q}=Z_{p,1}\odot Z_{1,q},$

   (4)

   where $\odot$ is the matrix element-wise product (Hadamard matrix product).

A few comments about this construction procedure are in order. The matrix $Z_{p,q}$ can be seen as a $p,q$ block of a large block matrix. The first block row of this matrix is constructed in Step 1. In Step 2, we construct the first block column, where every block is a transpose of the corresponding matrix in the first block row. Then, in step 3, we construct the remaining block matrices by simply multiplying the matrices belonging to the first block column with the matrices belonging to the first block row. The constructed 2D Walsh pattern matrices are shown in Fig. 1 for $M=3$.

Next, we explain the decomposition process of the matrix $W$, that is, we explain how to compute the coefficients $a_{p,q}$ in (1) for known $W$. First, we vectorize the equation (1). The vectorization process is done by introducing the vectorization operator $\text{vec}(\cdot)$ [49]. Let $X\in\mathbb{R}^{n\times n}$ be an arbitrary matrix with the column vectors denoted by $\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{n},\in\mathbb{R}^{n}$. The vectorization operator is defined by

That is, the vectorization operator produces the vector $\mathbf{x}\in\mathbb{R}^{n^{2}}$ obtained by stacking column vectors $\mathbf{x}_{i}\in\mathbb{R}^{n}$ of the matrix $X$ on top of each other. By applying the vectorization operator to (1), we obtain

where $\mathbf{w}=\text{vec}\big{(}W\big{)}$, $\mathbf{w}\in\mathbb{R}^{n^{2}}$, and $\mathbf{z}_{p,q}=\text{vec}\big{(}Z_{p,q}\big{)}$, $\mathbf{z}_{p,q}\in\mathbb{R}^{n^{2}}$. On the other hand, due to the fact that the pattern matrices $Z_{p,q}$ are formed on the basis of the Walsh functions, we have that

The property (8) comes from the fact that the Walsh basis functions form a complete orthogonal set of functions. This property can be used to retrieve the coefficients $a_{p,q}$. Namely, from (7) and (8), we have

The equation (9) enables us to construct a projection matrix that produces the coefficients of the expansion (1). Let the matrices $\Psi_{i}\in\mathbb{R}^{M\times n^{2}}$ and $\Pi\in\mathbb{R}^{M^{2}\times n^{2}}$, and the vectors $\mathbf{a}_{i}\in\mathbb{R}^{M}$ and $\mathbf{a}\in\mathbb{R}^{M^{2}}$, be defined by

The vector $\mathbf{a}$ groups all the coefficients of the expansion (1). Then, using this construction, from (9) we have

The equation (11) shows that to express the deformation matrix $W$ as a sum of the Walsh pattern matrices, we just need to vectorize this matrix and to multiply the result with the projection matrix $\Pi$. That is, to compute the expansion coefficients $a_{p,q}$, we need to perform a single vector-matrix operation.

In this section, we present the control method. The basic idea of the control method is to update both the DM control actions and the DM influence function model. To develop the control method, we use the Walsh basis function expansion presented in Section II.

Let $W_{D}\in\mathbb{R}^{n\times n}$ be a desired mirror shape that we want to produce, and let $\mathbf{w}_{D}=\text{vec}\big{(}W_{D}\big{)}$, $\mathbf{w}_{D}\in\mathbb{R}^{n^{2}}$. Using (11), we compute the desired set of coefficients as follows

where $\Pi$ is the projection matrix introduced in (10) and $\mathbf{a}_{D}\in\mathbb{R}^{M^{2}}$ is the vector grouping the desired coefficients.

We send control actions to the DM or observe its deformation response at discrete-time instants, denoted by $k\in\mathbb{N}_{0}$. Let $\mathbf{w}_{k}\in\mathbb{R}^{n^{2}}$ be the observed mirror deformation in the vectorized form, that is $\mathbf{w}_{k}=\text{vec}\big{(}W_{k}\big{)}$, where $W_{k}\in\mathbb{R}^{n\times n}$ is the observed mirror deformation matrix at the discrete-time instant $k$. By using the projection (11), we have

where $\mathbf{a}_{k}\in\mathbb{R}^{M^{2}}$ are the coefficients. We postulate the following DM model

where $Q_{k}\in\mathbb{R}^{M^{2}\times r}$ is the influence matrix, and the vector $\mathbf{g}_{k}\in\mathbb{R}^{r}$ is defined by

and where $u_{i,k}$, $i=1,2,\ldots,r$, is the control input applied to the $i$th actuator at the time instant $k$, $r$ is the number of DM actuators (in our case $r=140$), and $\beta=1.7420$ is an estimate of the constant of the exponential dependence between DM control actions and the observed deformation response. This estimate is obtained using the linearized least-squares method explained in [33]. The entries of the vector $\mathbf{u}_{k}\in\mathbb{R}^{r}$ are the control inputs $u_{i,k}$. The values of the control input $u_{i,k}$ are in the interval $[0,1]$, with zero corresponding to no control action, and $1$ corresponding to the maximal control action applied to the actuator $i$ (maximal control voltage applied to the actuator). The vector $\mathbf{d}_{k+1}\in\mathbb{R}^{M^{2}}$, that is not known a priori, takes into account the measurement noise and unmodeled mirror behavior that is not captured by the model $Q_{k}\mathbf{g}_{k}$.

The control method consists of two steps. In the first step, we estimate an initial value of the influence matrix and compute the initial values of control inputs. These values are used to initialize the adaptive control algorithm. The adaptive control algorithm combines recursive estimation of the influence matrix with a novel control algorithm. In the sequel, we explain the first step.

In this section, we present the experimental results.

We test the developed approach using a Boston Micromachines MEMS DM with $r=140$ actuators. The actuation grid is $12$ by $12$ with all $4$ corner actuators being inactive. The DM stroke is $3.5$ $[\mu m]$ and the pitch is $400$ $[\mu m]$. The behavior of this DM type has been analyzed in many papers, see for example [53, 54] and follow-up works. Consequently, due to paper brevity, we do not further summarize other mirror properties.

The produced mirror surface shape is sensed by a Partitioned Aperture Wave-front (PAW) sensor  [55, 56, 57]. This sensor has a large dynamic range, it is relatively fast, and it operates with uncollimated light sources. It has a relatively high resolution that is only limited by the used camera. In addition, this sensor is speckle-free, robust, and polarization-independent. Further details related to this sensor can be found in [56, 57]. The used experimental setup is the same as the experimental setup used to generate the results in our previous publication [28]. Consequently, in the interest of brevity, we just mention its basic configuration. The light source is a red LED (660 [nm], Thorlabs). A system of optical components is used to direct and shape the beam and to illuminate the DM surface. A monochrome camera (Blackfly BFS-U3-123S6M-C) is used as a PAW sensing element. The maximal size of the observed image (deformation matrix) is $1001$ by $999$. The DM and sensor are controlled using the MATLAB programming environment.

As explained in Section II, to decompose the observed deformation as a sum of Walsh pattern matrices, the deformation matrix size should be expressed as a power of $2$. On the other hand, the camera of the PAW sensor produces a deformation image size of $1001$ by $999$. This image covers an area that is larger than the active DM area. Taking all these things into consideration, we have at least two options for selecting the deformation matrix size. The first option is $256$ by $256$ matrix, and the second option is $512$ by $512$. The second option covers the centers of all the actuators (including $4$ corner inactive actuators). However, in the second option, a part of deformation caused by the edge actuators will not take part in the defined $512$ by $512$ deformation matrix. We investigate the performance of the developed method for both options.

In this paper, we have developed a novel approach for adaptive control of DMs. On the basis of the feedback information provided by the sensor, our method updates both the DM influence matrix and the control actions. In this way, we are able to generate RMS control accuracy of $14-40$ $[nm]$. These results can additionally be improved by tuning the parameters of the algorithm. Besides introducing a novel control method, we have also demonstrated a good potential of using Walsh basis functions for DM control. Walsh basis function can achieve its full potential for control of segmented DMs. Our approach can straightforwardly be applied to the control problems involving segmented DMs. In future work, we will focus on improving the performance of the developed approach by optimally tuning the control algorithm parameters, and on reducing the computational complexity of the developed approach.