Disturbance Modelling for Minimum Variance Control in Adaptive Optics Systems Using Wavefront Sensor Sampled-Data

Abstract

Modern large telescopes are built based on the effectiveness of adaptive optics systems in mitigating the detrimental effects of wavefront distortions on astronomical images. In astronomical adaptive optics systems, the main sources of wavefront distortions are atmospheric turbulence and mechanical vibrations that are induced by the wind or the instrumentation systems, such as fans and cooling pumps. The mitigation of wavefront distortions is typically attained via a control law that is based on an adequate and accurate model. In this paper, we develop a modelling technique based on continuous-time damped-oscillators and on the Whittle’s likelihood method to estimate the parameters of disturbance models from wavefront sensor time-domain sampled-data. On the other hand, when the model is not accurate, the performance of the minimum variance controller is affected. We show that our modelling and identification techniques not only allow for more accurate estimates, but also for better minimum variance control performance. We illustrate the benefits of our proposal via numerical simulations.

1. Introduction

The last decades have marked the advent of the extremely large telescopes epoch in ground-based astronomy, in which astronomical observatories have experienced a growth in the aperture of their telescopes, see e.g., [1]. Currently, the Large Binocular Telescope (LBT) in Arizona, USA, reaches a combined effective aperture of $11.9$ m, whilst the Gran Telescopio Canarias (GTC) in the Canary Islands, Spain, reaches an effective aperture of $10.4$ m. This kind of modern telescopes is used for deep universe exploration as well as sharp high-definition images. However, the attainment of such images is subject to the inherent effects of the atmosphere and the vibrations of different equipment that are part of the telescope, see e.g., [2].

Adaptive optics (AO) is an optical technique that is used to improve astronomical images by compensating the effect of wavefront distortions caused by atmospheric turbulences and mechanical vibrations [3,4]. AO systems are very sensitive to vibrations acting in the propagation of the light [2,5,6]. Mechanical vibrations are typically induced by wind or elements within the instrumentation of the system, such as fans and cooling pumps. Hence, the modelling and mitigation of disturbances are a subject of increasing importance in many observatories [2,7,8,9,10]. This mitigation is achieved by deforming a deformable mirror (DM) in order to compensate the optical aberrations that are measured by a wavefront sensor (WFS), which, in turn, implies the need for both an accurate model of the complete AO system and adequate controllers.

In order to improve the quality of ground astronomical images, it is essential to obtain accurate models that define the dynamics of the plant that is comprised of the cascade connection of the DM and the WFS [11,12,13] and the accurate model parameters of all sources of noise and disturbances (turbulence and vibrations) that allow for implementing effective control techniques. For large telescopes (>8 m), disturbances become even more relevant, since the larger structures of this type of telescopes naturally oscillate according to different modes. In fact, disturbance mitigation is one of the key challenges for the next generation of extreme large telescopes [14,15].

In this paper, we focus on obtaining accurate disturbance models for the design of a minimum variance controller (MVC) in AO systems to improve the performance of the AO system. We present an extension of the work in [16,17,18,19] using continuous-time autoregressive (CAR) damped oscillators to model the disturbance and the Whittle’s likelihood technique to estimate the continuous-time parameters of the damped oscillators. We focus on the estimation of the continuous-time model parameters since the continuous-time model allows using different sampling intervals for the identification and the control scheme, which is necessary for more advanced control techniques, see e.g., [20]. We also analyse the impact of the disturbance model accuracy on the control performance, using a discrete-time transfer function of the disturbances for the design of the controller. From our analysis, we show that the control performance can drastically increase when the MVC is obtained from more accurate models. In our simulations, the MVC control performance improves in about $19\%$ when the controller is designed using the model that was obtained from our proposed identification technique.

The structure of the paper is as follows: In Section 2, a general description of a typical AO system is presented. In Section 3, we present an equivalent AO system model and the sampled-data model for the disturbances in AO systems. In Section 4, the proposed sample-data model for the disturbances is presented. In Section 5, we show the disturbance identification using the non-linear least square (NLS) fitting method. In Section 5, we also introduce our proposed identification algorithm that is based on the Whittle’s likelihood method to estimate the parameters of the continuous-time disturbance model. In Section 6, we present the MVC design and evaluate the MVC performance subject to model error. Numerical results are presented in Section 7, where we show the benefits of our proposed model and identification algorithm in terms of the MVC performance. Finally, we present conclusions in Section 8.

2. AO Systems

Figure 1 shows the classical model of an AO system used in astronomy. This AO closed-loop system consists of a wavefront sensor, a deformable mirror $\left(M_c\left(s\right)\right)$, a controller $\left(K\left(z^{-1}\right)\right)$, and a zero order hold (ZOH) [6,21,22], where $s=\frac{d}{dt}$ should be understood as the derivative operator or the argument of the Laplace transform, and z is either the forward shift operator $(z{\phi }_k={\phi }_{k+1})$ or the argument of the Z-transform. The continuous-time signals ${\phi }^{\mathrm{Tot}}\left(t\right)$, ${\phi }^{\mathrm{cor}}\left(t\right)$ and ${\phi }^{\mathrm{res}}\left(t\right)$ are the total disturbance (the amplitude of the disturbed wavefront; turbulence+vibrations), the correction disturbance (that corrects the wavefront), and the residual disturbance, respectively. The disturbed wavefront ${\phi }^{\mathrm{Tot}}\left(t\right)$ is defined as the sum of the atmospheric turbulence effect ${\phi }^{\mathrm{tur}}\left(t\right)$ and the different vibration sources effect ${\phi }^{{\mathrm{vib}}_l}\left(t\right)$, $l=1,\dots ,m$ [5,6,22]. That is:

$${\phi }^{\mathrm{Tot}}\left(t\right)={\phi }^{\mathrm{tur}}\left(t\right)+{\phi }^{{\mathrm{vib}}_1}\left(t\right)+{\phi }^{{\mathrm{vib}}_2}\left(t\right)+\cdots +{\phi }^{{\mathrm{vib}}_m}\left(t\right),$$

(1)

where m is the number of vibration sources. The signal ${\eta }_k$ shown in Figure 1 is a discrete-time additive zero-mean white Gaussian noise with variance ${\sigma }_{\eta }^2$.

2.1. Wavefront Sensor

A wavefront sensor (Figure 2) is an optical device that is used to measure the aberrations of an image. It consists of an array of lenses (lenslets) focused on an CCD or CMOS array (image plane). The basic idea is to measure the displacement of the centroid with respect to the ideal wavefront (planar wavefront) centroid position. Afterwards, using this information we can obtain a reconstructed wavefront (a linear approximation). Because only tilts are measured, the Shack–Hartmann wavefront sensor is not capable of detecting discontinuous steps in the wavefront. The CCD o CMOS detector requires a period of time (integration time) to accumulate enough photons to sample the wavefront and obtain a centroid. This implies that the WFS integrates the incoming wavefront during the integration time.

2.2. Deformable Mirror

