Open AccessArticle

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

Departamento Electrónica, Universidad Técnica Federico Santa María (UTFSM), Av. España 1680, Valparaíso 2390123, Chile

Advanced Center for Electrical and Electronic Engiennering, AC3E, Av. Matta 222, Valparaíso 2580129, Chile

Departamento de Ingeniería Electrica, Facultad de Ingeniería, Universidad de Los Andes, Av. Alberto Carnevali, Mérida 5101, Venezuela

⁴

Instituto de Electricidad y Electrónica, Facultad de Ciencias de la Ingeniería, Universidad Austral de Chile (UACH), Genaral Lagos 2086, Valdivia 5111187, Chile

Author to whom correspondence should be addressed.

Sensors 2021, 21(9), 3054; https://doi.org/10.3390/s21093054

Submission received: 22 March 2021 / Revised: 21 April 2021 / Accepted: 24 April 2021 / Published: 27 April 2021

(This article belongs to the Section Optical Sensors)

Download

Browse Figures

Versions Notes

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.

Keywords:

adaptive optics; wavefront sensor; disturbances; modelling; identification; minimum variance controller; Whittle’s likelihood

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

(M_{c} (s))

, a controller

(K (z^{- 1}))

, and a zero order hold (ZOH) [6,21,22], where

s = \frac{d}{d t}

should be understood as the derivative operator or the argument of the Laplace transform, and z is either the forward shift operator

(z φ_{k} = φ_{k + 1})

or the argument of the Z-transform. The continuous-time signals

φ^{Tot} (t)

φ^{cor} (t)

and

φ^{res} (t)

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

φ^{Tot} (t)

is defined as the sum of the atmospheric turbulence effect

φ^{tur} (t)

and the different vibration sources effect

φ^{{vib}_{l}} (t)

l = 1, \dots, m

[5,6,22]. That is:

φ^{Tot} (t) = φ^{tur} (t) + φ^{{vib}_{1}} (t) + φ^{{vib}_{2}} (t) + \dots + φ^{{vib}_{m}} (t),

(1)

where m is the number of vibration sources. The signal

η_{k}

shown in Figure 1 is a discrete-time additive zero-mean white Gaussian noise with variance

