Open Access
Issue
A&A
Volume 650, June 2021
Article Number A41
Number of page(s) 19
Section Astronomical instrumentation
DOI https://doi.org/10.1051/0004-6361/202040216
Published online 03 June 2021

© V. Deo et al. 2021

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Thanks to its extreme sensitivity as compared to other general-purpose wavefront sensors (WFSs) for astronomical adaptive optics (AO), the pyramid WFS (PWFS; Ragazzoni 1996) has been the design choice included in all first-light AO systems for the three upcoming extremely large telescopes (ELTs; Tamai et al. 2018; Fanson et al. 2018; Liu & Sanders 2018) as well as for a number of high-performance systems on 8–10 m telescopes (Esposito et al. 2010; Guyon et al. 2010; Schatz et al. 2018).

In recent years, the community has shown interest in tackling the nonlinearity of the PWFS, the so-called optical gain (OG), which is an on-sky sensitivity reduction induced by the limited dynamic range of the sensor and can be modeled as a function of wavefront spatial frequency, with a magnitude depending on residual wavefront conditions, and therefore on ongoing turbulence statistics (Costa 2005; Korkiakoski et al. 2008a; Deo et al. 2018a). When unmeasured and not compensated for, OG is a significant hindrance that prevents the application of many algorithms relying upon the linearity of the servo-loop. Two extremely common but critical examples are (1) the explicit estimation of the temporal transfer function, which is necessary for applying modal gain optimization techniques (Gendron & Léna 1994; Dessenne et al. 1998), and (2) the proper subtraction of the noncommon path aberration setpoint (examined in, e.g., Esposito et al. 2020). For setpoint subtraction, the PWFS triggers a divergent positive feedback when attempting convergence to a wavefront setpoint whose gradient exceeds the modulation radius.

Our previously proposed modal OG compensation pipeline (Deo et al. 2019a) demonstrated significant performance improvements for different seeing conditions, but had some limitations: It required introducing probe signals on the deformable mirror (DM), it required a significant amount of preliminary modeling and computations, and it did not incorporate any modal variance minimization technique depending on noise level variations. We propose in the present paper a novel algorithm, a self-regulating and entirely automated real-time modal gain controller, that resolves these three shortcomings: the correlation-locking optimization scheme (CLOSE). The method is inspired by Montera et al. (2018), proposing to apply neural networks in particular for tip-tilt sensitivity tracking and compensation. Preliminary results with CLOSE were published in Deo et al. (2019b).

CLOSE monitors the integrator overshoot through the temporal autocorrelation of modal WFS measurements and drives modal gains through real-time multiplicative updates. This allows tracking and optimizing the sole temporal properties of the loop integrator. The OG compensation multipliers are automatically factored in the gain that is set into the command law and do not require separate explicit computations. The resulting steady-state command law may then be optimized regardless of the OG sensitivity reduction and can easily be tuned so as to compare to the minimum variance (MV) control of Gendron & Léna (1994).

While we started designing CLOSE in the context of the single-conjugate AO (SCAO) module of the MICADO instrument (Davies et al. 2018; Clénet et al. 2019), the technique naturally extends to a range of AO systems, namely for nonlinear sensors performing well under a locally valid, modal, linearized description, such as the PWFS. This includes other Fourier-like wavefront sensors (Fauvarque et al. 2019; Chambouleyron et al. 2020) as well as quad-cell Shack–Hartmanns or the adaptation of a diffraction-limited, point-source AO calibration into an adequate command law for extended natural sources or elongated laser stars, etc.

This paper is organized as follows: in Sect. 2 we present our formal model for analyzing the AO loop dynamics when using a WFS prone to optical gain. In Sect. 3 we present the rationale for the CLOSE servo-loop as well as semi-analytical demonstrations showing the achieved correspondence with minimum variance integrators. Then we present in Sect. 4 two proposed implementations, real-time and offline, and briefly lay out the computational requirements associated with the use of CLOSE. We show in Sect. 5 the results of end-to-end numerical simulations, showing convergence, improved performance, and the breadth of applicable conditions for CLOSE in the MICADO SCAO system. Finally, Sect. 6 includes some discussions of possible extensions and limitations of the scheme.

2. Modeling the AO loop

We show in Fig. 1 the schematic of a SCAO loop that is considered throughout this paper. This model that was designed to conveniently describe the OG effect on PWFSs operating in low Strehl regimes. We recall here the key elements relevant to this study. A more extensive description of this approach may be found in Deo et al. (2019a).

thumbnail Fig. 1.

General modeling of the SCAO loop showing the WFS, real-time computer, and DM. The WFS is modeled using the confusion matrix framework (Sect. 2.1). The confusion matrix A is a random variable depending upon the residual wavefront ϕRes.

Open with DEXTER

While the wavefronts defined in Fig. 1, ϕAtm, ϕDM, and ϕRes, are continuous functions over the telescope aperture, but are here implicitly meant as their decomposition over a control basis of the DM: (ϕ1, …, ϕN), where N is the number of controlled modes, plus an additional fitting component beyond the DM capabilities. The left half of Fig. 1 is represented in this modal space. The right half, where the measurement noise p is introduced, is in the WFS measurement space. For a Shack–Hartmann, this would be the centroid displacement space, or for a PWFS, the space of either the normalized pixels or the gradient-like slopes maps.

The WFS is represented by the OG-describing matrix A (Sect. 2.1) and the transformation from modal decomposition to measurement dWFS: the modal interaction matrix of the wavefront sensor, defined as the Jacobian of the WFS response around the flat wavefront, that is, computed using infinitesimal push-pulls. Wavefront reconstruction is performed by matrix-vector multiplication, with the matrix Rec computed as the generalized inverse of dWFS, assuming the latter is adequately conditioned. Finally, Fig. 1 shows that the temporal control of the loop is operated through a modal integrator using a gain vector G = [Gi],1 ≤ i ≤ N.

2.1. Optical gain: The confusion matrix model

In the general case, the small-signal response of a nonlinear sensor is modified by the presence of atmospheric wavefront residuals. The WFS Jacobian dWFS(ϕRes) around any nonzero setpoint ϕRes may be written as dWFS(ϕRes) = A ⋅ dWFS, as shown in Fig. 1. In this description, the modal space operator A describes a local mixing of model components around ϕRes as compared to the calibrated response, and we therefore call A the modal confusion matrix. The operator A generally is a random variable dependent on ϕRes. In the simplest case, for instance, a Shack–Hartmann WFS with uniform centroid gain, the matrix A is a scaled identity matrix.

For the PWFS case, it has been shown that the confusion matrix has some reasonable properties when described on an appropriate modal basis. These properties are the foundation of optical gain modal compensation for the PWFS (Korkiakoski et al. 2008b). In previous work (Deo et al. 2018a, 2019a), we performed a thorough numerical assessment of the fluctuations of A when the spatial power spectrum density (PSD) of ϕRes is stationary. These analyses were performed using a Karhunen–Loève (KL) basis orthonormalized on the DM (Ferreira et al. 2018a). This basis is made of modes containing an isotropic mix of spatial frequencies of a single norm, sorted by increasing frequency. The last ∼300 modes (out of 4300 total modes for the ELT), whose structure is affected by the DM cutoff, contain a variety of waffles. We use this basis for all purposes in this paper.

When a decomposition is used on our DM KL basis, we have previously demonstrated that (1) A is essentially diagonal for low-order modes, which bear most of the power of the atmospheric turbulence; (2) that its diagonal coefficients vary by no more than a few percent for a given set of wavefronts ϕRes of identical PSD, a property in accordance with convolutional PWFS descriptions (Fauvarque et al. 2019; Chambouleyron et al. 2020); and (3) that the off-diagonal portion of A, while non-negligible for high-order ϕi, is of negligible statistical average.

Using this basis, which statistically diagonalizes A, enables the modal gain compensation strategy, as shown in the Reconstruction block of Fig. 1. The modal gains set in G cover two functionalities: First, they include the compensation of A for ongoing turbulent conditions; and second, they define the transfer function gain such that the modal integrators exhibit appropriate rejection levels. Both of these terms are ambiguously mixed into the coefficients of G set into the controller. This ambiguity affects the dynamical modeling of the system, which we describe in Sect. 2.3, as well as the optimization strategies explained in Sect. 3. Its relation to proper noncommon path aberration compensation, which requires independently identifying the compensation of A, is discussed in Sect. 6.2. The stability of A against the PSD of ϕRes ensures that adequate Gi values vary with the same temporal scale as the descriptive statistical parameters (r0, L0, , …) of the turbulence.

2.2. Diagonalizing the loop model

Following through with the approximation of statistical diagonality of A, Fig. 1 can then be simplified to the flowchart shown in Fig. 2 because the reference interaction and reconstruction matrices dWFS and Rec reduce to an identity matrix. Figure 2 applies as one of N decoupled servo-loops for each of the controlled modes, with A reduced to its ith diagonal coefficient αi. We call αi the modal sensitivity reduction coefficient because 0 < αi ≤ 1 stands for a PWFS. With αi varying only with the statistical properties of the turbulence, it can be considered as a constant for the purposes of analyzing the model in Fig. 2 for stationary or slowly evolving atmospheric conditions. With these hypotheses, Fig. 2 is no more than a most classical linear integral feedback loop, but for the unknown and unmeasured αi.

thumbnail Fig. 2.

Diagonalization of the AO loop model (Fig. 1) for mode i under the approximation that A is diagonal and is thus reduced to its ith diagonal term αi.

Open with DEXTER

Furthermore, our choice of basis is descriptively similar to atmospheric KL modes (Dai 1996) for the modes that bear the most of the turbulence power. Thus, the N modal loops as shown in Fig. 2 may be approximated as uncorrelated, with their independent optimizations resulting in a global one.

2.3. Transfer functions

Each of the N modal servo-loops is entirely described by a small number of parameters: the temporal spectra of the turbulence and noise for the ith mode, and the αi (unknown) and Gi (known) scalars. The atmospheric temporal spectrum for a single KL mode is well described in the literature (Conan et al. 1995; Gendron 1995) for Kolmogorov or Von Kármán turbulence,

(1)

where f denotes temporal frequencies, temporal Fourier transforms, and fi is the mode cutoff frequency, given by the wind speed and the mode radial order. We use the spectrum of Eq. (1) as a reference template for the analyses described in Sect. 3.

The minimization objective we consider for integral control optimization is the variance of the residual,

(2)

where k is the temporal sample index. The value of the cannot be accessed directly, however, but only those of the measurements mi[k]. To perform semianalytical computations of how to achieve this minimization, as we do in Sect. 3.3, it is necessary to introduce the various transfer functions between the quantities involved.

These can easily be obtained (e.g., Madec 1999) from Fig. 2, and we here recall their expressions in a time-sampled framework. A discussion of continuous versus discrete time approaches to AO is proposed in Appendix A. We only note here that discrete-time is valid (Kulcsár et al. 2006) because (1) temporal aliasing is negligible, which is the case from the rapid decrease of Eq. (1); and (2) frequencies involved are well below Nyquist. We define the shorthand (Eq. (A.8)), where T is the sampling period and j2 = −1,

(3)

using which, we express the transfer functions, where the subscripts identify the input (either or the modal noise ni) and output ( or the modal measurement mi) considered,

(4)

(5)

and

(6)

(7)

These transfer functions are formally similar to those of an AO loop that is not affected by optical gain (αi = 1), except for an effective noise amplification, as seen from the −1/αi in Eq. (5). It follows from Eqs. (6) and (7) that the Fourier transform of the measurements can be written as

(8)

Equation (8) emphasizes that for a given turbulence PSD, is entirely defined by two parameters: the transfer function gain αiGi, and the sensitivity-adjusted signal-to-noise ratio (S/N), noted σi,

(9)

Minimum variance control involves setting Gi such that αiGi is the minimizer of the quantity expressed in Eq. (2) for the given adjusted S/N. This either requires measuring αi explicitly, as is done for instance in other optical gain compensation methods (Korkiakoski et al. 2008b; Esposito et al. 2015, 2020), or finding a proxy to indirectly infer the hidden value αiGi, and adjust Gi adequately.

3. Correlation-locking scheme

When the modal sensitivity reductions αi are unknown, the closed-loop measurements mi still contain sufficient information to retrieve h(f; αiGi) and control Gi toward optimal rejection. We intend to perform this control without computing an extensive spectrum of the measurements and without performing an explicit estimation of the system response, but leveraging the short-term temporal autocorrelation of the measurements as a proxy for the effective modal loop gain αiGi. We present in this section the correlation-locking optimization scheme (CLOSE), which drives the modal gains Gi in real-time toward a favorable solution for rejection. This section covers the foundation and steady-state solutions for the loop transfer functions achieved using CLOSE. We then present in Sect. 4 the implementation of CLOSE as a second-layer supervisory loop, taking the modal WFS measurements as inputs and operating on the gain vector.

3.1. Rationale: Using the loop resonance

For each mode, the loop described on Fig. 2 is a classical feedback loop with an integral controller of gain αiGi and total delay1τ + T. The transfer function h(f; g), as developed in Sect. 2.3, is a well-characterized high-pass filter, with a f+2 square modulus up to the roll-off, followed by a resonance peak located beyond. Some examples of h(f; g), using τ = 2T, are shown in Fig. 3 (top). h(f; g) is a stable filter up to a maximum gain value g = gcrit, which depends only on the normalized latency . As the gain g increases toward gcrit, the resonant peak sharpens and its amplitude increases. The central peak frequency increases and converges toward the critical frequency fcrit. In other terms, the lowest frequency pole of the Laplace transform associated with h(f; g) is displaced toward the imaginary axis and eventually intercepts it at s = 2jπfcrit for g = gcrit. Values of fcrit and gcrit with latency are provided in Table 1 for reference. The given general formula is demonstrated in Appendix A.

thumbnail Fig. 3.

Typical spectral components for our AO semianalytical computations. Top: input spectra for turbulence and noise, and the measurement power transfer function |hi(f; αiGi)|2 for three αiGi values. Bottom: corresponding measurement power spectra . For this example, τ = 2T, and fS = 500 Hz, yielding gcrit ≈ 0.61 and fcrit = 50 Hz.

Open with DEXTER

Table 1.

Parameters related to transfer function divergence depending on the latency of the system.

The behavior of the rejection peak as the transfer function approaches divergence is illustrated in Fig. 3, showing the PSDs for the atmospheric mode and the noise, and some examples of transfer functions and output spectra for three values of g = αiGi. The latency is τ = 2T, and the sensitivity adjusted S/N (Eq. (9)) was taken as σi = 10. Without loss of generality, we normalized the atmospheric spectrum at .

The atmospheric spectrum, filtered by the high-pass transfer function, results in the leftmost peak of seen on Fig. 3 (bottom), at the modal cutoff frequency fi = 1 Hz. Near 10 Hz, we see a transition from rejected turbulence, showing a f−11/3 spectrum, into a regime of f+2, as the input is dominated by the measurement noise. This regime is followed by the rejection peak, which amplifies and shifts toward fcrit = 50 Hz as the gain increases.

3.2. Correlation locking: Steady-state objective

We assume that the AO latency τ has been calibrated and is a fixed known parameter; consequently, so is fcrit. For reasonable τ values, fcrit lies in the noise floor of the spectrum, well beyond the turbulence cutoff frequency. In its simplest form, our philosophy is to note that the amplitude of the resonant peak may be used as a proxy for the effective modal gain αiGi.

We deemed it particularly inconvenient, however, to attempt real-time estimations of the peak amplitude or structure, which calls for extensive buffer acquisitions, explicit PSD estimations, etc. The latter remains well within technical reach and may be done at later stages of this research.

Instead, we propose to use the autocorrelation (AC) of the modal measurements, noted , because it provides another indirect measure of αiGi. For reference, the AC curves at small time-shifts δk corresponding to the measurement spectra shown in the bottom panel of Fig. 3 are shown in Fig. 4. As αiGi increases and as the loop response approaches divergence, an oscillation of half-period nearing

(10)

thumbnail Fig. 4.

Autocorrelation functions corresponding to the measurement PSDs shown in Fig. 3. δkcrit = 5 frames for .

Open with DEXTER

is superimposed on the typically wider bell-shaped AC curve. This oscillation is the correspondence in the AC domain of the resonance peak, growing in amplitude and converging to fcrit.

In particular, the first minimum of the AC curve for a time-shift of δkcrit is a strongly marked signal that allows measuring the amount of resonance in the loop with a minimal latency of δkcrit frames. This value is the one metric that is monitored by CLOSE to implicitly register the value of the loop gain αiGi,

(11)

which, as we show in Sect. 4, is also convenient to track with simple real-time estimators. Decreasing αiGi reduces the amplitude of the oscillation over the AC function, which tends to increase the value of , and the opposite occurs when αiGi is increased (within the stability limits).

Because the relation between αiGi and is monotonic, we can act on Gi in order to lock the δkcrit correlation value onto a steady-state solution,

(12)

where we call r the setpoint value. Without additional knowledge, the value of r ought to be adjusted (or defined per mode) to fit a performance-maximizing criterion in all useful situations that the AO would face and across the complete range of the effective modal S/N. The automatic driving of Gi to satisfy Eq. (12) is of course particularly convenient if and only if a unique or a small number of setpoint values may be found.

The δkcrit correlation value and the value r we target to lock it onto are empirically representative of the spectral energy ratio between the low-frequency atmospheric rejection residual and the overshoot peak near fcrit. This is shown in Figs. 3 and 4.

Generally speaking, using higher values for r imposes a more cautious and robust steady-state control solution with a lower loop gain, with a strong, near unit correlation at δkcrit separation. Lower values lead to more aggressive loop behaviors that might reach nearly divergent transfer functions, but with a maximized rejection of the low-frequency components. As αiGi approaches gcrit, goes toward −1, and the output of the system is a slowly dampened sinusoid of period 2δkcrit frames.

Furthermore, with the condition of Eq. (12) satisfied, CLOSE enforces a transfer function constraint that is independent of the sensitivity reduction αi of the WFS. In other words, it automatically compensates for the optical gain effect using modal compensation coefficients.

3.3. Correlation-locking and minimum variance: Semianalytical solutions

We have introduced the concept of a correlation-locking condition to adjust the controller gains and achieve a given control transfer function independently of the sensitivity reduction. We propose to compare the solutions obeying Eq. (12) to a minimum variance (MV) modal integrator, verifying Eq. (2).

For the example developed in Figs. 3 and 4, the gain αiGi minimizing the variance of the wavefront residual is gMV ≈ 0.231. At this gain value, we obtain = − 0.025. This numerical observation indicates that r = 0 should be pursued as a candidate setpoint.

The computations are easily generalized for values other than the sensitivity-corrected S/N σi = 10 used for the examples. We show in Fig. 5 the variations of and of the residual variance (Eq. (2)) with αiGi, for S/N values σi = 1, 10, and 100. We also plot the two components of the variance that are due to the turbulence residual and the noise propagation, that is, the two components that the CLOSE servo-loop seeks to balance optimally. The observed variations of confirm its monotonicity with the gain as well as the relationship inferred in Sect. 3.2: a negative observable (or setpoint r) relates to a high gain loop, and a positive to a low gain loop. For all three S/Ns, Fig. 5 show the approximate match between the gain yielding minimum variance and the intercept = 0. This further establishes that r = 0 should be pursued as a special value that can be used in a variety of situations.

thumbnail Fig. 5.

For effective S/N σi = 1, 10, and 100, closed-loop residual variances (top panels) and resulting δkcrit autocorrelation of the measurements (bottom panels) as a function of the loop gain αiGi from 0 up to gcrit. The latency is taken as τ = 2.T. Black lines show the total residual variance normalized to 1 at its minimum, blue represents the variance induced by turbulence , and green shows the variance induced by noise .

Open with DEXTER

Generalizing these analyses further, we show in Fig. 6 the values of αiGi resulting in setpoints r of −0.5, −0.1, 0, 0.1, and 0.5, together with the MV solution for effective S/N values 10−1–104. The gain gMV minimizing the variance for a given S/N is found following Gendron & Léna (1994) by numerically minimizing the joint Eqs. (4) and (5). CLOSE solutions are found by numerically solving =r for αiGi through Eqs. (6) and (7). Additional data for τ = 0 and τ = T are presented in Appendix B.

thumbnail Fig. 6.

Minimum variance gain gMV and solutions found using CLOSE with five different setpoint values, depending on the sensitivity-corrected S/N σi. These computations are performed using the modal input spectrum of Eq. (1) with fi = 1 Hz.

Open with DEXTER

First, we note that the correlation-locked solution for an empirical r = 0 is a remarkably close match to the MV solution gMV. The discrepancy at worst reaches 20% and 15% of gMV in the two knees of the curve, near σi = 3 and σi = 100, respectively. We also confirm the statement made earlier: Positive setpoints lead to lower values of αiGi for a given σi, and negative ones lead to a higher gain value. Interestingly, negative setpoints impose a gain floor even when the S/N is extremely poor. To achieve r < 0 over essentially white noise, the controller forcefully increases the gain so as to introduce a detectable oscillation of the amplified noise. Overall, r = 0 seems to provide an approximate near-MV solution over most of the σi range, with some room for improvement using a combination of r > 0 values for the higher end of the S/N range. The exact discrepancy between the gMV and the CLOSE solution at null setpoint of course depends on the actual spectrum of turbulence, in particular, on the cutoff frequency fi.

These semianalytical simulations altogether show that reaching the CLOSE steady-state solution using r = 0 can provide a near minimum variance solution throughout the entire range of σi. This makes it a relevant control technique for all the controlled DM modes for any choice of guide-star magnitude, for any amount of input turbulence, and for any sensitivity reduction αi because all these parameters merely factor in the computation of σi. With a unique setpoint value, semianalytical computations offer a simple criterion that would enable fully automatic (almost) minimum variance integral control even for sensors with poorly modeled sensitivity variations.

We validate these claims using end-to-end numerical simulations in Sect. 5. We propose in the next section a practical implementation of the steady-state equations presented above.

4. Practical implementation

4.1. Real-time

Modal WFS telemetry mi[k] is obtained at each time-step k as the product of the WFS output by the modal reconstructor Rec. Two online autocovariance estimators are built from the mi[k] using discrete integrators,

(13)

where k is the time-step index and p ∈ [0, 1] is a smoothing parameter. and are thus low-pass time-filtered estimates such that tracks the series variance Vark(mi[k]), and their ratio tracks the critical autocorrelation . After an empirical optimization of the parameter p, we opted for fast integrators with p = 0.3 for all simulations presented in Sect. 5. One objective achieved through this smoothing parameter is that the gap of δkcrit frames between the joint estimations of the variance and of is bridged.

After the AC estimation, the modal gains Gi are updated using multiplicative increments as follows:

(14)

The r parameter is the loop setpoint as defined in Sect. 3.

The q± learning factor may encompass two different values, with either q+ and q used depending on the sign of

q+ being used for Gi increases in Eq. (14), and q for decreases. This asymmetry is kept as an option to make the algorithm more reactive to overshooting transients (using q > q+) as compared to tracking gain increases due to a transfer function that is deemed too slow (using q+).

This asymmetric tracking is not used for results shown in this paper, but was used in preliminary numerical simulations (and throughout Deo et al. 2019b). It allows maintaining stability at higher q± values. However, using q+ ≠ q induces a bias in the mean value achieved for that then differs from r. We empirically observed that for q = 5q+, for example, r ∈ [ − 0.2, −0.1] should be used to retrieve the performance seen with r = 0 in the q+ = q case.

The q± learning factors determine the time constants associated with the convergence and tracking ability of the CLOSE servo-loop. We infer that for a real AO system, q± values in the range of 10−3 − 10−4 should be used (assuming 500 Hz frequency), hence providing typical time constants in the 2−20 s range. The ideal choice of q± will probably remain dependent on the system and will certainly require some adjustments accounting for robustness and responsiveness to variations of turbulence conditions, vibrations, or other transient events. While the theoretical derivations were most accurate using r = 0, we do not exclude that for each system, some sort of empirical tweaking of r may be necessary to either privilege consistent stability or aggressive rejection.

4.2. Computational strain

Implementing CLOSE in a real-time fashion is of course expected to increase the AO loop computational requirements. While the AC estimations and gain updates themselves (Eqs. (13) and (14)) are negligible compared to the matrix-vector multiplication (MVM), having the mi[k] available in real-time requires to perform the reconstruction in two successive MVM steps. The first MVM converts WFS measurements into modal values, with a computational burden nearly identical to the usual overall MVM from measurements to DM commands. The second step computes DM increments from modal values, with a nearly square matrix with a size of the number of actuators. While this two-step operation is not universal, it is worth noting that it is used routinely on some instruments (e.g., Guyon et al. 2018). While this two-step technique may become a computational bottleneck in particular for high-contrast systems on ELTs, it is currently baselined for the control of the deformable M4 on ESO’s ELT (Bonnet et al. 2019).

For a typical PWFS AO system with some spatial oversampling (the pixel projected size is smaller than the DM pitch), the number of pixels that is read out is typically 5–6 times the number of actuators because there are four pupil-like images, times the square of the oversampling factor. The number of outputs of the PWFS is therefore 2.5–6 times the number of modes, depending on whether slopes-map preprocessing is used or not. The algorithmic interest and computing cost of skipping the preprocessing is discussed in Deo et al. (2018b), for example. While the first MVM execution time depends on the WFS output dimensionality, the second MVM only depends on the number of controlled modes, with a smaller but not negligible computational footprint. As an example, we measured the real-time computer (RTC) computation time for the AO simulation setup we used in Sect. 5 (see Table 2): 881.6 ± 4.2 μs using a single MVM, against 1035.6 ± 3.6 μs using two cascaded MVMs. These timings where achieved using a single Nvidia Tesla P100 graphical processor.

Table 2.

AO numerical simulation parameters.

4.3. Offline implementation

If the RTC software cannot be altered on an existing system or if the additional strain is not acceptable within the RTC specifications, CLOSE can be implemented in a block-wise manner. All estimators, gain updates, and command matrix updates are performed in offline time, certainly in another process, and preferably on another machine, over batches of contiguously recorded measurements. This buffered strategy enables deploying CLOSE on nearly any existing AO system that reliably provides its WFS telemetry without excessive delays.

A time-continuous buffer of K WFS measurements is forwarded to the CLOSE process, which turns them into modal measurements mi, …,  mi[K − 1] using the modal reconstructor Rec. For each mode, the AC estimators of Eq. (13) are replaced by the direct computation of the normalized δkcrit-shifted AC term over the telemetry buffer,

(15)

The gain-updating equation can then be performed,

(16)

using q± factors adjusted for the longer integration time and the increased SNR on AC estimation. Typically, q± ought to be larger by a factor for a dynamical effect comparable to the real-time implementation. The new command matrix can then be computed accounting for the new Gi values, and when all side-tracked computations are completed, can be set into the RTC.

5. Numerical simulation results

This section describes some end-to-end numerical simulations that demonstrate the performance achieved with CLOSE when applied to the MICADO SCAO design (Clénet et al. 2019; Vidal et al. 2019). The main parameters of the system and the simulations are summarized in Table 2. All simulations were performed using the COMPASS platform (Ferreira et al. 2018b).

In Sect. 5.1 we analyze the convergence of modal gains when the AO loop is bootstrapped. In Sect. 5.2 we verify the steady-state performance achieved using r = 0 for various S/N levels. This is expanded in Sect. 5.3 by exploring the effect of different setpoint values. Finally, we show in Sect. 5.4 the dynamic behavior of CLOSE when it is exposed to sudden changes in conditions. Throughout the following sections, seeing conditions are referred to using the Fried parameter r0, always given at 500 nm; and guide-star brightnesses are identified by the R-band magnitude (MR), related to the photon flux per the zeropoint and system throughput given in Table 2. AO performance is most often measured in terms of H-band long-exposure (LE) Strehl ratios (SR), computed from simulated monochromatic point-spread functions at 1650 nm. SR comparisons given in percent are always in percentage point units, never relative variations.

5.1. Gain convergence upon closing the loop

We first investigate the dynamics of the modal gains upon closing the AO loop with CLOSE enabled. These simulations were all initialized identically, regardless of r0 or MR. We opted for the starting value Gi[k​ = ​0] = 0.5 for all modes. With the system latency τ = 2T simulated, the critical gain value was gcrit ≈ 0.61. Because the sensitivity reduction αi is always smaller than 1, this ensured that the loop was closed with a comfortable stability margin. From these initial 0.5 values, the Gi were driven by CLOSE to their steady-state values, accounting both for nonlinearity compensation and temporal variance minimization.

We show in Fig. 7 the temporal averages of the 4301 modal gains Gi for several time windows within the first two seconds after the AO loop is closed. These simulations were performed for four different cases, with r0 of 14.5 and 9.0 cm, and guide stars with brightnesses of MR = 0 and MR = 16.

thumbnail Fig. 7.

Convergence of CLOSE gains on the 2.0 s following closing of the AO loop for guide stars of MR = 0 and 16 and atmospheric r0 of 14.5 and 9.0 cm. All Gi are initialized to 0.5 at the start (blue line). Gi values are shown as averaged over the time windows given in the legend. Curves are smoothed along the i index for clarity. Final SRs are given in H band and are computed from the cumulative exposure over the last 200 ms (k ∈ [900, 1000]).

Open with DEXTER

For the bright cases, steady state is reached by frame k ≈ 500, that is, within one second. The process is slightly slower for the MR = 16 cases, with a continued convergence of the Gi between frames k = 500 and k = 1000. Simulations at MR = 0 are essentially performed with an infinite S/N. The dynamical gains αiGi therefore evolve nearby the maximum stability value gcrit, and the Gi coefficients reached in steady state essentially reflect the inverse of the PWFS sensitivity reduction. The Gi curves reached at the end of convergence in MR = 0 cases are in good accordance with the abacuses presented in previous work (Deo et al. 2018a, 2019a) that were obtained by directly measuring αi sensitivity reductions on static turbulence screens, with αi typically depending on mode index i as decreasing up to mode 30, which contains spatial frequencies corresponding to the modulation radius, then again increasing roughly as a power law up to the highest-order modes.

Some more insight into the bootstrapping of CLOSE can be gained by inspecting temporal series of modal gains, as shown in Fig. 8. These series correspond to four modes and are the same data set as in Fig. 7 for r0 = 14.5 cm, MR = 0. While the AO loop bootstraps in only a few frames, the convergence of CLOSE takes longer and induces some modal gain oscillations after the initial ramp-up. The amplitude and time constants of these oscillations certainly depend on a number of parameters, and importantly on the learning factor q±. In Fig. 8, various periods are observed, from ≈1 s down to smaller oscillations typically every ∼30 frames (60 ms).

thumbnail Fig. 8.

Time series of CLOSE gains for four select modes in the 2.0 s following closing of the AO loop for the simulation r0 = 14.5 cm, MR = 0 shown in Fig. 7.

Open with DEXTER

The randomness of the turbulence screens is certainly a factor in continuous variations of Gi, and the oscillations continue even as the AO integrator and CLOSE reach steady states, beyond 1–2 s. of runtime, but with little effect on the correction quality. Figures 11 and 12 show complementary data to those in Figs. 7 and 8, including longer time series and higher Strehl ratios.

Furthermore, the physical nature of the PWFS is certainly a contributor to the Gi fluctuations as compared to a purely linear sensor. When Gi reaches values that are too high, or if we had started from Gi values higher than the steady state that is ultimately reached, the optical gain phenomenon itself helps to maintain stable suboptimal control states while the controller performs the slow gain decrease. To summarize this effect, reduced atmospheric residuals (e.g., bootstrapping or improving conditions) induce an increase in WFS sensitivity. When the gains in Gi are overset, transfer functions become highly resonant or temporarily unstable. In turn, the added wavefront residual from transitory divergence of the loop or noise oscillations induces a reduction in the PWFS sensitivity, however. This regime is progressively stabilized as the adaptive filtering eventually adapts the command law to the ongoing conditions and reaches a steady-state regime. These transitory effects are further discussed in Sect. 5.4.

When we compare in Fig. 7 the behavior in the MR = 0 and MR = 16 cases, we observe the effect of the implicit optimization of the transfer function, with steady-state gain values dampened by typically 20–50% in the dim case relative to the bright one, depending on the mode number and r0. The results presented in Fig. 7 tend to validate that without any priors and regardless of the PWFS sensitivity reduction, CLOSE successfully drives the modal integrator gains to convergence in a period of 1–2 s.

5.2. Performance in stationary conditions for r = 0

In addition to the adaptive capability of CLOSE shown in Sect. 5.1, we investigated the steady-state AO performance. In order to perform this analysis, we generalized similar simulations as performed in Sect. 5.1 to a wider range of r0 values (based on statistics provided by ESO within the frame of the development of ELT instruments) and guide-star magnitudes (11.5–17.5).

The measured performances are shown in Fig. 9 (top), with the long-exposure H-band SR plotted against the star magnitude and computed for five different seeing conditions. For all results the SR was averaged over 2 s exposures, starting 3 s after the AO loop was closed, and with initial modal gains Gi = 0.5 as previously. These simulations were also performed with a manually optimized scalar integrator (Gi= constant, with one constant value per seeing and magnitude), and the performance is shown in Fig. 9 (bottom). This again shows the improvement that is achieved using fine-tuned modal control, which was demonstrated in previous research. In particular, using modal control enables us to (1) improve the performance in poor seeing, even with bright guide stars, because nonlinearity is then a dominant member of the error budget. It also enables us (2) to increase the performance at the faint end, with a gain of typically +0.5 mag for an identical objective.

thumbnail Fig. 9.

Long-exposure SR in H band obtained for end-to-end simulations of the MICADO SCAO setup with CLOSE (top), and with a manually optimized scalar integrator (bottom) for guide-star magnitudes MR = 11.5–17.5 and seeing conditions r0 = 8.9–21.5 cm.

Open with DEXTER

While the simulation setup is not a perfect simulation of the MICADO system, Fig. 9 (top) is indicative of the performance that could be achieved using CLOSE on such a system. If we exclude the worst seeing conditions r0 = 8.9 cm, a performance of 70% SR or better is always achieved for bright stars, it even reaches 89% in excellent seeings. 40% or better is achieved for MR up to 16. SRs better than 40% are also achieved for bright stars for r0 = 8.9 cm. This potentially enables us to salvage poor-seeing nights for a number of scientific cases with milder correction requirements. For the metric used here, CLOSE performs equally to other previously introduced modal compensation techniques for the PWFS (Korkiakoski et al. 2008b; Deo et al. 2019a; Chambouleyron et al. 2020), but with the added value of automation and adaptability, without the offline computations, seeing estimations, dithering signals, or modifications to the optical setup required by these methods.

5.3. Validation of the setpoint

We have shown throughout Sects. 3 and 4 that CLOSE is entirely configured with very few parameters, namely, δkcrit, which ought to be chosen depending on the system latency, and p, q±, and r being the choice of the operator. While the temporal filtering introduced with p and q± is easily apprehended and only affects the performance during transitory regimes, the determination of r is subject to more caution because it determines the final performance. Although semianalytical computations indicated r = 0 as an apparently universal choice, these computations were the conclusion of a number of modeling hypotheses and approximations, as described in Sect. 2.

Our objective here is to validate whether r = 0 remains an appropriate choice in most situations. We performed end-to-end simulations extending those presented in Sect. 5.2, now considering various r setpoints from −0.5 to 0.5 and wind speeds from 10 to 40 m s−1. The performances measured are shown in Fig. 10. Some minor discrepancies can be found between Figs. 9 and 10, which are explained by the different durations of simulated long exposures and the mismatched random turbulence screens. The SRs for Fig. 10 were computed on 600 ms exposures, following a 600 ms bootstrap for CLOSE and the AO. As compared to Fig. 9 (2 + 2 s), this was a necessary speedup given the large number of numerical simulations that were performed. We also note that LE SRs are generally determined to no better than 2–3% of standard deviation over the distribution of turbulence screens.

thumbnail Fig. 10.

Long-exposure SR in H band (color and text) for stationary simulations exploring different conditions regarding the seeing (outer rows), the guide-star magnitude (outer columns), the wind speed (inner rows), and the CLOSE setpoint (inner columns) for the MICADO SCAO simulation. Color scales are local to each of the subplots and do not match each other.

Open with DEXTER

In addition to the expected variations in SR with seeing, guide-star magnitude, and wind speed, we observe in Fig. 10 that the setpoint yielding the maximum SR, rmax thereafter, is almost always −0.1 or 0, except for two cases (MR = 12, r0 = 12.8 and 8.9 cm and 30 m s−1) where rmax = −0.2. When for a given r0, MR, and wind speed we obtain rmax ≠ 0, we observe that the difference with the corresponding performance at r = 0 is generally only 1–2%. The slight bias toward rmax < 0 may partly be explained by the reduced simulated time for CLOSE convergence. Cases with lower r would therefore have increased the modal gains more effectively in the allocated 300 frames because of the proportional effect introduced in Eq. (14).

The variations in SR with r are clearly determined from Fig. 10 as an asymmetric bell curve with a longer decrease on the r > 0 side. With r larger than the optimal value, the AO uses slower modal integrators with less turbulence rejection, which are more forgiving with regard to performance than r being too small. The latter cases introduce a buildup of noise that ultimately leads to diverging oscillations as r becomes too negative, hence a faster decrease in SR for simulations with r smaller than the optimum.

A few outlying cases show a significant performance gap between the maximum SR and the SR achieved for r = 0, up to a greatest difference of 6.3% (MR = 15, r0 = 12.8 cm, 20 m s−1). These cases, where r = 0 induces a noticeably suboptimal performance, are all found for wind speeds of 20 m s−1 or more. With the high latency of τ = 2 frames simulated here, the higher wind speed induces a narrow acceptable range of gain αiGi for each mode to achieve near-optimal rejection, and thus a narrow range of r values leading to this optimization with CLOSE. As discussed in Sect. 3.3, while the CLOSE solution achieved for the modal gain is an empirical close match to the optimum rejection solution, we may be reaching the limits of this approximation for cases with high latency and high wind speeds.

To investigate the usability of CLOSE with r ≈ 0 even further, we also performed similar simulations for different values of the system latency using τ = 0, 1 or 2 frames and using AO setups other than the MICADO SCAO, namely two SCAOs on an 8 m telescope, using a PWFS and a Shack–Hartmann (SH) WFS, respectively. The configuration of these simulations and the obtained results are shown in Appendix C. Beyond the results shown in Fig. 10, the extensive simulations exposed in Appendix C confirm that the setpoint r = 0 may be used as a baseline for a variety of AO systems in a large number of observation situations. The performance and steady-state behavior may eventually be fine-tuned by adjusting the r setpoint upon empirical criteria if deemed necessary.

5.4. Transients in observation conditions

Beyond the performances of correlation-locking observed in stationary conditions exposed in Sects. 5.2 and 5.3, we propose here to analyze the dynamical capabilities of CLOSE in situations where the PWFS sensitivity or the illumination vary rapidly. For this purpose, we designed two simulations with evolving conditions. The first simulation emulates a seeing burst in which atmospheric conditions degrade dramatically during a short period, and the second simulation emulates a passing cirrus, with an equivalent drop of 3 mag (−94% of flux) of the guide star.

The results for these simulations are shown in Figs. 11 and 12. We proceeded identically for both transients, simulating 8 s (4000 frames) of total runtime decomposed as 2 s of bootstrapping the AO loop from flat gains Gi(t = 0) = 0.5, followed by 4 s in degraded conditions, and finally 2 s after reverting to the initial conditions. Both Figs. 11 and 12 show the time series of the H-band SR and of the modal gains Gi for four select modes. The lower panels in Figs. 11 and 12 show a snapshot of the modal gain vector G over all 4301 DM modes at times t = 2, 6, and 8 s, that is, immediately before the changes in seeing or brightness, and at the end of the simulation.

thumbnail Fig. 11.

Top: time series of the H-band SR and modal gains Gi for a few select modes during the simulation of a burst of seeing; r0 = 21.5 cm for t < 2 s and t > 6 s and r0 = 8.9 cm for 2 < t < 6 s. Curves are smoothed using a 10 ms window. Bottom: snapshot of the modal gain vector G at time t = 2, 6, and 8 s.

Open with DEXTER

Seeing burst. For the seeing burst (Fig. 11), the turbulence screen was scaled up between t = 2 and 6 s to simulate a r0 change from 21.5 cm to 8.9 cm while keeping a guide star of brightness MR = 14. The sudden changes in r0 at t = 2 and 6 s induce discontinuities in the wavefront because we performed an instantaneous scaling across the entire aperture. This induced very short drops to 0% SR (a few frames). At t = 2 s, the modal gains optimized for r0 = 21.5 cm do not permit reaching optimal performance immediately now that r0 = 8.9 cm. For 2 < t < 3 s, CLOSE drives the Gi up by a factor 2–4 to accommodate the OG-induced sensitivity reduction, which is significantly stronger at r0 = 8.9 cm. This second long convergence to the new optimal state allows an SR improvement up to the 30–35% consistently measured for 3 < t < 6 s, which is the expected performance in these stationary conditions (r0 = 8.9 cm, MR = 14; as seen in Fig. 9).

Similarly, the second transition back to r0 = 21.5 cm at t = 6 s shows an immediate SR performance about 45% below the expectation, indicating that the system is out of tune because the Gi are now beyond reasonable values for r0 = 21.5 cm. The modal gains are automatically decreased, however, and the SR converges to 85% within 0.5 s. This is consistent with the magnitude of the time constants involved because we used a learning factor q± = 10−3 at a frame rate of 500 Hz.

Photometric variations. The “cirrus simulation” (Fig. 12) was run using the median seeing r0 = 14.5 cm, with a change of guide-star brightness from the initial MR = 14 down to MR = 17 during the perturbation interval 2 < t < 6 s. While such a flux attenuation qualitatively appears as a major change, the effect on the gain is rather mild when compared to the previous case. The reaction of CLOSE that can be observed in Fig. 12 is that an appropriate overall reduction in modal gains is progressively introduced during the transient and then reverted after t = 6 s. This behavior is typically expected from modal gain optimization with S/N. The SR temporal evolution shows that the progressive optimization of the transfer functions for MR = 17 allows for an improvement from 10 to 20% of SR within the 4 seconds spent at this brightness level.

thumbnail Fig. 12.

Similar to Fig. 11. Time series and modal gain snapshots during a cloud-like transient, simulated as a guide-star magnitude change: MR = 14 for t < 2 s and t > 6 s and MR = 17 for 2 < t < 6 s.

Open with DEXTER

It is worth noting that while the modal gain curve reached at t = 6 s only shows a modest reduction compared to t = 2 s, this reduction is in fact the product of two counterbalancing effects: First, a reduction of αiGi to improve the noise rejection, but coupled to an increase in the sensitivity compensation because the mean residual wavefront is notably increased, which induces a stronger OG sensitivity reduction in the PWFS response. These two effects are entirely implicit: As the flux drops by 3 mag, the corrective capability of the AO immediately decreases, the residual wavefront increases, and thus the sensitivity αi decreases as well. Factored in the transfer function, the gain αiGi is decreased immediately and compensates for the sudden shift in S/N. Afterward, CLOSE accounts for the fine-tuning of the Gi to more precisely set the rejection at the new sensitivity level. The opposite occurs at t = 6 s: The sharpening of the correction improves the sensitivity, and the dynamical gain αiGi immediately increases and provides more turbulence rejection, now with less noise to amplify. This effects explains why compared to the seeing burst, (1) the fluctuations in Gi are significantly milder and (2) the return-to-steady-state time is much shorter.

Finally, as with the seeing burst, we note that the performance in the final state t > 6 s is identical to the beginning t < 2 s. This indicates that no diverging signals were accumulated during the transitions and that the AO is capable of reverting cleanly after such transients.

Synthesis. From the two experiments performed with sudden changes in observing conditions, we may conclude that CLOSE efficiently and automatically accommodates sudden significant changes in S/N, whether from the seeing or the guide-star brightness. This robustness to transients is observed for improving or deteriorating conditions. This efficient convergence toward nominal performance, when the seeing improves in particular, underlines the robustness of the proposed AO control scheme. When the seeing suddenly improves, the AO is configured with gains that are significantly too high for the newly set turbulent conditions, but the control scheme allows reoptimizing the gains without the loop diverging.

6. Discussion

6.1. In general: An interaction matrix transform

In the most general sense, as developed in Sects. 2 and 3, CLOSE is an automatic minimum variance controller for servo-loops whose input is built of a low-frequency rumble and a white noise floor (Eq. (1), Fig. 3). When these conditions are met, as is the case with astronomical adaptive optics, and when independent modal loops can be decoupled, CLOSE is an excellent tool for optimizing the system regardless of miscalibrated sensitivity factors.

CLOSE is therefore adequate for any system that admits a description as in Fig. 1, where an eigenbasis can be found that properly describes the nonlinearity of the system, that is, diagonalizes (well enough) the modal confusion matrix A. This is the case for the PWFS, as has been proven experimentally (Deo et al. 2018a) and analytically (Fauvarque et al. 2019; Chambouleyron et al. 2020). This was extended to the entire class of Fourier-based wavefront sensors.

Furthermore, it has been demonstrated that a number of configuration changes of PWFSs are very well described by variations in the modal matrix A, while the interaction matrix dWFS can be kept to the unresolved, bright, on-axis, single-star matrix acquired in the laboratory (or simulated). These changes include online changes to the modulation radius, or tracking double stars of various separations and orientations relative to the pyramidal prism edges, as well as other extended objects (Titan, uncompensated atmospheric dispersion beacon; Vidal et al. 2019). For binary stars, it has in particular been shown that using CLOSE restores the performance of the servo-loop to near single-star levels, without the additional burden of biasing the loop setpoint to drive one binary component to the PWFS apex, which is the usual and nonetheless efficient alternative.

All of these transformations, similarly to optical gain, can be approximated by Fourier filters that alter the description of the sensor and thus are adequately tracked by CLOSE. Ideally, a Fourier basis would have to be used to properly diagonalize A. Our experience and numerical analyses demonstrate that a frequency-ordered KL basis is a sufficient approximation. We previously studied (Deo et al. 2018a) for the case of optical gain alone that the scaling applied by A on Fourier modes of identical norm is isotropic within a few percent. This proves that KL modes are adequate for such an operation.

6.2. Further required steps

Finally, we would like to point out a few elements that are currently beyond the scope of what CLOSE is capable of achieving. One of the most prominent issues with high-order PWFSs is the proper subtraction of non-common path aberrations (NCPAs) as the servo-loop setpoint, which issue is tightly coupled with the estimation of modal OG compensation coefficients. In the state presented in this paper, CLOSE does not enable an independent estimation of the transfer function gain αiGi and the OG compensation coefficient , as we observed in the cloud transient analyzed in Sect. 5.4. Obtaining is necessary for a proper OG-adjusted compensation of the setpoint. Other techniques, which include injecting probe signals into the system DM, allow an explicit retrieval of (Esposito et al. 2020). While CLOSE does not include this feature, it is very compatible with such techniques but is much more reactive to short-term transients. A joint implementation of several modal gain estimation techniques would be extremely beneficial for an all-in-one solution for PWFSs control, including minimum variance laws, proper NCPA compensation, shift and misalignment tracking (Héritier et al. 2018), etc.

One additional comment should be made about the potential limitations of the modeling described in Sect. 2. The proper variance minimization obtained with CLOSE relies on the clear separation of two frequency components from the turbulence and the noise in the resulting measurement PSD: . The requirements of this paper exceed this by requiring that the modal noise is white and that the turbulence spectrum obeys Eq. (1). We proved with the Shack–Hartmann results shown in Appendix C, however, that the decay exponent −17/3 of Eq. (1) may allow some tolerance because it is then actually −11/3 due to additional subaperture aliasing. Together with using different von Kármán outer scales L0, this is very likely to affect the quantitative conclusions of the semianalytical derivations described in Sect. 3.3, but we are confident that it will not compromise the monotonicity of the relationship between the δkcrit correlation and the loop gain αiGi. Thus, a simple reassessment of the loop setpoint r to try and achieve quasi-optimality is likely to be sufficient for situations like this.

Telescope vibrations are a less trivial modification of the input turbulent spectrum. Naturally, CLOSE cannot achieve any better vibration rejection than other nonpredictive integrators. As for other spectrum modifications, a reassessment of r could help in achieving near optimality, as good as this may be in the presence of narrow vibration peaks. Low-frequency (e.g., < 10 Hz) vibrations would implicitly be registered as a degradation in seeing, and would require biasing the setpoint toward r > 0 to prevent the corresponding gain increase. Conversely, vibrations between fcrit/2 and 3fcrit/2 would be registered as noise energy and require biasing the setpoint toward r < 0.

We foresee that vibrations at or near fcrit will cause specific issues. All of our simulations have shown that CLOSE, seeking its near-MV solutions, tends to push the loop gain very close to its critical value when the S/N is good enough. This could lead to transients with very high vibration amplification, for instance, as the vibration restarts after a pause during which αiGi was boosted near gcrit. These situations will require the implementation of specific fail-safes to avoid the sudden divergence of the AO loop. More generally, CLOSE in the presented version does not currently cope well with spurious divergences of the AO loop: The loss of WFS feedback associated with DM divergence (or clipping at maximum value) will induce CLOSE (if r ≤ 0) to continue increasing gains toward infinity, even though the loop has already diverged. It is now necessary to experiment with various modifications of the input spectra (correlated noise, vibrations, non von Kármán spectra) and the integral command law (leaky integrators, linear predictive methods) to investigate the extent to which and the versatility with which we may obtain satisfactory results with the version of CLOSE presented in this paper, and which further improvements are necessary.

7. Conclusion

We have presented a self-regulating adaptive filtering method that enables tracking the modal gains of the integrator controller in an AO system in hard or soft real-time. This technique allows keeping the system always close to peak performance against variable conditions (seeing, wind speed, and S/N) and the variable sensitivity of the WFS. Because it is entirely automatic and independent of WFS sensitivity fluctuations, the proposed scheme is particularly well suited for high-order PWFS designs on ELTs, which are prone to strong optical gain effects. This OG effect is well decoupled and stable against seeing statistics using a KL control basis of the DM as we use, however, thus enabling a separable modal analysis throughout.

Our method is called CLOSE, for correlation-locking optimization scheme. It leverages the monotonic relationship between the hidden actual transfer function gain and the loop resonance peak, whose amplitude is indirectly assessed through the temporal autocorrelation of the modal phase residuals computed at a single given time-shift δkcrit. The appropriate counter-reaction is then applied on the modal gain in order to lock the previously mentioned δkcrit correlation value to a chosen setpoint. Using semianalytical computations, we have demonstrated that for given conditions, there is an unequivocal relation between the transfer function gain and the value of the δkcrit correlation, and that correlation locking easily permits almost achieving a modal minimum variance criterion regardless of the modal S/N.

We have used extensive end-to-end AO simulations to show the versatility of CLOSE. Using the unique setpoint r = 0, we showed automatic convergence of the modal gains in a wide variety of simulated conditions using different wind speeds, guide-star magnitudes, and seeing parameters. We also demonstrated the automatic adaptability of the scheme to abrupt changes in observation conditions, hopping between different sets of optimal control gains in a matter of seconds. These simulations use the design parameters of the PWFS of the MICADO SCAO system, and compared to recent design studies for this system (Vidal et al. 2019), achieve nominal performance accounting both for optical gain modal compensation and modal transfer function optimization. This performance is achieved while entirely alleviating from extensive situation-dependent optimizations, database look-ups, or manual tune-ups, thus offering CLOSE as an all-in-one baseline strategy for AO modal control.

While CLOSE elegantly paves the way for solving variance minimization issues with nonlinear sensors, we cautioned in Sect. 6.2 that the scheme in the version described in this paper does not provide an explicit estimation of the sensitivity reduction coefficient nor a fully fledged solution to the NCPA subtraction issue for sensors of varying sensitivity.

Work is currently underway to refine and expand the work described in this paper to make it usable for on-sky operations. Importantly, we will design and deploy the necessary fail-safes in case CLOSE suffers from a loss of reliable feedback from the modal outputs, or if the AO falls into nonlinear measurement pits. In these cases, the modal gain update equations tend to increase the gains exponentially. These situations must be automatically identified and unambiguously distinguished from transient sensitivity reductions. If its capabilities, convenience, and versatility are confirmed, we envision CLOSE as a core control technique for future AO systems, embedded within a larger framework of leveraging real-time telemetry for AO control.


1

Expected value of the delay between the occurrence of a perturbation ϕAtm(t) and the mean time of its correction on the DM.

Acknowledgments

This research is performed in the frame of the development of MICADO, first light instrument of the ELT (ESO), with the support of ESO, INSU/CNRS and Observatoire de Paris. V. Deo is supported by NASA grant 80NSSC19K0336. Additional support was provided through the WOLF project, ANR-18-CE31-0018 of the French National Research Agency (ANR), and the OPTICON Program of the European Commission (H2020, Grant number 730890). The authors thank the COMPASS development and support team for their continued support.

References

  1. Bonnet, H., Spyromilio, J., Vérinaud, C., et al. 2019, in ELT Wavefront Control Strategy, 6th AO4ELT conference [Google Scholar]
  2. Chambouleyron, V., Fauvarque, O., Janin-Potiron, P., et al. 2020, A&A, 644, A6 [EDP Sciences] [Google Scholar]
  3. Clénet, Y., Buey, T., Gendron, É., et al. 2019, in MICADO-MAORY SCAO Preliminary Design, Development Plan and Calibration Strategies, 6th AO4ELT Conference [Google Scholar]
  4. Clergeon, C. 2014, PhD Thesis, Observatoire de Paris - Subaru NAOJ [Google Scholar]
  5. Conan, J.-M., Rousset, G., & Madec, P.-Y. 1995, J. Opt. Soc. Am. A, 12, 1559 [Google Scholar]
  6. Costa, J. B. 2005, Appl. Opt., 44, 60 [Google Scholar]
  7. Dai, G.-M. 1996, J. Opt. Soc. Am. A, 13, 1218 [Google Scholar]
  8. Davies, R., Alves, J., Clénet, Y., et al. 2018, Proc. SPIE, 10702, 107021S [Google Scholar]
  9. Deo, V., Gendron, E., Rousset, G., et al. 2018a, Proc. SPIE, 10703, 1070320 [Google Scholar]
  10. Deo, V., Gendron, E., Rousset, G., et al. 2018b, A&A, 619, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Deo, V., Gendron, E., Rousset, G., et al. 2019a, A&A, 629, A107 [CrossRef] [EDP Sciences] [Google Scholar]
  12. Deo, V., Rozel, M., Bertrou-Cantou, A., et al. 2019b, in CLOSE: a Self-regulating, Best-performance Tracker for Modal Integrator Based AO Loops, 6th AO4ELT Conference [Google Scholar]
  13. Dessenne, C., Madec, P. Y., & Rousset, G. 1998, Appl. Opt., 37, 4623 [Google Scholar]
  14. Esposito, S., Riccardi, A., Fini, L., et al. 2010, Proc. SPIE, 7736, 773609 [Google Scholar]
  15. Esposito, S., Pinna, E., Puglisi, A., et al. 2015, Non common path aberration correction with non linear WFSs, 4th AO4ELT Conference [Google Scholar]
  16. Esposito, S., Puglisi, A., Pinna, E., et al. 2020, A&A, 636, A88 [EDP Sciences] [Google Scholar]
  17. Fanson, J., McCarthy, P., Bernstein, R., et al. 2018, Proc. SPIE, 10700, 1070012 [Google Scholar]
  18. Fauvarque, O., Janin-Potiron, P., Correia, C., et al. 2019, J. Opt. Soc. Am. A, 36, 1241 [Google Scholar]
  19. Ferreira, F., Gendron, É., Rousset, G., & Gratadour, D. 2018a, A&A, 616, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Ferreira, F., Gratadour, D., Sevin, A., & Doucet, N. 2018b, in 2018 International Conference on High Performance Computing Simulation (HPCS), 180 [Google Scholar]
  21. Gendron, É. 1995, PhD Thesis, Univ. Paris Diderot, France [Google Scholar]
  22. Gendron, É., & Léna, P. 1994, A&A, 291, 337 [Google Scholar]
  23. Guyon, O. 2005, ApJ, 629, 592 [Google Scholar]
  24. Guyon, O., Martinache, F., Garrel, V., et al. 2010, Proc. SPIE, 7736, 773624 [Google Scholar]
  25. Guyon, O., Sevin, A., Gratadour, D., et al. 2018, Proc. SPIE, 10703, 107031E [Google Scholar]
  26. Héritier, C. T., Esposito, S., Fusco, T., et al. 2018, Proc. SPIE, 10703, 107034P [Google Scholar]
  27. Korkiakoski, V., Vérinaud, C., & Le Louarn, M. 2008a, Proc. SPIE, 7015, 701554 [Google Scholar]
  28. Korkiakoski, V., Vérinaud, C., & Le Louarn, M. 2008b, Appl. Opt., 47, 79 [Google Scholar]
  29. Kulcsár, C., Raynaud, H.-F., Petit, C., Conan, J.-M., & de Lesegno, P. V. 2006, Opt. Express, 14, 7464 [Google Scholar]
  30. Liu, F., & Sanders, G. 2018, Proc., SPIE, 10700, 1070013 [Google Scholar]
  31. Madec, P.-Y. 1999, ed. F. Roddier, Adaptive Optics in Astronomy (Cambridge: Cambridge University Press), 131 [Google Scholar]
  32. Montera, D. A., Brown, J. M., Buckman, M. D., et al. 2018, Proc. SPIE, 10703, 107031H [Google Scholar]
  33. Ragazzoni, R. 1996, J. Mod. Opt., 43, 289 [Google Scholar]
  34. Schatz, L. H., Males, J. R., Close, L. M., et al. 2018, Proc. SPIE, 10703, 1070321 [Google Scholar]
  35. Tamai, R., Koehler, B., Cirasuolo, M., et al. 2018, Proc. SPIE, 10700, 1070014 [Google Scholar]
  36. Vidal, F., Rozel, M., Deo, V., et al. 2019, in Tests and Characterisations of the ALPAO 64×64 Deformable Mirror, the MICADO-MAORY SCAO AIT Facility, 6th AO4ELT Conference [Google Scholar]

Appendix A: Expanding on AO transfer functions

We detail here the transfer functions of the AO loop as a complement to Sect. 2.3 and highlight the differences between a continuous-time description, as is necessary for a real physical system, and a discrete-time description that is adequate for AO simulators such as COMPASS. The material reviewed in Appendix A.1 is well known (e.g., Madec 1999; Kulcsár et al. 2006) but seldom presented, and we deem it necessary for a proper mathematical grounding of Sects. 2.3 and 3. Ultimately, this development lets us derive the values for fcrit and gcrit provided in Table 1.

A.1. Continuous- and discrete-time descriptions

Both the real continuous-time system (CTS) and the simulated discrete-time system (DTS) are clocked at the sampling period T, which is also the integration time of the WFS detector. In the DTS, the turbulence is not averaged over T, but sampled at a single instant. Likewise, the notion of the duration of action of the deformable mirror simply does not exist. Both systems, however, have in common the core of the control: they share the same numerical algorithms as are implemented in the real-time computer, which in both cases are numerical and in discrete time steps.

We introduce several transfer functions hereafter. The WFS performs an integration of the signal during T. In Fourier formalism, its transfer function can be written as

(A.1)

which is the Fourier transform of a rectangular window extending from t = 0 to t = T. We may see Eq. (A.1) as a smoothing by a zero-centered box: sinc(πfT), combined with a half-frame delay, exp(−jπfT), thus ensuring hwfs describes a causal system.

The WFS measurement is processed by the reconstructor matrix to be expressed in a modal space m[k], we drop here the mode index i for clarity, and is integrated with a gain g into the modal command v[k] using the classical discrete integrator,

(A.2)

where k is the loop iteration number. The direct translation of Eq. (A.2) into the frequency space gives the expression of the transfer function of the discrete time integrator,

(A.3)

It should be noted that hdigInt(f) is nonphysical, first because it shows no computational time delay, and second because it is periodically undefined for every multiple of the sampling frequency fS = 1/T: it is only appropriate for representing a discrete time system.

Now, in an real CTS, the RTC computation result v[k] is sent to a digital-to-analog converter (DAC), which holds the (modal) command on the DM during the period T. The DAC, a zero-order hold, has the same transfer function as the WFS,

(A.4)

also including a half-frame delay. With the association hwfs ⋅ hdac, a full frame delay thus naturally appears.

Equation (A.4) joins the digital world back into the physical reality of the CTS by transforming the Dirac comb of numerical samples into an analog signal. Noticeably, the association of the nonphysical numerical integrator with the DAC transfer function form the following product,

(A.5)

that is, the proper transfer function of a continuous-time integrator.

The final point of the discussion is that of the time delay τ. We defined the loop delay as the additional amount of time, on top of the causal process. For a CTS, this is the time spent between the end of the WFS integration and the beginning of the application of the command on the DM. In this case, the transfer function is

(A.6)

When τ = 0, the DM command v[k] exactly applies during the acquisition of m[k + 1]. This one-frame shift is built in the equations thanks to the two half-frames from the WFS and DAC. All the elements are introduced for the transfer function between the input turbulence and the modal measurements mi[k] in the case of a CTS,

(A.7)

For the DTS, there is neither WFS nor DAC windowing, but we need to introduce an apparently artificial one-frame delay exp(−2jπfT) to account for the fact that a command computed at a given iteration can at best only serve for the next one,

(A.8)

Expanding, we obtain

(A.9)

that is, the shorthand h(f; g) of Eq. (3).

With the closed-loop transfer functions with a CTS and a DTS representation of the AO, we now study the critical behavior of the system as the gain increases to approach divergence.

A.2. Critical point. Analog case

The critical frequency and gain (Sect. 3.1) are determined by finding the couple of gcrit and fcrit that zero the denominator of Eq. (A.7). The equation to be solved for g and f is

(A.10)

Equating for phase and modulus,

(A.11)

which reduce to

(A.12)

(A.13)

A.3. Critical point. Digital case

In the same way as for Appendix A.2, the equation to be solved is

(A.14)

Leveraging

we obtain the phase and modulus equations,

(A.15)

and ultimately,

(A.16)

(A.17)

which is the general formula given in Table 1. We finally note that for a CTS with latency τ ≥ 1, we have fcritT < 1/6 ≪ 1, and therefore Eqs. (A.13) and (A.17) both reduce to the same value gcrit ≈ 2πfcritT in a first-order approximation. The general discussions and derivations presented in Sect. 3 are in this case also applicable to a continuous-time description of the system. As a final note, Kulcsár et al. (2006) described the necessary hypotheses and proved that a discrete-time treatment of the minimization criterion in Eq. (2) is equivalent to a continous-time treatment.

Appendix B: Semianalytical CLOSE solutions for latencies of 0 and 1 frames

We propose in this appendix some additional data to expand the discussion in Sect. 3.3 by exploring the close match obtained for all sensitivity-corrected S/N σi between the minimum variance gain gMV and the steady-state gain yielded by CLOSE for a setpoint r = 0. While in Fig. 6 this comparison was restricted to a latency τ = 2T, we show in Fig. B.1 the same analysis for latencies τ = 0 and τ = T. All parameters are otherwise similar, in particular, the sample atmospheric spectrum (Eq. (1)) with a cutoff of 1 Hz. The time difference δkcrit at which the autocorrelation is locked to the setpoint r is taken from Table 1.

thumbnail Fig. B.1.

Minimum variance gain gMV and solutions found using CLOSE with five different setpoint values, depending on the sensitivity corrected S/N σi. These computations are performed using the modal input spectrum shown in Fig. 3.

Open with DEXTER

For all three latencies τ = 0, T, and 2T, we observe the same remarkable match between the minimum variance solution and CLOSE for a null setpoint. This was also simulated for integer latencies up to τ = 5T, although these cases are not necessarily relevant for a real AO system and are not reported here. While Fig. B.1 (top) reports a very satisfactory correspondence, the analysis we proposed reaches the limit of the discrete time approximation of real AO systems, as we showed in Appendix A. With fcrit = fS/2, the case τ = 0 is therefore only adequate for simulated systems.

As commented on in the main text, the same conclusions apply at other latencies: Positive setpoints induce an undersetting of the integrator gain, while negative setpoints produce an oversetting and impose a gain floor at low S/N. This concludes our analysis of the minimum-variance capability of CLOSE for various latencies. We are confident that this property would extend in between, to fractional latencies within that range and beyond.

Appendix C: Experimenting with the correlation-locking setpoint – additional data

In this appendix we propose supplementary simulations to those in Sect. 5.3 to investigate the final long-exposure SR that is produced in the system imager versus the setpoint r used to drive CLOSE. The experimental protocol is identical to that in Sect. 5.3, generalized to three different latency values τ/T = 0, 1, and 2 frames for three different AO system designs: the MICADO parameters used in the main text, and a PWFS and a Shack–Hartmann, both on an 8 m diameter telescope. The parameters for these simulations, thereafter referred to as PWFS8M and SH8M, are given in Table C.1 or are identical to those in Table 2 if unspecified.

Table C.1.

AO numerical simulation parameters for simulation setups PWFS8M and SH8M.

Results with varying latencies, guide-star magnitudes, seeing conditions, and wind speeds against the r setpoint are given in Tables C.2C.4, which also include results presented in Sect. 5.3 (see Table C.2, “Latency τ = 2.T”). In the tables presented here, we provide the maximum SR in H band for each given set of conditions [AO setup, τ, r0, MR, wind speed V] that was achieved over r values from −0.5 to 0.5. This value is noted SR(rmax), and the argument at which it is attained is noted rmax. For rmax ≠ 0, we provide the value that was attained for the same parameters at our favorite all-in-one solution r = 0. This value is noted SR(r = 0).

Table C.2.

Summary of long-exposure SRs in H band for stationary simulations exploring the latency, seeing conditions, guide-star magnitude, and wind speed against the CLOSE r setpoint for the MICADO simulation setup.

Table C.3.

Identical to Table C.2 for the PWFS8M simulation setup as defined in Table C.1.

Table C.4.

Identical to Table C.2 for the SH8M simulation setup as defined in Table C.1.

Conclusions for the MICADO simulations (Table C.2) are essentially identical to the analysis provided in Sect. 5.3: minor fluctuations in the rmax value, with r = −0.1 or r = 0 being almost always optimal within the LE SR error bar of a few percent. These simulations are therefore conclusive in showing the versatility of CLOSE for different system latencies without additional complexities other than changing the value of δkcrit.

The simulations performed with the PWFS8M setup (Table C.2) show very similar trends to the MICADO setup, which indicates that our conclusion that r = 0 almost always optimally works may very well apply to a wide range of PWFS designs. A performance improvement and the automation of the integrator tune-up can also be a benefit for AO systems smaller than an ELT. The PWFS8M simulations at τ = 0, MR = 12 allow us to note the broadening of the bell function of the SR as a function of r, as was visible in Fig. 10. Some rmax values are shifted toward positive values, up to +0.3, but the difference between SR(rMax) and SR(r = 0) is negligible. It is lower than 1%. This behavior is expected because the margin for setting αiGi, thus r, decreases with increased latency, as the trade-off between low frequency rejection and noise amplification becomes tighter. Thus, for very short latencies, the precision of αiGi for stable control and variance minimization is of little importance, and it follows that the setpoint r also benefits from an increased margin.

The simulations with the SH8M setup reported in Table C.4 are comparable to PWFSs, except for the brighter cases MR = 8. For these, and most particularly for MR = 8, τ = 0, the ideal setpoint rmax is noticeably pushed to positive values, even hitting the extremum of our parameter space r = +0.5, although based on the recorded SRs, we may safely assume that no significant improvement would be reached for r in the 0.5–1.0 range. However, the difference in performance between rmax ≥ 0.5 and r = 0 is again insignificant, except for the cases with a wind speed of 40 m s−1.

The cases discussed here (SH8M, τ = 0, MR = 8) are the one simulation in which modal optimization is least important, if at all: the sensor is linear, the S/N is extremely good, and the bandwidth effect is minimum. The need for balance between noise amplification and turbulence rejection is quite different from the design case explored in the main text. For all other cases with the SH8M setup, that is, where variance minimization by the rejection or noise trade-off is valid, we validated the near-minimum variance performance of CLOSE with a null setpoint, extending the conclusions seen for PWFSs. In particular, this also implies that the method we proposed is sound beyond certain hypotheses of the numerical analysis of Sect. 3, and in particular beyond the spectrum proposed in Eq. (1). The power law in Eq. (1) ought to be changed to −11/3 when accounting for the aliasing with a Shack–Hartmann sensor (Conan et al. 1995).

All Tables

Table 1.

Parameters related to transfer function divergence depending on the latency of the system.

Table 2.

AO numerical simulation parameters.

Table C.1.

AO numerical simulation parameters for simulation setups PWFS8M and SH8M.

Table C.2.

Summary of long-exposure SRs in H band for stationary simulations exploring the latency, seeing conditions, guide-star magnitude, and wind speed against the CLOSE r setpoint for the MICADO simulation setup.

Table C.3.

Identical to Table C.2 for the PWFS8M simulation setup as defined in Table C.1.

Table C.4.

Identical to Table C.2 for the SH8M simulation setup as defined in Table C.1.

All Figures

thumbnail Fig. 1.

General modeling of the SCAO loop showing the WFS, real-time computer, and DM. The WFS is modeled using the confusion matrix framework (Sect. 2.1). The confusion matrix A is a random variable depending upon the residual wavefront ϕRes.

Open with DEXTER
In the text
thumbnail Fig. 2.

Diagonalization of the AO loop model (Fig. 1) for mode i under the approximation that A is diagonal and is thus reduced to its ith diagonal term αi.

Open with DEXTER
In the text
thumbnail Fig. 3.

Typical spectral components for our AO semianalytical computations. Top: input spectra for turbulence and noise, and the measurement power transfer function |hi(f; αiGi)|2 for three αiGi values. Bottom: corresponding measurement power spectra . For this example, τ = 2T, and fS = 500 Hz, yielding gcrit ≈ 0.61 and fcrit = 50 Hz.

Open with DEXTER
In the text
thumbnail Fig. 4.

Autocorrelation functions corresponding to the measurement PSDs shown in Fig. 3. δkcrit = 5 frames for .

Open with DEXTER
In the text
thumbnail Fig. 5.

For effective S/N σi = 1, 10, and 100, closed-loop residual variances (top panels) and resulting δkcrit autocorrelation of the measurements (bottom panels) as a function of the loop gain αiGi from 0 up to gcrit. The latency is taken as τ = 2.T. Black lines show the total residual variance normalized to 1 at its minimum, blue represents the variance induced by turbulence , and green shows the variance induced by noise .

Open with DEXTER
In the text
thumbnail Fig. 6.

Minimum variance gain gMV and solutions found using CLOSE with five different setpoint values, depending on the sensitivity-corrected S/N σi. These computations are performed using the modal input spectrum of Eq. (1) with fi = 1 Hz.

Open with DEXTER
In the text
thumbnail Fig. 7.

Convergence of CLOSE gains on the 2.0 s following closing of the AO loop for guide stars of MR = 0 and 16 and atmospheric r0 of 14.5 and 9.0 cm. All Gi are initialized to 0.5 at the start (blue line). Gi values are shown as averaged over the time windows given in the legend. Curves are smoothed along the i index for clarity. Final SRs are given in H band and are computed from the cumulative exposure over the last 200 ms (k ∈ [900, 1000]).

Open with DEXTER
In the text
thumbnail Fig. 8.

Time series of CLOSE gains for four select modes in the 2.0 s following closing of the AO loop for the simulation r0 = 14.5 cm, MR = 0 shown in Fig. 7.

Open with DEXTER
In the text
thumbnail Fig. 9.

Long-exposure SR in H band obtained for end-to-end simulations of the MICADO SCAO setup with CLOSE (top), and with a manually optimized scalar integrator (bottom) for guide-star magnitudes MR = 11.5–17.5 and seeing conditions r0 = 8.9–21.5 cm.

Open with DEXTER
In the text
thumbnail Fig. 10.

Long-exposure SR in H band (color and text) for stationary simulations exploring different conditions regarding the seeing (outer rows), the guide-star magnitude (outer columns), the wind speed (inner rows), and the CLOSE setpoint (inner columns) for the MICADO SCAO simulation. Color scales are local to each of the subplots and do not match each other.

Open with DEXTER
In the text
thumbnail Fig. 11.

Top: time series of the H-band SR and modal gains Gi for a few select modes during the simulation of a burst of seeing; r0 = 21.5 cm for t < 2 s and t > 6 s and r0 = 8.9 cm for 2 < t < 6 s. Curves are smoothed using a 10 ms window. Bottom: snapshot of the modal gain vector G at time t = 2, 6, and 8 s.

Open with DEXTER
In the text
thumbnail Fig. 12.

Similar to Fig. 11. Time series and modal gain snapshots during a cloud-like transient, simulated as a guide-star magnitude change: MR = 14 for t < 2 s and t > 6 s and MR = 17 for 2 < t < 6 s.

Open with DEXTER
In the text
thumbnail Fig. B.1.

Minimum variance gain gMV and solutions found using CLOSE with five different setpoint values, depending on the sensitivity corrected S/N σi. These computations are performed using the modal input spectrum shown in Fig. 3.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.