A deformable mirror (DM) is an opto-mechanical reflective device that is capable of deforming the shape of its surface (see Figure 3). The shape of the DM can be controlled to correct the optical aberrations that are measured by a wavefront sensor in an AO system. This correction depends of the number of actuators, the actuator pitch (distance between actuator centres), the actuator stroke (maximum possible actuator displacement), the unfluence function i(characteristic surface shape corresponding to a single actuator), and the actuator coupling (displacement of the neighbours actuators).

2.3. AO Controller

The controller in an AO system sends the adequate signal to the deformable mirror to correct the disturbed wavefront that was measured by the WFS. Because the computation of the controller is demanding, the wavefront aberrations are represented by a linear combination of orthogonal basis elements. This linear combination allows for the separation of the control problem into independent modal controls, such as tip, tilt, and defocus [23,24].

The most used basis functions to describe aberrated wavefronts in optics are the Zernike polynomials, from which we obtain several orthogonal modes [23]. On the other hand, the varied nature of disturbances are perceived in different ways. For example, the mechanical vibration that is caused by excitations of on-site instruments usually exhibits narrow-band high-frequency components, the structural responses of large telescopes are within a low frequency band of natural resonances, and the atmospheric turbulences exhibit a (roughly) constant component below a corner frequency and a roll-off at high frequency on a Zernike basis [25]. In addition, the compensation of tip and tilt are one of the main tasks of AO [2,23,26]. Hence, an adequate identification of the disturbance source model is of great interest, which, in turn, allows for the development and implementation of effective control techniques.

A proportional-integral (PI) controller is a widely used control algorithm in AO [27,28,29]. However, with the construction of extremely large telescopes, the demanding requirements and challenging features increase, which causes the new generation of AO systems to require the implementation of more sophisticated identification and control techniques [1], such as linear-quadratic Gaussian (LQG) [1,5,7,29], minimum variance control (MVC) [26,30], or model predictive control (MPC) [31]. Notice that, in [26], the authors show that MVC is an equivalent representation of the typically used LQG controller for AO systems.

3. Disturbance Model in AO Systems

Typically, in the AO literature, the modelling and identification of disturbances in AO systems have been addressed using a second-order auto-regressive (AR(2)) discrete-time model with both time-domain data [26,29,32] and frequency-domain data [6,7,33,34,35]. In particular, in [35], an identification approach using a NLS fitting method was presented. This approach was successfully used to design a control strategy to mitigate the vibrations. Nevertheless, this model can present low accuracy at high frequencies (in the range $[4F_s/10,F_s/2]$, where $F_s$ is the sampling frequency). This, in turn, can lead to an unsatisfactory control performance in this high frequency range, especially for large telescopes. This behaviour is not desirable, since, for instance, high-order aberrations can be induced by misaligned components in an the AO system [23].

3.1. Equivalent AO System Model

In AO systems for astronomical observation, the WFS integrates the residual phase ${\phi }^{\mathrm{res}}\left(t\right)$ during a time interval $\left[t_{k-1},t_k\right]$, where $t_k$ is a sampling time, accumulating photons over a time interval $\Delta$ [22]. Thus, the discrete-time residual phase is given by:

$${\phi }_k^{\mathrm{res}}=\frac{1}{\Delta }{\int }_{t_{k-1}}^{t_k}{\phi }^{\mathrm{res}}\left(t\right)dt; t_k=k\Delta ,$$

(2)

where $k\in {\mathbb{N}}_0$ and $\Delta$ represent the sampling period. The corresponding discrete-time transfer function of the WFS in Equation (2) is given by [22]:

$$D\left(z^{-1}\right)=D_0{z}^{-\mu },$$

(3)

where $D_0$ is the gain and $\mu$ is the measurement delay.

Typically, the transfer function of the DM is modelled as the following discrete-time transfer function [5,11,22,26,35]

$$M\left(z^{-1}\right)=M_0{z}^{-\tau },$$

(4)

where $M_0$ is the gain and $\tau$ is the correction delay.

On the other hand, we consider that all of the elements in Figure 1 are linear. Subsequently, ${\phi }_k^{\mathrm{res}}$ is a linear function of the past values of ${\eta }_k$. The output of the system presented in Figure 1 is given by

$$\begin{array}{c}\hfill y_k=D_0{\phi }_{k-\mu }^{\mathrm{res}}+{\eta }_k,\end{array}$$

(5)

where ${\phi }_k^{\mathrm{res}}$ is the discrete-time residual phase that is given by:

$${\phi }_k^{\mathrm{res}}={\phi }_k^{\mathrm{Tot}}-{\phi }_k^{\mathrm{cor}},$$

(6)

and ${\eta }_k$ is a zero-mean white Gaussian noise. Subsequently, the variance of the output $y_k$ is given by:

$$\mathrm{var}\left\{y_k\right\}={\sigma }_{\eta }^2+D_0^2\mathrm{var}\left\{{\phi }_{k-\mu }^{\mathrm{res}}\right\},$$

(7)

thus, minimising the variance of $y_k$ is the same as minimising the variance of ${\phi }_k^{\mathrm{res}}$. Based on the equivalent block diagram for the AO system shown in Figure 1 that was presented in [26], a simple control system theory interpretation of this diagram is shown in Figure 4, where $K\left(z^{-1}\right)$ is the discrete-time transfer function of the controller, $M\left(z^{-1}\right)$ is the DM discrete-time transfer function, $D\left(z^{-1}\right)$ is the WFS discrete-time transfer function, $u_k$ is the output of controller, and $y_k$ is the system output. Note that $G\left(z^{-1}\right)$ is the equivalent discrete-time transfer function of the plant.

When considering the equivalent block diagram shown in Figure 4, the output signal is given by:

$$y_k=-D_0{\phi }_{k-\mu }^{\mathrm{cor}}+{\chi }_k,$$

(8)

where

$${\chi }_k={\eta }_k+D_0{\phi }_{k-\mu }^{\mathrm{Tot}},$$

(9)

corresponds to the disturbances of the AO system. Subsequently,

$$y_k=D_0\left[{\phi }_{k-\mu }^{\mathrm{Tot}}-{\phi }_{k-\mu }^{\mathrm{cor}}\right]+{\eta }_k.$$

(10)

Note that ${\phi }_k^{\mathrm{Tot}}$ can be modelled in different ways. In Section 3.2 we present the typical model utilised in AO systems. In Section 4, we present our proposed model for the disturbances.

3.2. Classical Sampled-Data Model for Disturbances in AO Systems

Typically the continuous-time AO disturbance model is expressed, as follows [22]:

$$\begin{array}{cc}\hfill {\phi }^{\mathrm{Tot}}\left(t\right)& =\sum _{l=0}^m\frac{{\beta }_l}{s^2+4{\zeta }_l\pi {\alpha }_{l}s+{\left(2\pi {\alpha }_l\right)}^2}{\stackrel{·}{\nu }}_l\left(t\right),\hfill \\ \hfill & =\sum _{l=0}^m\frac{{\beta }_l}{s^2+2{\zeta }_l{\varpi }_{l}s+{\varpi }_l^2}{\stackrel{·}{\nu }}_l\left(t\right),\hfill \end{array}$$

(11)

where ${\beta }_l$ is the gain, ${\alpha }_l$ (Hz) is the natural frequency, ${\zeta }_l$ is the damping coefficient, and ${\stackrel{·}{\nu }}_l\left(t\right)$ is a continuous-time zero-mean white Gaussian noise with variance ${\sigma }_l^2=1{\delta }_D\left(t\right),\forall l$ (${\delta }_D\left(t\right)$ is the Dirac delta). We assume that the noises ${\stackrel{·}{\nu }}_l\left(t\right),l=0,\cdots ,m$ are jointly uncorrelated. We use $({\varpi }_l=2\pi {\alpha }_l)$ for simplicity in the presentation.

In AO, the disturbances in discrete-time are typically modelled as:

$${\phi }_k^{\mathrm{Tot}}=\sum _{l=0}^m{H}_l\left(z^{-1}\right){\epsilon }_{l,k},$$

(12)

where ${\epsilon }_{l,k}$ is a zero-mean white Gaussian noise with variance ${\sigma }_{{\epsilon }_{l,k}}^2=1,\forall l$, jointly independent, $H_0\left(z^{-1}\right)$ is the discrete-time transfer function that is associated with the turbulence term and $H_l\left(z^{-1}\right)$, $l=1,\cdots ,m$, and the transfer function associated with the vibrations terms. The discrete-time transfer function $H_l\left(z^{-1}\right)$ is modelled using the following approximated auto-regressive (AR(2)) discrete-time representation [5,6,7,11,26,34,35]:

$$H_l\left(z^{-1}\right)=\frac{{\gamma }_l}{1-a_{1_l}{z}^{-1}+a_{2_l}{z}^{-2}},$$

(13)

where ${\gamma }_l$, $a_{1_l}$ and $a_{2_l}$ are the parameters for the autoregressive model. For the vibrations models, the parameters are defined by:

$$a_{1_l}=2e^{-{\zeta }_l{\varpi }_l\Delta }\cos \left({\varpi }_l\sqrt{1-{\zeta }_l^2}\Delta \right),$$

(14)

$$a_{2_l}=-e^{-2{\zeta }_l{\varpi }_l\Delta }.$$

(15)

Remark 1.

In the approximated model, the mathematical expression for the gains ${\gamma }_l$ in (13) are not provided in the AO literature. Thus, they have to be estimated from the experimental data.

The atmospheric turbulence is typically assumed to be time-invariant and statistically stationary for a fixed period of time [36,37,38]. This condition is referred to as the frozen-law approximation. It is widely understood that the model of the turbulence is very complex and a very accurate model will lead to a high computational load for AO systems. It has been shown in [6,39] that an adequate relaxation of the model can be achieved by fitting models that are simple enough for understanding (modelling) the corresponding phase distortion and carrying out fast computations (of the controllers), such as AR(1). In [40], it was shown that the AR(2) models better describe second order statistics of the atmospheric turbulence. In addition, in [20,34], an AR(2) model for the turbulences mode was proposed with the following parameters:

$$a_{1_0}=2e^{-{\zeta }_0{\varpi }_0\mathcal{T}}\cos \left({\varpi }_0\sqrt{1-{\zeta }_0^2}\mathcal{T}\right),$$

(16)

$$a_{2_0}=-e^{-2{\zeta }_0{\varpi }_0\mathcal{T}},$$

(17)

$${\varpi }_0=0.6\pi (\mathfrak{n}\left(i\right)+1)\frac{V_0}{\mathcal{D}},$$

(18)

where $\mathfrak{n}\left(i\right)$ is the radial order of the Zernike mode i, $\mathcal{T}$ is the exposure time of the system, $\mathcal{D}$ is the aperture diameter, and $V_0$ is the wind constant velocity. Hence, in this paper, we model the atmospheric turbulence as a second order system.

Remark 2.

In the case that the atmosphere can be modelled as a sum of second order systems, our identification technique could be used to adequately estimate such a model provided the numbers of second order systems is sufficiently large enough.

4. Proposed Modelling for Disturbances

In this section, we present a methodology to obtain an exact sampled-data model of the continuous-time disturbance model in Equation (11). The model in Equation (11) can be represented in state-space form as a CAR system:

$$\stackrel{·}{x}\left(t\right)=A_{c}x\left(t\right)+{\kappa }_c{\mathsf{\Psi }}_c\left(t\right),$$

(19)

$${\phi }^{\mathrm{Tot}}\left(t\right)=\stackrel{·}{b}\left(t\right)=C_{c}x\left(t\right),$$

(20)

where $A_c\in {\mathbb{R}}^{2m\times 2m}$, $x\left(t\right)\in {\mathbb{R}}^{2m\times 1}$, ${\kappa }_c\in {\mathbb{R}}^{2m\times m}$, ${\mathsf{\Psi }}_c\in {\mathbb{R}}^{m\times 1}$, $C_c\in {\mathbb{R}}^{1\times 2m}$, m is the number of damping oscillators used, and the noise variance is $Q_c={\kappa }_c{\kappa }_c^T$. These continuous-time system matrices are given by:

$$A_c=\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ -{\varpi }_1^2& -2{\zeta }_1{\varpi }_1& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 0& 1\\ 0& 0& \cdots & -{\varpi }_m^2& -2{\zeta }_m{\varpi }_m\end{array}\right],$$

(21)

$${\kappa }_c=\left[\begin{array}{ccccc}0& \cdots & \cdots & \cdots & 0\\ {\beta }_1& 0& \cdots & \cdots & ⋮\\ 0& 0& \ddots & \ddots & ⋮\\ ⋮& {\beta }_2& \ddots & \ddots & ⋮\\ ⋮& 0& \ddots & \ddots & ⋮\\ ⋮& ⋮& 0& \ddots & ⋮\\ ⋮& ⋮& ⋮& \ddots & 0\\ 0& 0& 0& 0& {\beta }_m\end{array}\right],\hspace{1em}\hspace{1em}{\mathsf{\Psi }}_c\left(t\right)=\left[\begin{array}{c}{\stackrel{·}{\omega }}_1\left(t\right)\\ {\stackrel{·}{\omega }}_2\left(t\right)\\ ⋮\\ {\stackrel{·}{\omega }}_m\left(t\right)\end{array}\right],$$

(22)

$$C_c=\left[ 1\hspace{1em}0\hspace{1em}|\hspace{1em}1\hspace{1em}0\hspace{1em}|\hspace{1em}\cdots \hspace{1em}|\hspace{1em}1\hspace{1em}0 \right],$$

(23)

where the vertical lines that are shown in Equation (23) represent the grouping of duple elements, corresponding to each term of the sum in Equation (11).

Note that it is necessary to define the auxiliary variable $\stackrel{·}{b}\left(t\right)$ to obtain the exact discrete-time model [41].

We consider that the WFS integrates the residual phase in order to obtain the sampled-data model for the continuous-time disturbance model in Equation (11):

$${\phi }_k^{\mathrm{res}}=\frac{1}{\Delta }{\int }_{t_{k-1}}^{t_k}{\phi }^{\mathrm{res}}\left(t\right)dt,$$

(24)

this could be understood as an averaging anti-aliasing filter (AAF) that is used before taking samples (see [17,42]). For this kind of sampling, we define a system with an extended state $\mathcal{X}\left(t\right)={\left[x^T\left(t\right)b\left(t\right)\right]}^T$ [41]. Subsequently, the continuous-time state-space model is as follows:

$$\stackrel{·}{\mathcal{X}}\left(t\right)=\mathcal{A}\mathcal{X}\left(t\right)+\left[\begin{array}{c}{\kappa }_c\\ 0\end{array}\right]{\mathsf{\Psi }}_c\left(t\right),$$

(25)

where

$$\mathcal{A}=\left[\begin{array}{cc}{A}_c& 0\\ C_c& 0\end{array}\right].$$

(26)

Thus, the corresponding discrete-time state equation model is given by [41]:

$${\mathcal{X}}_{k+1}=e^{\mathcal{A}\Delta }{\mathcal{X}}_k+n_{k+1},$$

(27)

if we consider that $A_c$ is invertible, we have the following result [43]:

$$e^{\mathcal{A}\Delta }=\left[\begin{array}{cc}{e}^{A_c\Delta }& 0\\ C_c{A}_c^{-1}(e^{A_c\Delta }-I)& \mathcal{I}\end{array}\right],$$

(28)

where $\mathcal{I}$ is the identity matrix and $n_{k+1}$ is a correlated discrete-time zero-mean white Gaussian noise with variance Q [41,44],

$$n_{k+1}={\int }_0^{\Delta }{e}^{\begin{array}{c}\hfill \left[\begin{array}{cc}{A}_c& 0\\ C_c& 0\end{array}\right]\xi \end{array}}\left[\begin{array}{c}{\kappa }_c\\ 0\end{array}\right]{\mathsf{\Psi }}_{k+1-\eta }d\xi .$$

(29)

The variance Q and the exponential matrix can be obtained using that presented in [45]:

$$e^{\begin{array}{c}\hfill \left[\begin{array}{cc}-\mathcal{A}& Q_c^{*}\\ 0& {\mathcal{A}}^T\end{array}\right]\Delta \end{array}}=\left[\begin{array}{cc}{P}_{11}& P_{12}\\ 0& P_{22}\end{array}\right],$$

(30)

where:

$$Q_c^{*}=\left[\begin{array}{c}{\kappa }_c\\ 0\end{array}\right]{\left[\begin{array}{c}{\kappa }_c\\ 0\end{array}\right]}^T,$$

(31)

then,

$$\begin{array}{cc}\hfill e^{\mathcal{A}\Delta }& =P_{22}^T,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill Q& ={\int }_0^{\Delta }{e}^{\mathcal{A}\xi }{Q}_c^{*}{e}^{{\mathcal{A}\xi }^T}d\xi =P_{22}^T{P}_{12},\hfill \\ \hfill & =\left[\begin{array}{cc}{Q}_{11}& Q_{12}\\ Q_{12}^T& Q_{22}\end{array}\right],\hfill \end{array}$$

(33)

where $Q_{11}\in {\mathbb{R}}^{2m\times 2m}$, $Q_{12}\in {\mathbb{R}}^{1\times 2m}$ and $Q_{22}\in {\mathbb{R}}^{1\times 1}$.

Thus, we have that

$$\left[\begin{array}{c}{x}_{k+1}\\ b_{k+1}\end{array}\right]=\left[\begin{array}{cc}{e}^{A_c\Delta }& 0\\ C_c{A}_c^{-1}(e^{A_c\Delta }-I)& \mathcal{I}\end{array}\right]\left[\begin{array}{c}{x}_k\\ b_k\end{array}\right]+n\left(t_{k+1}\right),$$

(34)

$$\begin{array}{c}\hfill x_{k+1}=A{x}_k+w_{k+1},\\ \hfill {\phi }_{k+1}^{\mathrm{Tot}}\Delta =b_{k+1}-b_k=C{x}_k+v_{k+1},\end{array}$$

(35)

where:

$$\begin{array}{cc}\hfill A& =e^{A_c\Delta },\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill C& =C_c{A}_c^{-1}\{A-\mathcal{I}\},\hfill \end{array}$$

(37)

and the signals $w_{k+1}$ and $v_{k+1}$ are correlated discrete-time zero-mean white Gaussian noises with covariance matrix Q.

Remark 3.

If we represent the plant $(G\left(z^{-1}\right)=M\left(z^{-1}\right)D\left(z^{-1}\right))$ in state-space form as:

$$\begin{array}{cc}\hfill {\stackrel{·}{x}}_p\left(t\right)& =A_{p_c}x\left(t\right)+B_{p_c},\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill {\stackrel{·}{b}}_p\left(t\right)& =C_{p_c}{x}_p\left(t\right),\hfill \end{array}$$

(39)

where a ZOH and an averaging AAF are utilised, then the corresponding discrete-time system is given by [41]:

$$\begin{array}{cc}\hfill x_{p_{k+1}}& =A_p{x}_{p_k}+B_p{u}_k,\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill \frac{b_{p_{k+1}}-b_{p_k}}{\Delta }& =C_p{x}_{p_k}+D_p{u}_k,\hfill \end{array}$$

(41)

where:

$$\begin{array}{cc}\hfill A_p& =e^{A_{p_c}\Delta },\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill B_p& ={\int }_0^{\Delta }{e}^{A_{p_c}\xi }{B}_{p_c}d\xi ,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill C_p& =\frac{1}{\Delta }{\int }_0^{\Delta }{e}^{A_{p_c}\xi }d\xi ,\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill D_p& =\frac{1}{\Delta }{C}_{p_c}{\int }_0^{\Delta }{\int }_0^{\xi }{e}^{A_{p_c}\varphi }d\varphi d\xi .\hfill \end{array}$$

(45)

On the other hand, the design of controllers, such as MVC, typically require only one discrete-time transfer function with one noise signal. In order to obtain this, we need to compute the spectral factorization $H\left(z^{-1}\right)$ [46] via numerical approximation, solving a Riccati equation [46,47,48]. Subsequently, the proposed model of the disturbances is as follows (for more details on spectral factorization, see [46] (pp. 71–78)):

$${\phi }_k^{\mathrm{Tot}}=H\left(z^{-1}\right)e_k,$$

(46)

where $e_k$ is a discrete-time zero-mean white Gaussian noise with variance ${\lambda }^2$. Note that both $H\left(z^{-1}\right)$ and ${\lambda }^2$ depend on the continuous-time parameters.

In order to develop an identification algorithm we define the discrete-time power spectral density (PSD) as follows [34,35]:

$$S\left(e^{\mathrm{i}{\omega }_j\Delta }\right)={\left|H\left(e^{\mathrm{i}{\omega }_j\Delta }\right)\right|}^2{\lambda }^2,$$

(47)

where $\mathrm{i}=\sqrt{-1}$, $j=1,2,\cdots ,N$ and $0<{\omega }_j\le 2\pi (F_s/2)$, and ${\omega }_j=\left(2\pi j\right)/\Delta$.

Finally, we summarize the procedure to obtain the discrete-time PSD shown in Algorithm 1.

Algorithm 1 Discrete-time PSD

1: Obtain the continuous-time state-space model of the disturbances using Equations (19) and (20) 2: Obtain the continuous-time state-space model of the extended system in Equation (25) 3: Compute the discrete-time state equation of the extended system in Equation (27) 4: Compute the discrete-time system matrices, A and C, of the disturbances using Equations (36) and (37) 5: Obtain the spectral factorization $H\left(z^{-1}\right)$ and ${\lambda }^2$ 6: Obtain the discrete-time PSD in Equation (47)

5. Identification of Disturbances

In this section, we present two forms to obtain the disturbance model parameters. First, we show the method presented in [35], in which the typical AR(2) discrete-time disturbance model is used (see Equation (13)). Next, we propose to use the Whittle’s likelihood method to obtain the continuous-time disturbance parameters using the proposed model from the previous section.

5.1. Nonlinear Least Square Fitting Method

In [35], the authors estimated the disturbance frequencies, the damping coefficients, and the variances of the discrete-time noises $\left(\epsilon \right)$ signals assuming the model in Equation (12) and that the discrete gains were all equal to one. Using this method, in this paper, we assume the noise variance ${\sigma }_{\epsilon }^2=1$ for all disturbance components, and we identify the disturbance frequencies, the damping coefficients, and the discrete gains. Despite these assumptions differing from the original approach in [35], it can be successfully implemented, since they are just a normalization procedure.

In order to formulate the NLS estimation algorithm, we define the vector of parameters to be estimated, ${\overrightarrow{\theta }}_{NLS}$, as follows:

$${\overrightarrow{\theta }}_{NLS}={\left[\begin{array}{cc}{\overrightarrow{\zeta }}^T& {\overrightarrow{\gamma }}^T\end{array}\right]}^T,$$

(48)

where $\overrightarrow{\zeta }={\left[{\zeta }_1{\zeta }_2\cdots {\zeta }_m\right]}^T$ is the vector that contains all of the damping coefficients, and $\overrightarrow{\gamma }={\left[{\gamma }_1{\gamma }_2\cdots {\gamma }_m\right]}^T$ is the vector that contains all of the discrete-time gains.

Following [35], we use the periodogram of the sampled-data and fit it with model ${\overline{S}}_l\left(e^{\mathrm{i}{\omega }_j\Delta }\right)$ in Equation (49). The different frequencies (Hz) in $\overrightarrow{\alpha }={\left[{\alpha }_1{\alpha }_2\cdots {\alpha }_m\right]}^T$ were found in the periodogram as the frequency peaks.

The discrete-time PSD is given by:

$${\overline{S}}_l\left(e^{\mathrm{i}{\omega }_j\Delta }\right)={\left|\frac{{\gamma }_l}{1-a_{1_l}{e}^{-\mathrm{i}{\omega }_j\Delta }-a_{2_l}{e}^{-2\mathrm{i}{\omega }_j\Delta }}\right|}^2,$$

(49)

and the NLS cost function $F_l\left({\theta }_{NLS}\right)$ is obtained from

$${ϵ}_l\left(e^{\mathrm{i}{\omega }_j\Delta }\right)=10\log \left[\frac{I_l\left(e^{\mathrm{i}{\omega }_j\Delta }\right)}{{\overline{S}}_l\left(e^{\mathrm{i}{\omega }_j\Delta }\right)}\right],$$

(50)

where $I\left(e^{\mathrm{i}{\omega }_j\Delta }\right)$ is the periodogram of the sampled-data series ${\phi }^{\mathrm{Tot}}$, and it is given by:

$$I\left(e^{\mathrm{i}{\omega }_j\Delta }\right)={\left|\frac{1}{\sqrt{N}}\sum _{t=1}^N{\phi }_t^{\mathrm{Tot}}{e}^{-\mathrm{i}{\omega }_j\Delta }\right|}^2,$$

(51)

then,

$$F_l\left({\overrightarrow{\theta }}_{NLS}\right)=\frac{1}{2}{ϵ}_l{\left(e^{\mathrm{i}{\omega }_j\Delta }\right)}^T{ϵ}_l\left(e^{\mathrm{i}{\omega }_j\Delta }\right).$$

(52)

Note that the periodogram must be divided in m sections according to the frequencies that were previously obtained, since each AR(2) system is fitted independently (not jointly) of the others.

Finally, the estimation problem can be written as the following optimisation problem:

$${\widehat{\theta }}_{NL{S}_l}=\arg \underset{{\overrightarrow{\theta }}_{NL{S}_l}}{\min }\left(F_l\left({\overrightarrow{\theta }}_{NLS}\right)\right),\mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}0<{\zeta }_l<1\\ {\gamma }_l>0\end{array}\right.,l=0,\cdots ,m.$$

(53)

5.2. Whittle’s Likelihood

The Whittle’s likelihood function is a widely used frequency-domain approximation of the traditional likelihood function [49,50,51].

It is well known that the (time-domain) Whittle log-likelihood function is given by:

$${\ell }_N\left(\overrightarrow{\theta }\right)=-\frac{1}{N}\sum _{j=1}^N\log S\left(e^{\mathrm{i}{\omega }_j\Delta }\right)-\frac{1}{N}\sum _{j=1}^N\frac{I\left(e^{\mathrm{i}{\omega }_j\Delta }\right)}{S\left(e^{\mathrm{i}{\omega }_j\Delta }\right)},$$

(54)

where N is the data length and $\overrightarrow{\theta }$ is the vector of parameters to be estimated that is given by:

$$\overrightarrow{\theta }={\left[\begin{array}{ccc}{\overrightarrow{\alpha }}^T& {\overrightarrow{\zeta }}^T& {\overrightarrow{\beta }}^T\end{array}\right]}^T,$$

(55)

$\overrightarrow{\alpha }={\left[{\alpha }_1{\alpha }_2\cdots {\alpha }_m\right]}^T$ is the vector that contain all the natural frequencies (in Hz), $\overrightarrow{\zeta }={\left[{\zeta }_1{\zeta }_2\cdots {\zeta }_m\right]}^T$ is the vector that contains all the damping coefficients, and $\overrightarrow{\beta }={\left[{\beta }_1{\beta }_2\cdots {\beta }_m\right]}^T$ is the vector that contains all the gains. In Equation (54) $S\left(e^{\mathrm{i}{\omega }_j\Delta }\right)$ is the discrete-time PSD of the disturbances given in Equation (47) and $I\left(e^{\mathrm{i}{\omega }_j\Delta }\right)$ is the periodogram that is given in Equation (51).

Notice that the Whittle log-likelihood is a function of the continuous-time system parameters (to be estimated) $\overrightarrow{\theta }$, since the spectrum $S\left(e^{\mathrm{i}{\omega }_j\Delta }\right)$ depends on them.

Afterwards, the ML estimation problem is given by:

$${\widehat{\theta }}_{\mathrm{ML}}=\arg \underset{\overrightarrow{\theta }}{\min }-{\ell }_N\left(\overrightarrow{\theta }\right),\mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}0<\overrightarrow{\alpha }<F_s/2\\ 0<\overrightarrow{\zeta }<1\\ \overrightarrow{\beta }>0\end{array}\right..$$

(56)

Finally, we summarize our identification technique in Algorithm 2.

Algorithm 2 Identification algorithm

1: Compute the periodogram $I\left(e^{\mathrm{i}{\omega }_j\Delta }\right)$ using Equation (51) with ${\omega }_j=\left(2\pi j\right)/\Delta$, and $j=1,2,\cdots ,N$ 2: Choose an initial guess ${\widehat{\theta }}^{\left(0\right)}$ 3: Obtain the discrete-time spectrum $S\left(e^{\mathrm{i}\omega \Delta }\right)$ in function of the vector of continuous-time parameters $\overrightarrow{\theta }$ using the Algorithm 1 4: Define the Whittle’s likelihood using Equation (54) 5: Solve the ML estimation problem in Equation (56) e.g., using Matlab® fmincon function

6. MVC Performance in AO Systems

6.1. Minimum Variance Control Design

In classical AO, the goal is to mitigate the wavefront distortions that are caused by atmospheric turbulences and mechanical vibrations in the astronomical images. That is equivalent to minimising the continuous-time residual phase variance, $\mathrm{var}\left\{{\phi }^{\mathrm{res}}\left(t\right)\right\}$[22]. However, in [22], the authors demonstrated that minimising the continuous-time residual phase variance is equivalent to minimising the discrete-time residual phase variance, $\mathrm{var}\left\{{\phi }_k^{\mathrm{res}}\right\}$.

The goal of MVC is to minimise the variance of the output signal $\mathrm{var}\left\{y_k\right\}$[42]. For simplicity of the presentation, we consider that the measurement noise, ${\eta }_k$, is an independent and identically distributed (iid) signal. We also assume that the atmospheric turbulence ${\phi }^{\mathrm{tur}}$, the vibrations ${\phi }^{\mathrm{vib}}$, and the measurement noise $\eta$, are mutually independent. In this setup, the measurement noise ${\eta }_k$ is assumed negligible. The atmospheric turbulence and vibrations will be modelled using the same model structure with different parameters. Thus the output of the AO loop is given by:

$$y_k=z^{-d}{D}_0{M}_0{u}_k+H\left(z^{-1}\right)e_k,$$

(57)

where $d=\mu +\tau$ and $u_k$ is the output of the controller. $H\left(z^{-1}\right)$ is the discrete-time transfer function obtained from the spectral factorization [46] of the spectrum of the signal ${\chi }_k$ in Equation (9). The signal $e_k$ is a zero-mean white noise with variance ${\lambda }^2$ that is also obtained from the spectral factorization procedure.

In order to obtain the MVC [42,52] we rewrite $H\left(z^{-1}\right)$ as

$$H\left(z^{-1}\right)=F\left(z^{-1}\right)+z^{-d}R\left(z^{-1}\right),$$

(58)

where $F\left(z^{-1}\right)$ is a finite impulse response (FIR) filter that corresponds to the first d samples of the impulse response of $H\left(z^{-1}\right)$, given by:

$$F\left(z^{-1}\right)=1+f_1{z}^{-1}+f_2{z}^{-2}+\cdots +f_{d-1}{z}^{-(d-1)},$$

(59)

where $f_{\iota }$ is the $\iota$th coefficient of the FIR filter $F\left(z^{-1}\right)$, and $R\left(z^{-1}\right)$ is a causal filter that can be obtained from Equation (58):

$$R\left(z^{-1}\right)=\left(H(z^{-1})-F(z^{-1})\right)z^d.$$

(60)

Finally the discrete-time transfer function of the minimum variance controller is given by [42,52]:

$$\begin{array}{c}\hfill K\left(z^{-1}\right)=\frac{R\left(z^{-1}\right)}{F\left(z^{-1}\right)D_0{M}_0},\end{array}$$

(61)

which corresponds to the controller in Figure 4.

6.2. Performance of MVC Subject to Model Error

It is well known that the performance of MVC is greatly affected by the accuracy of the model utilised in the design of the controller [42]. In order to analyse the effect of model error in its performance, we model the estimates of the disturbance transfer function in terms of the true value and an error term as:

$$\begin{array}{c}\hfill \widehat{H}\left(z^{-1}\right)=H\left(z^{-1}\right)\left(1+H_{\mathcal{E}}(z^{-1})\vartheta \right),\end{array}$$

(62)

where $H\left(z^{-1}\right)$ is the true discrete-time transfer function and $H_{\mathcal{E}}\left(z^{-1}\right)\vartheta$ is the multiplicative model error. The subscript $\mathcal{E}$ indicates multiplicative errors, whilst the parameter $\vartheta$ represents the gain of the multiplicative model error. When carrying out an identification procedure, the gain of the multiplicative model error can be set to $\vartheta =1$ in order to obtain and subsequently analyse $H_{\mathcal{E}}\left(z^{-1}\right)$. In some cases, a different set-up (e.g., a different measurement noise variance) could lead to the same $H_{\mathcal{E}}\left(z^{-1}\right)$, but with a different gain.

In

Table 1. The values used for six continuous-time damped oscillators.

${\mathit{\alpha }}_{\mathit{l}}\left(\mathrm{Hz}\right)$	2	14	20	29	43	60
${\mathbf{\zeta }}_{\mathbf{l}}$	0.9	0.05	0.05	0.05	0.05	0.05
${\mathbf{\beta }}_{\mathbf{l}}(\times {\mathbf{10}}^{\mathbf{4}})$	1.439	3.493	3.734	4.253	5.540	9.711

, we present the parameters of a continuous-time disturbance model that is comprised of 6 oscillators. To evaluate the performance of the MVC for this model, it is necessary to obtain the spectral factorization of ${\chi }_k$ by solving a Riccati equation [46,47,48].

In order to illustrate the effect of the error term on the estimated disturbance transfer function, we consider the continuous-time system (see Equation (11)) that is defined by the parameters in Table 1. Utilizing Algorithm 1, we obtain the corresponding PSD. We estimate the system model from the NLS method that is presented in [35] and, by fixing $\vartheta =1$, we obtain the multiplicative model error $\left(H_{\mathcal{E}}\left(z^{-1}\right)\right)$ using Equation (62). Finally, we select different values of $\vartheta$ and obtain the frequency response of the different multiplicative model errors. Figure 5 shows these frequency responses. We observe that there is uncertainty in the magnitude of the frequency peaks, since they correspond to the different natural frequencies of the oscillators in the model. We consider this uncertainty as multiplicative model error since the natural frequency of the disturbance $\left({\varpi }_l\right)$ are “easy” to obtain (e.g., using the power spectral density), but the damping coefficients $\left({\zeta }_l\right)$ and the gains $\left({\gamma }_l\right)$ are not necessarily easy to obtain.

For a more detailed analysis, in Figure 6 we present the true frequency response of the true disturbance discrete-time transfer function for this example and the estimated discrete-time transfer function using the method in [35] (with $\vartheta =1$). The true system discrete-time transfer function was obtained from the continuous-time equation in Equation (11), the parameters in Table 1, and utilising Algorithm 1. We observe that, for some frequencies, the phase exhibits a very different behaviour when compared to the true frequency response, especially at low frequencies.

On the other hand, the performance of MVC can be obtained using the sensitivity function ($\mathbb{S}$) [23,42,48]. The sensitivity function represents the closed-loop transfer function between the external perturbations and system output. The sensitivity function of the AO system in Figure 4 is defined as:

$$\mathbb{S}\left(z^{-1}\right)=\frac{1}{1+K\left(z^{-1}\right)G\left(z^{-1}\right)},$$

(63)

thus, when considering the disturbance model with uncertainty $\widehat{H}\left(z^{-1}\right)=\widehat{F}\left(z^{-1}\right)+z^{-d}\widehat{R}\left(z^{-1}\right)$, we obtain

$$\mathbb{S}\left(z^{-1}\right)=\frac{\widehat{F}\left(z^{-1}\right)}{\widehat{F}\left(z^{-1}\right)+z^{-d}\widehat{R}\left(z^{-1}\right)}=\frac{\widehat{F}\left(z^{-1}\right)}{\widehat{H}\left(z^{-1}\right)}.$$

(64)

If the true $H\left(z^{-1}\right)$ was known $(\widehat{H}\left(z^{-1}\right)=H\left(z^{-1}\right))$, the output of the system can be obtained from the definition of the sensitivity as [48]:

$$\begin{array}{c}\hfill {\widehat{y}}_k=\mathbb{S}\left(z^{-1}\right)H\left(z^{-1}\right)e_k=F\left(z^{-1}\right)e_k.\end{array}$$

(65)

From the last expression, the variance of the output signal is given by:

$${\mathrm{var}}_{y_0}=\mathrm{var}\left\{y_k\right\}=\left(1+f_1^2+f_2^2+\cdots +f_{d-1}^2\right){\lambda }^2.$$

(66)

However, since we only have an estimation of $H\left(z^{-1}\right)$, defined as $\widehat{H}\left(z^{-1}\right)$, the expression for the output signal is given by:

$$\begin{array}{c}\hfill {\widehat{y}}_k=\mathbb{S}\left(z^{-1}\right)H\left(z^{-1}\right)e_k=\frac{\widehat{F}\left(z^{-1}\right)H\left(z^{-1}\right)}{\widehat{H}\left(z^{-1}\right)}{e}_k.\end{array}$$

(67)

Notice that the model mismatch between the true and estimated transfer functions has the effect of increasing the output signal variance, Therefore, the identification of the AO system is an important issue in the attainment of good control performance, which, in turn, allows for improving the quality of astronomical images.

6.3. Control Performance under Model Mismatch

In order to evaluate the impact of model accuracy on the control performance, we define the following coefficient as a performance metric:

$$E=\left(\frac{{\mathrm{var}}_{\widehat{y}}}{{\mathrm{var}}_{y_0}}-1\right)100\%,$$

(68)

where ${\mathrm{var}}_{\widehat{y}}$ is the output variance using the estimated disturbance transfer function $\widehat{H}\left(z^{-1}\right)$ and ${\mathrm{var}}_{y_0}$ is output variance obtain when the disturbance model matches the true disturbance, $\widehat{H}\left(z^{-1}\right)=H\left(z^{-1}\right)$.

Finally, in order to illustrate the performance of MVC, we use a mirror gain $M_0=1$ and WFS gain $D_0=1$. We also select a measurement delay $\mu =1$ and correction delay $\tau =1$, i.e., $d=2$. Subsequently, using the true disturbance transfer function $\left(H\left(z^{-1}\right)\right)$ and the estimated disturbance transfer function $\widehat{H}\left(z^{-1}\right)$ we can obtain the sensitivity function and the output variance to evaluate the performance of MVC.

Table 2. Performance of MVC for different disturbance model.

$\mathit{\vartheta }$	0	0.5
$\mathit{E}(\%)$	0	3.85

shows the performance (E) for different values of $\vartheta$, where it is clear that $\vartheta =0$ is the ideal case, i.e., the estimated model without uncertainty. The $E(\%)$ denotes the percentage variation of the output variance with respect to the ideal case. We observe that the model accuracy has a big impact on the performance of MVC.

7. Numerical Example

In this example, we consider the six damped oscillators that are shown in Table 1 (which corresponds to the numerical example presented in [26]), and we consider that this is the true system. Figure 7 shows the discrete-time PSD of the disturbances, where we can clearly observe the disturbance frequencies. Note that this example has a similar shape and structure as the ones obtained during the 2014b observing run by the University of Arizona MagAO Team on the night of 31 October 2014, at the Clay Telescope in Las Campanas Observatory [26], which we include in Figure 8 for the completeness of the presentation. However, the PSD in Figure 8 corresponds to the data that were obtained in closed-loop, while the simulated data in Figure 7 are in an open-loop. Because of this fact, both of the PSDs have the same resonance peaks, but different magnitudes.

7.1. Disturbance Identification

We identify the continuous-time system parameters from a set of simulated data (generated with the six damped oscillators shown in Table 1) using the proposed system model in Equation (46) with the Whittle’s likelihood identification technique, and the system model in Equation (13) with the NLS fitting method proposed in [35]. We utilise a sampling period of $\Delta =5$ ms, 100 Monte-Carlo (MC) simulations, and two data lengths, namely $N=500$ and $N=1000$. The estimation problems in Equations (53) and (56) are solved with a local optimisation algorithm using the Matlab^® function fmincon, setting the interior-point algorithm, the continuous-time, and discrete-time constraints.

Figure 9 shows the estimation results from the MC simulations of the system with six oscillators for a) Whittle’s likelihood method and b) NLS fitting method. Each identification approach yields a different $\left(H_{\mathcal{E}}\left(z^{-1}\right)\right)$. In this Figure, the solid blue line represents the average of all the periodograms, the dashed black line represents the true discrete-time PSD, the dotted red line represents the average of all the estimated discrete-time PSDs, and the gray shaded region represents the area in which all of the estimated spectra lie. Despite that Whittle’s likelihood is an asymptotic approximation $(N\to \infty )$ to the traditional likelihood function, in Figure 9 we can observe that the estimation with this approach exhibits good accuracy for small data length, namely $N=500$. Moreover, we observe that, the greater the N, the smaller the gray-shaded region. In contrast, the NLS fitting method does not exhibit an accurate estimation.

Table 3. Estimated parameters of six disturbance for Whittle’s likelihood method, using $N=1000$.

	$\widehat{\mathit{\alpha }}$ (Hz)	$\widehat{\mathit{\zeta }}$	$\widehat{\mathit{\beta }}(\times {10}^4)$
1	$2.143\pm 0.320$	$0.861\pm 0.152$	$1.559\pm 0.331$
2	$14.024\pm 0.194$	$0.059\pm 0.014$	$3.775\pm 0.316$
3	$20.026\pm 0.296$	$0.055\pm 0.012$	$3.807\pm 0.287$
4	$29.010\pm 0.385$	$0.054\pm 0.010$	$4.307\pm 0.301$
5	$43.020\pm 0.473$	$0.050\pm 0.009$	$5.503\pm 0.305$
6	$60.085\pm 0.411$	$0.050\pm 0.006$	$9.678\pm 0.416$

and

Table 4. Estimated parameters of six disturbance for NLS fitting method, using $N=1000$.

	$\mathit{\alpha }$ (Hz)	$\mathit{\zeta }$	${\widehat{\mathit{\gamma }}}_{\mathit{d}}$
1	$0.704\pm 0.424$	$0.999\pm 0.009$	$3.665\pm 0.697$
2	$13.964\pm 0.473$	$0.077\pm 0.024$	$14.567\pm 1.595$
3	$19.930\pm 0.540$	$0.088\pm 0.097$	$16.340\pm 7.776$
4	$28.878\pm 0.813$	$0.108\pm 0.054$	$19.363\pm 5.393$
5	$42.912\pm 0.969$	$0.142\pm 0.053$	$24.132\pm 4.714$
6	$59.794\pm 1.082$	$0.118\pm 0.033$	$23.031\pm 3.630$

show the average and standard deviation of all the estimated parameters. The true values are presented in Table 1. It is clear that, for the proposed method, the estimates are very accurate, exhibiting a mean value that is almost identical to the true values and small standard deviations.

In addition, the dashed-red line frequency response plot shown in Figure 6 corresponds to the estimated model of the disturbances that were obtained with the average of all the estimated parameters using the NLS fitting method in this example.

7.2. Performance of MVC in AO System

In order to illustrate the performance of MVC, we consider a mirror gain $M_0=1$ and WFS gain $D_0=1$. We also select a measurement delay $\mu =1$ and correction delay $\tau =1$, i.e. $d=2$. The MVC algorithm used in this section is summarised in Algorithm 3.

Algorithm 3 MVC algorithm

1: Compute the disturbance continuous-time model using the Algorithm 2 or other identification method 2: Obtain the discrete-time model of disturbance in state-space consider that the sampling period can be different for control 3: Compute the spectral factorization 4: Compute the coefficients of the finite impulse response filter $F\left(z^{-1}\right)$ using Equation (59) 5: Compute the transfer function $R\left(z^{-1}\right)$ using Equation (58) 6: Compute the controller $K\left(z^{-1}\right)$ using Equation (61) 7: Compute de output variance ${\mathrm{var}}_y$ using Equation (66)

On the other hand,

Table 5. Performance of MVC $\left(E\right)$ for different identifications techniques of disturbance, and different values of $\vartheta$. We use $\Delta =5$ ms, and $N=1000$.

	$\widehat{\mathit{H}}\left({\mathit{z}}^{-1}\right)$
$\mathbf{\vartheta }$	Proposed Method	NLS [35]
0	0	0
$0.5$	$1.56\%$	$3.85\%$
1	$0.870\%$	$16.52\%$
$1.5$	$2.29\%$	$44.79\%$

shows the performance of MVC (E), for $\Delta =5$ ms, when using the estimated model that was obtained from both our proposal and the NLS method from [35]. For illustrative purposes, we modified the value of $\vartheta$, from 1 to $0.5$ and $1.5$, in order to verify how the performance of MVC varies for different values of $\vartheta$. The results show that our method outperforms [35] in every case. In addition, the performance of the MVC does not change significantly when designing the MVC using the estimated model from our proposal for different values of $\vartheta$. The same could not be said about the method in [35]. From these results, we can conclude that our proposal yields more accurate estimates that, in turn, can result in more effective controllers. Note that the performance of MVC for NLS fitting method is the same as that shown in Table 2, i.e., using typical AR(2) model.

In Figure 10, we show the discrete-time PSD of the controlled output using the true model $S^{y_0}\left(e^{\mathrm{i}\omega \Delta }\right)$ that is represented by the solid orange line with square marks, the model obtained by Whittle’s likelihood method $\left(S^y\left(e^{\mathrm{i}\omega \Delta }\right)\right)$ represented by the dashed blue line, and the discrete-time PSD of the controlled output using the AR(2) model that is obtained by the NLS fitting method $\left({\overline{S}}^y\left(e^{\mathrm{i}\omega \Delta }\right)\right)$ represented by the dotted black line. Note that these controlled outputs are obtained in closed-loop. We also present the discrete-time PSD of the disturbance signal $\left(S^{{\phi }^{\mathrm{Tot}}}\left(e^{\mathrm{i}\omega \Delta }\right)\right)$ in open-loop represented by the densely dotted red line, corresponding to the input signal. We observe that the controlled output using the proposed modelling and identification methods exhibits better performance than using the NLS method. It is clear that the controlled output obtained with our proposal and the true controlled output are almost identical. To see the latter, we compute the normalized mean square error (NMSE), as follows:

$$\mathrm{NMSE}=\frac{\sqrt{{\sum }_{k=1}^N{\left(|S^{y_0}\left(e^{\mathrm{i}{\omega }_k\Delta }\right)|-|S^y\left(e^{\mathrm{i}{\omega }_k\Delta }\right)|\right)}^2}}{\sqrt{{\sum }_{k=1}^N{\left(|S^{y_0}\left(e^{\mathrm{i}{\omega }_k\Delta }\right)|\right)}^2}}=0.0067,$$

(69)

where $|·|$ denotes the magnitude. We observe that the NMSE value is very small, indicating that both outputs are almost the same and corroborating the benefits of our identification technique on MVC performance.

Note that the controlled output using the estimated AR(2) model and NLS fitting method do not exhibit the disturbance peaks for low frequencies, but, for high frequencies, the disturbances are amplified. This is a consequence of the so-called waterbed effect [48], where the magnitude of the frequency response is reduced in one part of the spectrum, but it is increased in other parts of the spectrum.

Finally we evaluate the effect of the accuracy the estimated model in the control performance when we use a different sampling period, namely $\Delta =10$ ms. The results are the following: $E=2.83\%$ for the proposed method and $E=18.08\%$ for the NLS. It is clear that the control performance does not exhibit a great variation with respect to the results with $\Delta =5$ ms. Moreover, the new estimates were directly obtained by replacing $\Delta =0.010$ in Equation (46), without the need of running new experiments. In contrast, the approach that is presented in [35] requires carrying out the identification method again with a new sampling period and a new set of samples to obtain the new estimates and then the performance of MVC.

8. Discussion and Conclusions

In this paper we addressed the problem of disturbance modelling for minimum variance control in AO systems. We analysed the impact of the accuracy of the disturbance model in the control performance. We showed that the model mismatch (from the parameter estimation) has an impact on the control performance when we used a typical auto-regressive discrete-time disturbance model. We presented a discrete-time model (arising from a continuous-time autoregressive model) for the disturbance that is more accurate at high frequencies.

We obtained a sampled-data state-space model that is more compact and of smaller order than the one that is presented in [30]. We show that the model of the disturbance and the plant can be obtained separately, thus simplifying the design of the controller.

From our simulations we can conclude that Whittle’s likelihood method allows for the identification of continuous-time parameters by directly optimising an explicit function of them. The corresponding cost function is similar to the typical cost function that is obtained from the AR(2) model. However, our proposal yields more accurate estimates whilst maintaining a low computational cost.

We showed that the output of the AO closed-loop system, when using our modelling and estimation proposal, has the best performance, eliminating all of the disturbance peaks for the frequency range of interest. We observed that, when we used the MVC with the correct identification, we improved the performance of controller in $19\%$ with respect to the MVC with the estimated model utilising the classical estimation technique shown in [35].