σ_{η}^{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

[4 F_{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

φ^{res} (t)

during a time interval

[t_{k - 1}, t_{k}]

, where

t_{k}

is a sampling time, accumulating photons over a time interval

Δ

[22]. Thus, the discrete-time residual phase is given by:

φ_{k}^{res} = \frac{1}{Δ} \int_{t_{k - 1}}^{t_{k}} φ^{res} (t) d t; t_{k} = k Δ,

(2)

where

k \in N_{0}

and

Δ

represent the sampling period. The corresponding discrete-time transfer function of the WFS in Equation (2) is given by [22]:

D (z^{- 1}) = D_{0} z^{- μ},

(3)

where

D_{0}

is the gain and

μ

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 (z^{- 1}) = M_{0} z^{- τ},

(4)

where

M_{0}

is the gain and

τ

is the correction delay.

On the other hand, we consider that all of the elements in Figure 1 are linear. Subsequently,

φ_{k}^{res}

is a linear function of the past values of

η_{k}

. The output of the system presented in Figure 1 is given by

\begin{matrix} y_{k} = D_{0} φ_{k - μ}^{res} + η_{k}, \end{matrix}

(5)

where

φ_{k}^{res}

is the discrete-time residual phase that is given by:

φ_{k}^{res} = φ_{k}^{Tot} - φ_{k}^{cor},

(6)

and

η_{k}

is a zero-mean white Gaussian noise. Subsequently, the variance of the output

y_{k}

is given by:

var \{y_{k}\} = σ_{η}^{2} + D_{0}^{2} var \{φ_{k - μ}^{res}\},

(7)

thus, minimising the variance of

y_{k}

is the same as minimising the variance of

φ_{k}^{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 (z^{- 1})

is the discrete-time transfer function of the controller,

M (z^{- 1})

is the DM discrete-time transfer function,

D (z^{- 1})

is the WFS discrete-time transfer function,

u_{k}

is the output of controller, and

y_{k}

is the system output. Note that

G (z^{- 1})

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} φ_{k - μ}^{cor} + χ_{k},

(8)

where

χ_{k} = η_{k} + D_{0} φ_{k - μ}^{Tot},

(9)

corresponds to the disturbances of the AO system. Subsequently,

y_{k} = D_{0} [φ_{k - μ}^{Tot} - φ_{k - μ}^{cor}] + η_{k} .

(10)

Note that

φ_{k}^{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{matrix} φ^{Tot} (t) & = \sum_{l = 0}^{m} \frac{β_{l}}{s^{2} + 4 ζ_{l} π α_{l} s + {(2 π α_{l})}^{2}} {\overset{\cdot}{ν}}_{l} (t), \\ = \sum_{l = 0}^{m} \frac{β_{l}}{s^{2} + 2 ζ_{l} ϖ_{l} s + ϖ_{l}^{2}} {\overset{\cdot}{ν}}_{l} (t), \end{matrix}

(11)

where

β_{l}

is the gain,

α_{l}

(Hz) is the natural frequency,

ζ_{l}

is the damping coefficient, and

{\overset{\cdot}{ν}}_{l} (t)

is a continuous-time zero-mean white Gaussian noise with variance

σ_{l}^{2} = 1 δ_{D} (t), \forall l

(

δ_{D} (t)

is the Dirac delta). We assume that the noises

{\overset{\cdot}{ν}}_{l} (t), l = 0, \dots, m

are jointly uncorrelated. We use

(ϖ_{l} = 2 π α_{l})

for simplicity in the presentation.

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

φ_{k}^{Tot} = \sum_{l = 0}^{m} H_{l} (z^{- 1}) ε_{l, k},

(12)

where

ε_{l, k}

is a zero-mean white Gaussian noise with variance

σ_{ε_{l, k}}^{2} = 1, \forall l

, jointly independent,

H_{0} (z^{- 1})

is the discrete-time transfer function that is associated with the turbulence term and

H_{l} (z^{- 1})

l = 1, \dots, m

, and the transfer function associated with the vibrations terms. The discrete-time transfer function

H_{l} (z^{- 1})

is modelled using the following approximated auto-regressive (AR(2)) discrete-time representation [5,6,7,11,26,34,35]:

H_{l} (z^{- 1}) = \frac{γ_{l}}{1 - a_{1_{l}} z^{- 1} + a_{2_{l}} z^{- 2}},

(13)

where

γ_{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}} = 2 e^{- ζ_{l} ϖ_{l} Δ} \cos (ϖ_{l} \sqrt{1 - ζ_{l}^{2}} Δ),

(14)

a_{2_{l}} = - e^{- 2 ζ_{l} ϖ_{l} Δ} .

(15)

Remark 1.

In the approximated model, the mathematical expression for the gains

γ_{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}} = 2 e^{- ζ_{0} ϖ_{0} T} \cos (ϖ_{0} \sqrt{1 - ζ_{0}^{2}} T),

(16)

a_{2_{0}} = - e^{- 2 ζ_{0} ϖ_{0} T},

(17)

ϖ_{0} = 0.6 π (n (i) + 1) \frac{V_{0}}{D},

(18)

where

n (i)

is the radial order of the Zernike mode i,

T

is the exposure time of the system,

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:

\overset{\cdot}{x} (t) = A_{c} x (t) + κ_{c} Ψ_{c} (t),

(19)

φ^{Tot} (t) = \overset{\cdot}{b} (t) = C_{c} x (t),

(20)

where

A_{c} \in R^{2 m \times 2 m}

x (t) \in R^{2 m \times 1}

κ_{c} \in R^{2 m \times m}

Ψ_{c} \in R^{m \times 1}

C_{c} \in R^{1 \times 2 m}

, m is the number of damping oscillators used, and the noise variance is

Q_{c} = κ_{c} κ_{c}^{T}

. These continuous-time system matrices are given by:

A_{c} = [\begin{matrix} 0 & 1 & 0 & \dots & 0 \\ - ϖ_{1}^{2} & - 2 ζ_{1} ϖ_{1} & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ 0 & 0 & \dots & 0 & 1 \\ 0 & 0 & \dots & - ϖ_{m}^{2} & - 2 ζ_{m} ϖ_{m} \end{matrix}],

(21)

κ_{c} = [\begin{matrix} 0 & \dots & \dots & \dots & 0 \\ β_{1} & 0 & \dots & \dots & ⋮ \\ 0 & 0 & ⋱ & ⋱ & ⋮ \\ ⋮ & β_{2} & ⋱ & ⋱ & ⋮ \\ ⋮ & 0 & ⋱ & ⋱ & ⋮ \\ ⋮ & ⋮ & 0 & ⋱ & ⋮ \\ ⋮ & ⋮ & ⋮ & ⋱ & 0 \\ 0 & 0 & 0 & 0 & β_{m} \end{matrix}], Ψ_{c} (t) = [\begin{matrix} {\overset{\cdot}{ω}}_{1} (t) \\ {\overset{\cdot}{ω}}_{2} (t) \\ ⋮ \\ {\overset{\cdot}{ω}}_{m} (t) \end{matrix}],

(22)

C_{c} = [1 0 | 1 0 | \dots | 1 0],

(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

\overset{\cdot}{b} (t)

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):

φ_{k}^{res} = \frac{1}{Δ} \int_{t_{k - 1}}^{t_{k}} φ^{res} (t) d t,

(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

X (t) = {[x^{T} (t) b (t)]}^{T}

[41]. Subsequently, the continuous-time state-space model is as follows:

\overset{\cdot}{X} (t) = A X (t) + [\begin{matrix} κ_{c} \\ 0 \end{matrix}] Ψ_{c} (t),

(25)

where

A = [\begin{matrix} A_{c} & 0 \\ C_{c} & 0 \end{matrix}] .

(26)

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

X_{k + 1} = e^{A Δ} X_{k} + n_{k + 1},

(27)

if we consider that

A_{c}

is invertible, we have the following result [43]:

e^{A Δ} = [\begin{matrix} e^{A_{c} Δ} & 0 \\ C_{c} A_{c}^{- 1} (e^{A_{c} Δ} - I) & I \end{matrix}],

(28)

where

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}^{Δ} e^{\begin{matrix} [\begin{matrix} A_{c} & 0 \\ C_{c} & 0 \end{matrix}] ξ \end{matrix}} [\begin{matrix} κ_{c} \\ 0 \end{matrix}] Ψ_{k + 1 - η} d ξ .

(29)

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

e^{\begin{matrix} [\begin{matrix} - A & Q_{c}^{*} \\ 0 & A^{T} \end{matrix}] Δ \end{matrix}} = [\begin{matrix} P_{11} & P_{12} \\ 0 & P_{22} \end{matrix}],

(30)

where:

Q_{c}^{*} = [\begin{matrix} κ_{c} \\ 0 \end{matrix}] {[\begin{matrix} κ_{c} \\ 0 \end{matrix}]}^{T},

(31)

then,

\begin{matrix} e^{A Δ} & = P_{22}^{T}, \end{matrix}

(32)

\begin{matrix} Q & = \int_{0}^{Δ} e^{A ξ} Q_{c}^{*} e^{{A ξ}^{T}} d ξ = P_{22}^{T} P_{12}, \\ = [\begin{matrix} Q_{11} & Q_{12} \\ Q_{12}^{T} & Q_{22} \end{matrix}], \end{matrix}

(33)

where

Q_{11} \in R^{2 m \times 2 m}

Q_{12} \in R^{1 \times 2 m}

and

Q_{22} \in R^{1 \times 1}

Thus, we have that

[\begin{matrix} x_{k + 1} \\ b_{k + 1} \end{matrix}] = [\begin{matrix} e^{A_{c} Δ} & 0 \\ C_{c} A_{c}^{- 1} (e^{A_{c} Δ} - I) & I \end{matrix}] [\begin{matrix} x_{k} \\ b_{k} \end{matrix}] + n (t_{k + 1}),

(34)

\begin{matrix} x_{k + 1} = A x_{k} + w_{k + 1}, \\ φ_{k + 1}^{Tot} Δ = b_{k + 1} - b_{k} = C x_{k} + v_{k + 1}, \end{matrix}

(35)

where:

\begin{matrix} A & = e^{A_{c} Δ}, \end{matrix}

(36)

\begin{matrix} C & = C_{c} A_{c}^{- 1} {A - I}, \end{matrix}

(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 (z^{- 1}) = M (z^{- 1}) D (z^{- 1}))

in state-space form as:

\begin{matrix} {\overset{\cdot}{x}}_{p} (t) & = A_{p_{c}} x (t) + B_{p_{c}}, \end{matrix}

(38)

\begin{matrix} {\overset{\cdot}{b}}_{p} (t) & = C_{p_{c}} x_{p} (t), \end{matrix}

(39)

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

\begin{matrix} x_{p_{k + 1}} & = A_{p} x_{p_{k}} + B_{p} u_{k}, \end{matrix}

(40)

\begin{matrix} \frac{b_{p_{k + 1}} - b_{p_{k}}}{Δ} & = C_{p} x_{p_{k}} + D_{p} u_{k}, \end{matrix}

(41)

where:

\begin{matrix} A_{p} & = e^{A_{p_{c}} Δ}, \end{matrix}

(42)

\begin{matrix} B_{p} & = \int_{0}^{Δ} e^{A_{p_{c}} ξ} B_{p_{c}} d ξ, \end{matrix}

(43)

\begin{matrix} C_{p} & = \frac{1}{Δ} \int_{0}^{Δ} e^{A_{p_{c}} ξ} d ξ, \end{matrix}

(44)

\begin{matrix} D_{p} & = \frac{1}{Δ} C_{p_{c}} \int_{0}^{Δ} \int_{0}^{ξ} e^{A_{p_{c}} ϕ} d ϕ d ξ . \end{matrix}

(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 (z^{- 1})

[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)):

φ_{k}^{Tot} = H (z^{- 1}) e_{k},

(46)

where

e_{k}

is a discrete-time zero-mean white Gaussian noise with variance

λ^{2}

. Note that both

H (z^{- 1})

and

λ^{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 (e^{i ω_{j} Δ}) = {| H (e^{i ω_{j} Δ}) |}^{2} λ^{2},

(47)

where

i = \sqrt{- 1}

j = 1, 2, \dots, N

and

0 < ω_{j} \leq 2 π (F_{s} / 2)

, and

ω_{j} = (2 π j) / Δ

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 (z^{- 1})$ and $λ^{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

(ε)

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

σ_{ε}^{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,

{\vec{θ}}_{N L S}

, as follows:

{\vec{θ}}_{N L S} = {[\begin{matrix} {\vec{ζ}}^{T} & {\vec{γ}}^{T} \end{matrix}]}^{T},

(48)

where

\vec{ζ} = {[ζ_{1} ζ_{2} \dots ζ_{m}]}^{T}

is the vector that contains all of the damping coefficients, and

\vec{γ} = {[γ_{1} γ_{2} \dots γ_{m}]}^{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

{\bar{S}}_{l} (e^{i ω_{j} Δ})

in Equation (49). The different frequencies (Hz) in

\vec{α} = {[α_{1} α_{2} \dots α_{m}]}^{T}

were found in the periodogram as the frequency peaks.

The discrete-time PSD is given by:

{\bar{S}}_{l} (e^{i ω_{j} Δ}) = {|\frac{γ_{l}}{1 - a_{1_{l}} e^{- i ω_{j} Δ} - a_{2_{l}} e^{- 2 i ω_{j} Δ}}|}^{2},

(49)

and the NLS cost function

F_{l} (θ_{N L S})

is obtained from

ϵ_{l} (e^{i ω_{j} Δ}) = 10 \log [\frac{I_{l} (e^{i ω_{j} Δ})}{{\bar{S}}_{l} (e^{i ω_{j} Δ})}],

(50)

where

I (e^{i ω_{j} Δ})

is the periodogram of the sampled-data series

φ^{Tot}

, and it is given by:

I (e^{i ω_{j} Δ}) = {|\frac{1}{\sqrt{N}} \sum_{t = 1}^{N} φ_{t}^{Tot} e^{- i ω_{j} Δ}|}^{2},

(51)

then,

F_{l} ({\vec{θ}}_{N L S}) = \frac{1}{2} ϵ_{l} {(e^{i ω_{j} Δ})}^{T} ϵ_{l} (e^{i ω_{j} Δ}) .

(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:

{\hat{θ}}_{N L S_{l}} = \arg \min_{{\vec{θ}}_{N L S_{l}}} (F_{l} ({\vec{θ}}_{N L S})), s . t . \{\begin{matrix} 0 < ζ_{l} < 1 \\ γ_{l} > 0 \end{matrix}, l = 0, \dots, 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:

ℓ_{N} (\vec{θ}) = - \frac{1}{N} \sum_{j = 1}^{N} \log S (e^{i ω_{j} Δ}) - \frac{1}{N} \sum_{j = 1}^{N} \frac{I (e^{i ω_{j} Δ})}{S (e^{i ω_{j} Δ})},

(54)

where N is the data length and

\vec{θ}

is the vector of parameters to be estimated that is given by:

\vec{θ} = {[\begin{matrix} {\vec{α}}^{T} & {\vec{ζ}}^{T} & {\vec{β}}^{T} \end{matrix}]}^{T},

(55)

\vec{α} = {[α_{1} α_{2} \dots α_{m}]}^{T}

is the vector that contain all the natural frequencies (in Hz),

\vec{ζ} = {[ζ_{1} ζ_{2} \dots ζ_{m}]}^{T}

is the vector that contains all the damping coefficients, and

\vec{β} = {[β_{1} β_{2} \dots β_{m}]}^{T}

is the vector that contains all the gains. In Equation (54)

S (e^{i ω_{j} Δ})

is the discrete-time PSD of the disturbances given in Equation (47) and

I (e^{i ω_{j} Δ})

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)

\vec{θ}

, since the spectrum

S (e^{i ω_{j} Δ})

depends on them.

Afterwards, the ML estimation problem is given by:

{\hat{θ}}_{ML} = \arg \min_{\vec{θ}} - ℓ_{N} (\vec{θ}), s . t . \{\begin{matrix} 0 < \vec{α} < F_{s} / 2 \\ 0 < \vec{ζ} < 1 \\ \vec{β} > 0 \end{matrix} .

(56)

Finally, we summarize our identification technique in Algorithm 2.

Algorithm 2 Identification algorithm

1:: Compute the periodogram $I (e^{i ω_{j} Δ})$ using Equation (51) with $ω_{j} = (2 π j) / Δ$ , and $j = 1, 2, \dots, N$
2:: Choose an initial guess ${\hat{θ}}^{(0)}$
3:: Obtain the discrete-time spectrum $S (e^{i ω Δ})$ in function of the vector of continuous-time parameters $\vec{θ}$ 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,

var \{φ^{res} (t)\}

[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,

var \{φ_{k}^{res}\}

The goal of MVC is to minimise the variance of the output signal

var \{y_{k}\}

[42]. For simplicity of the presentation, we consider that the measurement noise,

η_{k}

, is an independent and identically distributed (iid) signal. We also assume that the atmospheric turbulence

φ^{tur}

, the vibrations

φ^{vib}

, and the measurement noise

η

, are mutually independent. In this setup, the measurement noise

η_{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 (z^{- 1}) e_{k},

(57)

where

d = μ + τ

and

u_{k}

is the output of the controller.

H (z^{- 1})

is the discrete-time transfer function obtained from the spectral factorization [46] of the spectrum of the signal

χ_{k}

in Equation (9). The signal

e_{k}

is a zero-mean white noise with variance

λ^{2}

that is also obtained from the spectral factorization procedure.

In order to obtain the MVC [42,52] we rewrite

H (z^{- 1})

H (z^{- 1}) = F (z^{- 1}) + z^{- d} R (z^{- 1}),

(58)

where

F (z^{- 1})

is a finite impulse response (FIR) filter that corresponds to the first d samples of the impulse response of

H (z^{- 1})

, given by:

F (z^{- 1}) = 1 + f_{1} z^{- 1} + f_{2} z^{- 2} + \dots + f_{d - 1} z^{- (d - 1)},

(59)

where

f_{ι}

is the

ι

th coefficient of the FIR filter

F (z^{- 1})

, and

R (z^{- 1})

is a causal filter that can be obtained from Equation (58):

R (z^{- 1}) = (H (z^{- 1}) - F (z^{- 1})) z^{d} .

(60)

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

\begin{matrix} K (z^{- 1}) = \frac{R (z^{- 1})}{F (z^{- 1}) D_{0} M_{0}}, \end{matrix}

(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{matrix} \hat{H} (z^{- 1}) = H (z^{- 1}) (1 + H_{E} (z^{- 1}) ϑ), \end{matrix}

(62)

where

H (z^{- 1})

is the true discrete-time transfer function and

H_{E} (z^{- 1}) ϑ

is the multiplicative model error. The subscript

E

indicates multiplicative errors, whilst the parameter

ϑ

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

ϑ = 1

in order to obtain and subsequently analyse

H_{E} (z^{- 1})

. In some cases, a different set-up (e.g., a different measurement noise variance) could lead to the same

H_{E} (z^{- 1})

, but with a different gain.

In Table 1, 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

χ_{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

ϑ = 1

, we obtain the multiplicative model error

(H_{E} (z^{- 1}))

using Equation (62). Finally, we select different values of

ϑ

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

(ϖ_{l})

are “easy” to obtain (e.g., using the power spectral density), but the damping coefficients

(ζ_{l})

and the gains

(γ_{l})

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

ϑ = 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 (

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:

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

(63)

thus, when considering the disturbance model with uncertainty

\hat{H} (z^{- 1}) = \hat{F} (z^{- 1}) + z^{- d} \hat{R} (z^{- 1})

, we obtain

S (z^{- 1}) = \frac{\hat{F} (z^{- 1})}{\hat{F} (z^{- 1}) + z^{- d} \hat{R} (z^{- 1})} = \frac{\hat{F} (z^{- 1})}{\hat{H} (z^{- 1})} .

(64)

If the true

H (z^{- 1})

was known

(\hat{H} (z^{- 1}) = H (z^{- 1}))

, the output of the system can be obtained from the definition of the sensitivity as [48]:

\begin{matrix} {\hat{y}}_{k} = S (z^{- 1}) H (z^{- 1}) e_{k} = F (z^{- 1}) e_{k} . \end{matrix}

(65)

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

{var}_{y_{0}} = var \{y_{k}\} = (1 + f_{1}^{2} + f_{2}^{2} + \dots + f_{d - 1}^{2}) λ^{2} .

(66)

However, since we only have an estimation of

H (z^{- 1})

, defined as

\hat{H} (z^{- 1})

, the expression for the output signal is given by:

\begin{matrix} {\hat{y}}_{k} = S (z^{- 1}) H (z^{- 1}) e_{k} = \frac{\hat{F} (z^{- 1}) H (z^{- 1})}{\hat{H} (z^{- 1})} e_{k} . \end{matrix}

(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 = (\frac{{var}_{\hat{y}}}{{var}_{y_{0}}} - 1) 100 %,

(68)

where

{var}_{\hat{y}}

is the output variance using the estimated disturbance transfer function

\hat{H} (z^{- 1})

and

{var}_{y_{0}}

is output variance obtain when the disturbance model matches the true disturbance,

\hat{H} (z^{- 1}) = H (z^{- 1})

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

μ = 1

and correction delay

τ = 1

, i.e.,

d = 2

. Subsequently, using the true disturbance transfer function

(H (z^{- 1}))

and the estimated disturbance transfer function

\hat{H} (z^{- 1})

we can obtain the sensitivity function and the output variance to evaluate the performance of MVC. Table 2 shows the performance (E) for different values of

ϑ

, where it is clear that

ϑ = 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

Δ = 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

(H_{E} (z^{- 1}))

. 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 and Table 4 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

μ = 1

and correction delay

τ = 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 (z^{- 1})$ using Equation (59)
5:: Compute the transfer function $R (z^{- 1})$ using Equation (58)
6:: Compute the controller $K (z^{- 1})$ using Equation (61)
7:: Compute de output variance ${var}_{y}$ using Equation (66)

On the other hand, Table 5 shows the performance of MVC (E), for

Δ = 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

ϑ

, from 1 to

0.5

and

1.5

, in order to verify how the performance of MVC varies for different values of

ϑ

. 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

ϑ

. 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}} (e^{i ω Δ})

that is represented by the solid orange line with square marks, the model obtained by Whittle’s likelihood method

(S^{y} (e^{i ω Δ}))

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

({\bar{S}}^{y} (e^{i ω Δ}))

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

(S^{φ^{Tot}} (e^{i ω Δ}))

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:

NMSE = \frac{\sqrt{\sum_{k = 1}^{N} {(| S^{y_{0}} (e^{i ω_{k} Δ}) | - | S^{y} (e^{i ω_{k} Δ}) |)}^{2}}}{\sqrt{\sum_{k = 1}^{N} {(| S^{y_{0}} (e^{i ω_{k} Δ}) |)}^{2}}} = 0.0067,

(69)

where

| \cdot |

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

Δ = 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

Δ = 5

ms. Moreover, the new estimates were directly obtained by replacing

Δ = 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].

Author Contributions

Each author played essential roles in the development of this research. Conceptualization, M.C. and J.C.A.; formal analysis, M.C., R.C. and J.C.A.; funding acquisition, M.C., R.C. and J.C.A.; investigation, M.C., R.C., P.E. and J.C.A.; methodology, M.C. and J.C.A.; resources, J.C.A.; software, M.C.; validation, P.E.; visualization, P.E.; optics analysis, P.E.; writing—original draft preparation, M.C., R.C., P.E., and J.C.A.; writing—review and editing, M.C., R.C., P.E., J.C.A. All authors have read and agreed to the published version of the manuscript.

Funding

PIIC-004/2019 Universidad Técnica Federico Santa María; Comisión Nacional de Investigación Científica y Tecnológica (CONICYT) (CONICYT-PFCHA/Doctorado Nacional/2017- 21170804); ANID-FONDECYT 1211630 and 11201187; Advanced Center for Electrical and Electronic Engineering, AC3E, Basal Project FB0008, ANID, Chile.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analysed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

References

Conan, J.M.; Raynaud, H.F.; Kulcsár, C.; Meimon, S. Are integral controllers adapted to the new era of ELT adaptive optics? In Proceedings of the AO for ELT 2011—2nd International Conference on Adaptive Optics for Extremely Large Telescopes, Victoria, BC, Canada, 25–30 September 2011; p. 57. [Google Scholar]
Hayward, T.; Ripp, M.; Bonnet, H.; Cavedoni, C.; Galvez, R.; Gausachs, G.; Cho, M. Characterizing the Vibration Environments of the Gemini Telescopes. In Proceedings of the SPIE Astronomical Telescopes + Instrumentation, Ground-based and Airborne Telescopes VI, Edinburgh, UK, 26 June–1 July 2016; pp. 1956–1963. [Google Scholar]
Tyson, R. Adaptive Optics Engineering Handbook, 1st ed.; Marcel Dekker Inc.: New York, NY, USA, 2000. [Google Scholar]
Kulcsár, C.; Raynaud, H.; Conan, J.; Correia, C.; Petit, C. Control design and turbulent phase models in adaptive optics: A state-space interpretation. In Proceedings of the Adaptive Optics: Methods, Analysis and Applications, San José, CA, USA, 13–15 October 2009. AOWB1. [Google Scholar]
Petit, C.; Conan, J.M.; Kulcsár, C.; Raynaud, H.; Fusco, T. First laboratory validation of vibration filtering with LQG control law for adaptive optics. Opt. Express 2008, 16, 87–97. [Google Scholar] [CrossRef]
Sivo, G.; Kulcsár, C.; Conan, J.M.; Raynaud, H.; Éric, G.; Basden, A.; Vidal, F.; Morris, T.; Meimon, S.; Petit, C.; et al. First on-sky SCAO validation of full LQG control with vibration mitigation on the CANARY pathfinder. Opt. Express 2014, 22, 23565–23591. [Google Scholar] [CrossRef] [PubMed]
Guesalaga, A.; Neichel, B.; Rigaut, F.; Osborn, J.; Guzman, D. Comparison of vibration mitigation controllers for adaptive optics systems. Appl. Opt. 2012, 51, 4520–4535. [Google Scholar] [CrossRef]
Garcés, J.; Zúñiga, S.; Close, L.; Males, J.; Morzinski, K.; Escárate, P.; Castro, M.; Marchioni, J.; Rojas, D. Vibrations in MagAO: Resonance Sources Identification and First Approaches for Modeling and Control. In Proceedings of the SPIE Astronomical Telescopes + Instrumentation, Adaptive Optics Systems V, Edinburgh, UK, 26 June–1 July 2016; pp. 1079–1095. [Google Scholar]
Zúñiga, S.; Garcés, J.; Close, L.; Males, J.; Morzinski, K.; Escárate, P.; Castro, M.; Marchioni, J.; Zagals, D.R. Vibrations in MagAO: Frequency-Based Analysis of on-Sky Data, Resonance Sources Identification, and Future Challenges in Vibrations Mitigation. In Proceedings of the SPIE Astronomical Telescopes + Instrumentation, Adaptive Optics Systems V, Edinburgh, UK, 26 June–1 July 2016; pp. 1123–1128. [Google Scholar]
Escárate, P.; Christou, J.; Rahmer, G.; Miller, D.; Hill, J. Understanding the vibration environment for LBT/AO. In Proceedings of the AO4ELT5 Conference, Tenerife, Spain, 25–30 June 2017. [Google Scholar]
Glück, M.; Pott, J.; Sawodny, O. Model predictive control of multi-Mirror adaptive optics systems. In Proceedings of the 2018 IEEE Conference on Control Technology and Applications (CCTA), Copenhagen, Denmark, 21–24 August 2018; pp. 909–914. [Google Scholar]
Haber, A.; Verhaegen, M. Modeling and state-space identification of deformable mirrors. Opt. Express 2020, 28, 4726–4740. [Google Scholar] [CrossRef] [Green Version]
Mocci, J.; Quintavalla, M.; Chiuso, A.; Bonora, S.; Muradore, R. PI-shaped LQG control design for adaptive optics systems. Control Eng. Pract. 2020, 102, 104528. [Google Scholar] [CrossRef]
Sedghi, B.; Müller, M.; Jakob, G. E-elt vibration modeling, simulation, and budgeting. In Proceedings of the Integrated Modeling of Complex Optomechanical Systems II, Varenna, Italy, 7–9 October 2015; pp. 1–6. [Google Scholar]
Sedghi, B.; Müller, M.; Dimmler, M. Analyzing the impact of vibrations on e-elt primary segmented mirror. In Proceedings of the SPIE Astronomical Telescopes + Instrumentation, Modeling, Systems Engineering, and Project Management for Astronomy VII, Edinburgh, UK, 26–28 June 2016. 991111. [Google Scholar]
Escárate, P.; Coronel, M.; González, K.; Carvajal, R.; Agüero, J.C. Vibration model identification using the Maximum Likelihood method. In Proceedings of the SPIE Adaptive Optics Systems VI, Austin, TX, USA, 10–15 June 2018. [Google Scholar]
González, K.; Coronel, M.; Carvajal, R.; Escárate, P.; Agüero, J.C. Maximum Likelihood identification of a continuous-time oscillator utilizing sampled data. IFAC-PapersOnLine 2018, 51, 712–717. [Google Scholar] [CrossRef]
Coronel, M.; Escárate, P.; Christou, J.; Hill, J.; Rahmer, G.; Carvajal, R.; Agüero, J.C. Vibrations modelling at LBTO utilizing telemetry data. In Proceedings of the Adaptive Optics for Extremely Large Telescopes, Quebec City, QC, Canada, 9–14 June 2019. [Google Scholar]
Soto-Muñoz, N.; Langedijk, C.; Escárate, P.; Carvajal, R.; Agüero, J.C. Identification and control of an experimental Adaptive Optics setup. In Proceedings of the 2019 IEEE CHILEAN Conference on Electrical, Electronics Engineering, Information and Communication Technologies (CHILECON), Valparaíso, Chile, 13–27 November 2019. [Google Scholar]
Glück, M.; Pott, J.; Sawodny, O. Investigations of an accelerometer-based disturbance feedforward control for vibration suppression in adaptive optics of large telescopes. Astron. Soc. Pac. 2017, 129, 976. [Google Scholar] [CrossRef]
Kubin, G.; Lainscsek, C.; Rank, E. Identification of nonlinear oscillator models for speech analysis and synthesis. In Nonlinear Speech Modeling and Applications: Advanced Lectures and Revised Selected Papers; Chollet, G., Esposito, A., Faundez-Zanuy, M., Marinaro, M., Eds.; Springer: Berlin/Heidelberg, Germany, 2005; pp. 74–113. [Google Scholar]
Kulcsár, C.; Raynaud, H.F.; Petit, C.; Conan, J.M.; de Lesegno, P.V. Optimal control, observers and integrators in adaptive optics. Opt. Express 2006, 14, 7464–7476. [Google Scholar] [CrossRef]
Tyson, R. Principles of Adaptive Optics, 4th ed.; Series in Optics and Optoelectronics; CRC Press: Boca Raton, FL, USA, 2016. [Google Scholar]
Kulcsár, C.; Raynaud, H.F.; Petit, C.; Conan, J.M. Minimum variance prediction and control for adaptive optics. Automatica 2012, 48, 1939–1954. [Google Scholar] [CrossRef]
Conan, J.M.; Rousset, G.; Madec, P.Y. Wave-front temporal spectra in high-resolution imaging through turbulence. J. Opt. Soc. Am. A 1995, 12, 1559–1570. [Google Scholar] [CrossRef]
Escárate, P.; Carvajal, R.; Close, L.; Males, J.; Morzinski, K.; Agüero, J.C. Minimum variance control for mitigation of vibrations in adaptive optics systems. Appl. Opt. 2017, 56, 5388–5397. [Google Scholar] [CrossRef] [PubMed]
Sedghi, B.; Müller, M.; Bonnet, H.; Dimmler, M.; Bauvir, B. Field stabilization (tip/tilt control) of E-ELT. In Proceedings of the SPIE Astronomical Telescopes + Instrumentation, Ground-based and Airborne Telescopes III, San Diego, CA, USA, 27 June–2 July 2010; pp. 1402–1413. [Google Scholar]
Rodriguez, I.; Neichel, B.; Guesalaga, A.; Rigaut, F.; Guzman, D. Kalman and H-infinity controllers for GeMS. In Proceedings of the Adaptive Optics: Methods, Analysis and Applications 2011, Toronto, ON, Canada, 10–14 July 2011; p. JWA32. [Google Scholar]
Juvénal, R.; Kulcsár, C.; Raynaud, H.; Conan, J.M.; Petit, C.; Leboulleux, L.; Sivo, G.; Garrel, V. Tip-tilt modelling and control for GeMS: A performance comparison of identification techniques. In Proceedings of the Adaptive Optics for Extremely Large Telescopes 4, Lake Arrowhead, CA, USA, 26–30 October 2015. [Google Scholar]
Raynaud, H.F.; Correia, C.; Kulcsár, C.; Conan, J.M. Minimum-variance control of astronomical adaptive optic systems with actuator dynamics under synchronous and asynchronous sampling. Int. J. Robust Nonlinear Control 2011, 21, 768–789. [Google Scholar] [CrossRef]
Coronel, M.; Soto, N.; Carvajal, R.; Escárate, P.; Agüero, J.C. Identification and model predictive control of an experimental adaptive optics setup utilizing Kautz basis functions. In Proceedings of the SPIE Astronomical Telescopes + Instrumentation, Adaptive Optics Systems VII, Virtual Conference, 14–18 December 2020; pp. 487–498. [Google Scholar] [CrossRef]
Kulcsár, C.; Massioni, P.; Sivo, G.; Raynaud, H. Vibration mitigation in adaptive optics control. In Proceedings of the SPIE Astronomical Telescopes + Instrumentation, Adaptive Optics Systems III, Amsterdam, The Netherlands, 1–6 July 2012; p. 84470Z. [Google Scholar]
Fedrigo, E.; Muradore, R.; Zilio, D. High performance adaptive optics system with fine tip/tilt control. Control Eng. Pract. 2009, 17, 122–135. [Google Scholar] [CrossRef]
Meimon, S.; Petit, C.; Fusco, T.; Kulcsar, C. Tip–tilt disturbance model identification for Kalman-based control scheme: Application to XAO and ELT systems. J. Opt. Soc. Am. A 2010, 27, A122–A132. [Google Scholar] [CrossRef]
Yang, K.; Yang, P.; Chen, S.; Wang, S.; Wen, L.; Dong, L.; He, X.; Lai, B.; Yu, X.; Xu, B. Vibration identification based on Levenberg–Marquardt optimization for mitigation in adaptive optics systems. Appl. Opt. 2018, 57, 2820–2826. [Google Scholar] [CrossRef]
Gavel, D.T.; Wiberg, D. Toward Strehl-optimizing adaptive optics controllers. In Proceedings of the Adaptive Optical System Technologies II, Waikoloa, HI, USA, 7 February 2003; Volume 4839, pp. 890–901. [Google Scholar]
Piatrou, P.; Roggemann, M.C. Performance study of Kalman filter controller for multiconjugate adaptive optics. Appl. Opt. 2007, 46, 1446–1455. [Google Scholar] [CrossRef]
Correia, C.; Jackson, K.; Véran, J.P.; Andersen, D.; Lardière, O.; Bradley, C. Static and predictive tomographic reconstruction for wide-field multi-object adaptive optics systems. J. Opt. Soc. Am. A 2014, 31, 101–113. [Google Scholar] [CrossRef] [PubMed]
Correia, C.; Véran, J.P.; Herriot, G.; Ellerbroek, B.; Wang, L.; Gilles, L. Increased sky coverage with optimal correction of tilt and tilt-anisoplanatism modes in laser-guide-star multiconjugate adaptive optics. J. Opt. Soc. Am. A 2013, 30, 604–615. [Google Scholar] [CrossRef] [Green Version]
Jackson, K.; Correia, C.; Lardière, O.; Andersen, D.; Bradley, C. Linear prediction of atmospheric wave-fronts for tomographic adaptive optics systems: Modelling and robustness assessment. Opt. Lett. 2015, 40, 143–146. [Google Scholar] [CrossRef] [Green Version]
Goodwin, G.C.; Agüero, J.C.; Cea-Garrido, M.E.; Salgado, M.; Yuz, J. Sampling and sampled-data models: The interface between the continuous world and digital algorithms. IEEE Control Syst. Mag. 2013, 33, 34–53. [Google Scholar]
Åström, K.; Wittenmark, B. Computer-Controlled Systems Theory and Design, 3rd ed.; Prentice Hall: Hoboken, NJ, USA, 1997. [Google Scholar]
Bernstein, D. Matrix Mathematics: Theory, Facts, and Formulas, 2nd ed.; Princeton University Press: Woodstock, UK, 2009. [Google Scholar]
Garnier, H.; Wang, L. Identification of Continuous-Time Models from Sampled Data; Springer: London, UK, 2008. [Google Scholar]
Van Loan, C. Computing integrals involving the matrix exponential. IEEE Trans. Autom. Control 1978, 23, 395–404. [Google Scholar] [CrossRef] [Green Version]
Söderström, T. Discrete-Time Stochastic Systems: Estimation and Control, 2nd ed.; Springer-Verlag New York, Inc.: Secaucus, NJ, USA, 2002. [Google Scholar]
Anderson, B.; Moore, J. Optimal Filtering; Prentice-Hall Inc.: Englewood Cliffs, NJ, USA, 1979. [Google Scholar]
Goodwin, G.C.; Graebe, S.; Salgado, M. Control System Design, 1st ed.; Prentice Hall PTR: Upper Saddle River, NJ, USA, 2000. [Google Scholar]
Hannan, E. Time Series Analysis; John Wiley & Sons, Inc.: New York, NY, USA, 1960. [Google Scholar]
Pawitan, Y. Whittle Likelihood. In Encyclopedia of Statistical Sciences, 2nd ed.; Kotz, S., Balakrishnan, N., Read, C., Vidakovic, B., Eds.; John Wiley & Sons, Inc.: Somerset, NJ, USA, 2004; Volume 15, pp. 9136–9138. [Google Scholar]
Palma, W. Long-Memory Time Series Theory and Methods; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2007. [Google Scholar]
Agüero, J.C.; Goodwin, G.C. Choosing Between Open- and Closed-Loop Experiments in Linear System Identification. IEEE Trans. Autom. Control 2007, 52, 1475–1480. [Google Scholar] [CrossRef]

Figure 1. Block diagram for an AO closed-loop system.

Figure 2. Wavefront sensor.

Figure 3. Operation of a deformable mirror.

Figure 4. Equivalent block diagram model for AO system. The auxiliary variable

ρ_{k} = - φ_{k}^{cor}

Figure 4. Equivalent block diagram model for AO system. The auxiliary variable

ρ_{k} = - φ_{k}^{cor}

Figure 5. Frequency response of the multiplicative model error

H_{E} (z^{- 1}) ϑ

for different values of

ϑ

Figure 5. Frequency response of the multiplicative model error

H_{E} (z^{- 1}) ϑ

for different values of

ϑ

Figure 6. Frequency Response of the true discrete-time disturbance transfer function and the estimated discrete-time disturbance transfer function. The solid blue line represents the frequency response of the true disturbance model,

H (z^{- 1})

, and the dashed red line represents the frequency response of the estimated model with

ϑ = 1

H (z^{- 1})

, and the dashed red line represents the frequency response of the estimated model with

ϑ = 1

Figure 7. Discrete-time PSD

(| H (e^{i ω_{j} Δ}) |^{2} λ^{2})

of the disturbance signal.

Figure 7. Discrete-time PSD

(| H (e^{i ω_{j} Δ}) |^{2} λ^{2})

of the disturbance signal.

Figure 8. Discrete-time PSD of the turbulence-plus-vibrations perturbation signal in the Clay Telescope at the Las Campanas Observatory. 2014b observing run by the University of Arizona MagAO Team on the night of 31 October 2014.

Figure 9. Monte-Carlo simulation results to discrete-time PSD of six disturbance for (a) Whittle’s likelihood method and (b) NLS fitting method. We use for simulation the values show in Table 1,

Δ = 5

ms, data length

N = 500

and

N = 1000

. 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 PSDs. The gray shaded region represents the area in which all the estimated spectra lie.

Δ = 5

ms, data length

N = 500

and

N = 1000

Figure 10. Discrete-time PSD of the disturbance (input signal) and controlled outputs. We use the estimated values that are shown in Table 3 and Table 4,

Δ = 5

ms, and

N = 1000

. The densely dotted red line represents the discrete-time PSD of disturbance or input signal, the solid orange line with square mark represents discrete-time PSD of the controlled output using the true model, the dashed blue line represents discrete-time PSD of the controlled output using the model obtained by Whittle’s likelihood method, the dotted black line represents discrete-time PSD of the controlled output using the model obtained by NLS fitting method.

Figure 10. Discrete-time PSD of the disturbance (input signal) and controlled outputs. We use the estimated values that are shown in Table 3 and Table 4,

Δ = 5

ms, and

N = 1000

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

$α_{l} (Hz)$	2	14	20	29	43	60
$ζ_{l}$	0.9	0.05	0.05	0.05	0.05	0.05
$β_{l} (\times 10^{4})$	1.439	3.493	3.734	4.253	5.540	9.711

Table 2. Performance of MVC for different disturbance model.

$ϑ$	0	0.5
$E (%)$	0	3.85

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

N = 1000

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

N = 1000

	$\hat{α}$ (Hz)	$\hat{ζ}$	$\hat{β} (\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$

Table 4. Estimated parameters of six disturbance for NLS fitting method, using

N = 1000

Table 4. Estimated parameters of six disturbance for NLS fitting method, using

N = 1000

	$α$ (Hz)	$ζ$	${\hat{γ}}_{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$

Table 5. Performance of MVC

(E)

for different identifications techniques of disturbance, and different values of

ϑ

. We use

Δ = 5

ms, and

N = 1000

Table 5. Performance of MVC

(E)

for different identifications techniques of disturbance, and different values of

ϑ

. We use

Δ = 5

ms, and

N = 1000

	$\hat{H} (z^{- 1})$
$ϑ$	Proposed Method	NLS [35]
0	0	0
$0.5$	$1.56 %$	$3.85 %$
1	$0.870 %$	$16.52 %$
$1.5$	$2.29 %$	$44.79 %$

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Coronel, M.; Carvajal, R.; Escárate, P.; Agüero, J.C. Disturbance Modelling for Minimum Variance Control in Adaptive Optics Systems Using Wavefront Sensor Sampled-Data. Sensors 2021, 21, 3054. https://doi.org/10.3390/s21093054

AMA Style

Coronel M, Carvajal R, Escárate P, Agüero JC. Disturbance Modelling for Minimum Variance Control in Adaptive Optics Systems Using Wavefront Sensor Sampled-Data. Sensors. 2021; 21(9):3054. https://doi.org/10.3390/s21093054

Chicago/Turabian Style

Coronel, María, Rodrigo Carvajal, Pedro Escárate, and Juan C. Agüero. 2021. "Disturbance Modelling for Minimum Variance Control in Adaptive Optics Systems Using Wavefront Sensor Sampled-Data" Sensors 21, no. 9: 3054. https://doi.org/10.3390/s21093054

APA Style

Coronel, M., Carvajal, R., Escárate, P., & Agüero, J. C. (2021). Disturbance Modelling for Minimum Variance Control in Adaptive Optics Systems Using Wavefront Sensor Sampled-Data. Sensors, 21(9), 3054. https://doi.org/10.3390/s21093054

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

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

Abstract

1. Introduction

2. AO Systems

2.1. Wavefront Sensor

2.2. Deformable Mirror

2.3. AO Controller

3. Disturbance Model in AO Systems

3.1. Equivalent AO System Model

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

4. Proposed Modelling for Disturbances

5. Identification of Disturbances

5.1. Nonlinear Least Square Fitting Method

5.2. Whittle’s Likelihood

6. MVC Performance in AO Systems

6.1. Minimum Variance Control Design

6.2. Performance of MVC Subject to Model Error

6.3. Control Performance under Model Mismatch

7. Numerical Example

7.1. Disturbance Identification

7.2. Performance of MVC in AO System

8. Discussion and Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI