First on-sky SCAO validation of full LQG control with vibration mitigation on the CANARY pathfinder

: Adaptive optics provides real time correction of wavefront disturbances on ground based telescopes. Optimizing control and performance is a key issue for ever more demanding instruments on ever larger tele-scopes affected not only by atmospheric turbulence, but also by vibrations, windshake and tracking errors. Linear Quadratic Gaussian control achieves optimal correction when provided with a temporal model of the disturbance. We present in this paper the ﬁrst on-sky results of a Kalman ﬁlter based LQG control with vibration mitigation on the CANARY instrument at the Nasmyth platform of the 4 . 2-m William Herschel Telescope. The results demonstrate a clear improvement of performance for full LQG compared with standard integrator control,


Introduction
Adaptive optics (AO) is a technique that allows a deformable mirror (DM) to perform realtime compensation of the wavefront (WF) distortion induced by the atmospheric turbulence encountered during astronomical observations.AO has been successfully demonstrated and made available for astronomers since the early 1990's [1].Its success is such that the next generation of Extremely Large Telescopes (ELTs) cannot be conceived without it.The associated future astrophysical programs rely on sophisticated AO systems able to fulfill ever demanding requirements.These systems can be globally split into two categories.On one hand wide-field AO (WFAO) uses a large field for wavefront sensing and/or for correction.On the other hand single conjugated AO (SCAO) uses an on-axis guide star both for wavefront sensing and correction.Extreme AO corresponds to very high performance SCAO feeding specific high contrast instruments.
The present paper is focused on SCAO, with on-sky experiments conducted using the SCAO mode of the CANARY demonstrator [2].This demonstrator has been initially designed to validate a specific WFAO concept: multi-object AO (MOAO), required for the study of primordial galaxies, where several objects of small angular dimensions and spread over a large field of view are observed simultaneously [3].The CANARY MOAO pathfinder is the demonstrator for the future MOAO instrument EAGLE [4] and the first on-sky demonstration of MOAO was performed successfully in September 2010 [5].
Control and performance optimization of complex AO systems are a key issue.Linear Quadratic Gaussian (LQG) control is an appealing strategy for AO since it provides globally optimal control with respect to an explicit minimum variance performance criterion [6,7].Numerical studies [8] and laboratory demonstrations [9] have shown that LQG control can give better performance than integral controllers due to the use of spatial and temporal priors.Theoretical studies [10,11] have also shown that an advantage of LQG controllers over integrator was the possible relaxation of system constraints (actuator saturation, sampling frequency, noise requirements) which may be of particular interest for ELT systems.A major issue in AO is also the correction of non-turbulent perturbations such as windshake or vibrations [12][13][14][15].LQG control, thanks to its flexibility, can provide a very efficient compensation of these effects, as shown for vibrations through laboratory tests [16], and through numerical simulations for windshake [17].
LQG control has been implemented for tip-tilt correction on the SPHERE instrument, the European extrasolar planet finder [18].It has also been considered during preliminary studies on GPI using a Kalman filter with predictive Fourier-domain control [19].LQG control has been demonstrated on-lab for WFAO systems [20].Note that tip tilt optimal control has been tested on-sky using an AO system installed at the McMath-Pierce solar telescope [21].The next step was to demonstrate full LQG control on-sky, full meaning with optimal LQG control of all modes, and not only tip and tilt.
This paper presents a detailed analysis of the first SCAO on-sky validation of a full LQG control with vibration mitigation, performed in July and August 2012 [22,23].We compare performance obtained with integrator versus LQG, with and without vibration mitigation.Vibration models for tip and tilt modes are identified on line from WFS data, using the frequency-domain identification method selected for the instrument SPHERE [17], hence validating on sky and for the first time its tip-tilt identification and control strategy.Analysis of models relevance is conducted, and the advantage of LQG control is also discussed on the basis of power spectral densities and rejection transfer functions.
We organize the paper as follows.In Section 2 we briefly present the CANARY bench and its SCAO mode.Section 3 introduces the formalism of the LQG control for SCAO.In Section 4 we discuss practical implementation issues such as model identification or real-time complexity for the CANARY instrument.Then in Section 5 we present the first on-sky results with a detailed data analysis.This is followed in Section 6 by a frequency-domain and a robustness analysis for integrator and LQG, and Section 7 concludes this paper.

The single conjugated adaptive optics mode of CANARY
CANARY is a European project led by Durham University and Observatoire de Paris.The CANARY bench is set up at the Nasmyth focus of the 4.2-m William Herschel Telescope at the summit of the Roque de los Muchachos mountain on La Palma, Canary Islands, Spain [24,25].CANARY can perform AO in several possible configurations including SCAO, MOAO using a four Natural Guide Stars (NGSs) asterism (three off-axis stars and one on-axis for performance evaluation and scientific IR images) and laser based MOAO with four additional Laser Guide Stars (LGSs) in the technical Field of View (FoV).In this paper we focus on the SCAO mode of CANARY.In SCAO, an on-axis WF Sensor (WFS) (so-called Truth Sensor (TS)) provides closed-loop 7 × 7 Shack Hartmann WFS data on the on-axis star.Correction is performed by a high-speed tip-tilt mirror developed by Observatoire de Paris, and a CILAS 8 × 8 piezostack array DM conjugated to the telescope's pupil.The CANARY Real-Time Computer (RTC), called DARC, is developed by Durham University [26,27] and IR images are recorded thanks to the near-IR camera CAMICAz, developed by Observatoire de Paris.CAMICAz has a NICMOS-3 array composed of 256 × 256 pixels.The camera has several spectral filters.We observe in the two spectral bands H and K.The central wavelengths are λ H = 1.67µm and λ K = 2.18µm with ∆ λ H = 0.14µm and ∆ λ K = 0.15µm.The plate scale of the IR camera is p = 28.8mas.A telescope simulator developed by UK ATC allows performance evaluation on internal sources with emulated turbulence being introduced by rotating phase plates.

LQG control for SCAO
LQG control can be derived when the system is described by a linear state space model.This approach provides a general framework based on explicit models of the AO system including wavefront perturbations [8].We present here a full LQG on-sky SCAO demonstration, i.e.LQG control is applied on all modes (tip, tilt and higher orders).For tip and tilt modes, vibratory components are accounted for as in [16,17].
In this Section we explain how to establish the dynamical models including vibrations, and we derive an optimal control based on a criterion that simplifies system calibrations.We first present the minimum variance performance criterion and the optimal solution in Section 3.1.We then present the modeling of the disturbance stochastic process (Section 3.2) with a specific treatment for tip-tilt components (Section 3.3).Section 3.4 establishes the corresponding SCAO state space representation on which is based the Kalman filter used for control computation (Section 3.5).

Performance criterion and optimal solution
The perturbation phase is a continuous time process.Since the DM is applied through a zero order hold and the DM's time response is small compared to the loop sampling period T , the continuous time minimum-variance criterion can be equivalently replaced by a discrete time criterion for control optimization [6].We briefly recall here this construction in the simple 'synchronous' case where controls are applied at the same sampling times t = kT , and where one sampling interval is allocated for WFS measurement processing and control calculation.In this case, the control objective is therefore to minimize the following discrete-time infinitehorizon quadratic performance criterion J(u), which corresponds to the variance of the residual where the correction phase φ cor k is related to the command vector u through and where u k−1 is the control vector applied at time (k − 1)T , N being the DM influence matrix.
In formula (1), the discrete-time perturbation φ k is the averaged perturbation over the sampling interval [(k − 1)T, kT [, defined as All phase variables are defined similarly as temporal averages over one frame of duration T .When the phase is described by Zernike coefficients, the influence matrix N is the operator that connects command vector u to Zernike modes.The optimal solution u * k that minimizes J(u) is given through an LQG approach if these two hypothesis are satisfied: • the system can be described by a linear state-space model with additional noise inputs; • the noises are white and Gaussian, i.e. uncorrelated sequences of Gaussian vectors; covariance matrices are known.
More details on LQG control can be found in [28].When the control u k is computed from measurements y k , . . ., y 0 , the stochastic separation principle allows to solve this problem in two steps: • compute the estimated phase as the conditional expectation of φ k+1 knowing all measurements y k , . . ., y 0 available at time step k: • project the estimated phase onto the DM: where P u is the projector defined as pseudo-inverse of N: The optimal prediction φ k+1|k is computed recursively as the output of the Kalman filter based on a state-space representation of the complete system (see Section 3.5).The measurement equation is defined as which can be written using Eq.(2) as where D is the WFS matrix that characterizes the WFS response in terms of the input phase, and w is the measurement noise supposed to be white Gaussian and independent of φ , with covariance matrix Σ w .Note that this measurement equation features a two-frame loop delay between DM control u and WFS output y, also called 'plant delay' in control literature.This loop delay induces an important limitation of the closed-loop performance.
A difficulty may arise for the calculation of the projector N † in Eq. (6).It requires the knowledge of the influence matrix N, which depends on the actual system geometry.N is not only related to the DM influence functions but also to the relative position of the actuator grid with respect to the pupil.Moreover it has to be oriented with respect to the WFS geometry.Rotations, magnifications, and misregistrations may occur between these grids which imply specific calibration procedures [29,30] that can also be avoided by modifying the control criterion J(u) as explained below.
Since the interaction matrix defined as M int DN is easily and regularly recorded on the bench, we propose to minimize the alternative quadratic criterion L defined as This modified criterion is the residual slopes variance as opposed to the residual phase variance J(u).Minimizing L (u) instead of J(u) leads to the same optimal phase estimate and the same control structure, but with a different projector onto the DM, namely: Computing this new projector only requires the knowledge of the WFS matrix D and of the command matrix M com = (DN) † , i.e. the pseudo-inverse of the command-to-slopes interaction matrix M int = DN.Although sub-optimal with respect to the residual phase variance criterion J(u), this control has the advantage to be very robust to system mis-alignments and avoids tricky calibration procedures for the determination of N. The relative geometry between WFS and DM is indeed accounted for in M int , therefore avoiding an explicit determination of magnification, rotation and decentering between WFS subaperture grid and DM actuator grid.Note that the same 'slopes-oriented' approach lies at the heart of the APPLY MOAO control strategy, where open-loop on-axis WFS slopes are reconstructed from off-axis measurements using a static minimum-variance estimator [31].
The WFS matrix D is needed to compute the slopes-to-control projector P u .Because it appears in the measurement equation (8), it is also needed to compute the Kalman filter presented in Section 3.5.In practice, D was computed theoretically from a geometrical model and scaled using straightforward calibrations done on internal sources, a procedure which enables to determine both the sensitivity of the WFS and the ordering of the subapertures.The interaction matrix M int = DN was also calibrated on the bench.

Disturbance stochastic dynamical model
As explained above, the computation of the optimal prediction φ k+1|k relies on a state-space representation.It can be established when a stochastic dynamical model of the perturbation is available.Tip and tilt (TT) are the most energetic turbulent modes and frequently include additional perturbations due to vibrations in the telescope and system structure, windshake and residual tracking errors.These factors led us to separate TT modeling from that of the higher orders (HO).
Vibrations are additive perturbations and independent from the turbulence, so that the complete perturbation phase can be written as: Assuming that vibrations only affect the tip and tilt modes, the TT and HO components of the disturbance φ , written as can be decomposed as φ TT = φ tur,TT + φ vib,TT φ HO = φ tur,HO .
In the case of CANARY we have checked that there were no noticeable vibration peaks on temporal PSDs for modes above tip and tilt, hence the restriction of the specific treatment to these two modes.Note however that the above equations can be easily extended to the case where vibrations may affect higher orders than tip and tilt, as expected on large structures such as ELTs.In this case, the superscript TT has just to be replaced with LO (for low orders) and the corresponding vectors will contain as many components as needed.
A widely used class of models to describe the phase temporal correlation is the class of vector-valued autoregressive (AR) processes.First studies have shown that an AR model of order 2 (AR2) for each turbulent phase mode leads to good performance, see [32] and in the case of CANARY see [33].We can thus model the turbulence process as a discrete-time stochastic process in the form where A tur 1 and A tur 2 are diagonal matrices and v tur {v tur k } k≥0 is a generating Gaussian white noise.

Tip-tilt modeling
The vibrations affecting the system can also be modeled by an AR2 [16] (resonant order 2 model).Tip and tilt are treated separately and considered as decoupled.For either tip or tilt components each particular vibration peak will be associated with a scalar vibration model in the form where φ vib,•,i k is a scalar variable and • stands for tip or tilt.The scalar parameters a vib,•,i 1 and a vib,•,i 2 depend on damping coefficient and resonant frequency, and v vib,•,i {v vib,•,i k } k≥0 is a white Gaussian excitation noise with variance σ 2 v •,i .The global resulting phase in the pupil can be written as: with n v the total number of vibrations on tip or tilt.As explained in the introduction, these vibration models will be identified on-sky from WFS data (see Section 4.1).
As already mentioned, the extension to vibrations on other modes than tip and tilt is straightforwardly obtained by considering that • stands for all the concerned modes, and with superscripts tip and tilt replaced by a mode number.The right hand term in Eq. ( 17) will therefore have as many lines as the number of modes affected by vibrations.

SCAO state-space representation
The state vector can be simply chosen as the concatenation of turbulence and vibration phases at two consecutive time steps: (18) The models presented above can then be re-arranged into a compact state-space model in the form: The last line in Eq. ( 19) expresses the perturbation as a sum of contributors stored in the state x.
Both the WFS measurement y k and the disturbance φ k are thus obtained as outputs of the state model.The other vectors and matrices in the state space representation are defined as where v {v k } k≥0 and w {w k } k≥0 are independent Gaussian white noises with variancecovariance matrices Var(v k ) Σ v and Var(w k ) Σ w , and where To obtain the perturbation as the sum of turbulence and vibrations, the matrix C φ is defined as with Likewise, to obtain Cx k = Dφ k−1 , the matrix C should be defined as with

Kalman filter and control computation
The optimal prediction x k+1|k , which enters the optimal control expression in Eq. ( 5), is computed recursively as an output of the Kalman filter based on the state-space model (19).Since the control criterion is defined on an infinite horizon, we restrict ourselves, without any performance degradation, to the steady-state time invariant version of the Kalman filter [34].
Real-time implementation boils down to the so-called update and prediction equations: where H ∞ is the asymptotic Kalman gain computed off-line as Σ ∞ being the unique solution of the algebraic Riccati equation Finally we are able to compute the control vector applied to the DM: = P u ( φ tur k+1|k + φ vib k+1|k ) where P u depends on the control criterion: Eq. ( 6) for residual phase variance and Eq.(10) for residual slopes variance.The control gain that relates the state estimate xk|k to the optimal control is thus simply P u C φ A.

Model identification for turbulence and vibrations
We describe in this Section the identification procedures used to derive the parameters of our AR2 models.As explained before, we apply a separate treatment for the TT and the HO: • TT components are identified thanks to the method presented in [17], which allows identification of a TT turbulent model together with vibration components.This method uses open-loop slopes to compute the mean slopes PSD.The superimposition of several AR2 models is fitted to this PSD: a low frequency AR2 model for turbulence, and n v vibration components.For each identified vibration the method estimates the resonant frequency, damping coefficient and power (i.e.variance of the excitation noise).The number of vibration components n v to be identified is set to 10 both for tip and tilt.Note that pseudo-open loop slopes could be used as well if residual slopes only were available.
• As mentioned above, HO modes are not affected by vibrations, so that their model can be based on the heuristic formula in [8,35]  PSD structure [36].It relates the AR coefficients of the i th mode to wind speed norm and radial order and is an extension of the method proposed in [8]: where V wind is the average wind speed norm value in the pupil, D pup is the pupil diameter, n(i) the Zernike radial order of mode i, ξ is the damping factor set to 0.9, T is the sampling period.Note that the model noise variance σ 2 v i is derived from the HO covariance matrix deduced from the Fried parameter.Wind speed and Fried parameter are estimated thanks to the LEARN algorithm [31,37].The maximum radial order has been chosen equal to n rad = 14.

Complexity, real-time control and delays
From Eq. ( 18), the state vector x includes the values of all phase variables at two consecutive time steps.The phase value at each time step is expressed using n rad (n rad + 3)/2 = 119 Zernike modes.When filtering vibrations with 10 vibrations peaks both on tip and tilt, the size of the state x is therefore n x = 2 × (119 + 2 × 10) = 278.The number of measurements is equal to n y = 36 × 2 = 72.We have 52 actuators on the DM and 2 actuators for the TT mirror, so that the size of u is n a = 54 components.To speed computation, the LQG implementation of the DARC RTC has been optimised to remove unnecessary operations (such as products of zero or identity matrices).Consequently, the number of floating multiplications performed in real time is equal to n x /2 + n 2 x /2 + 3n x n a /2 + n x n y ≈ 52k per frame for this configuration.The Kalman gain computation, which centrally involves solving the Riccati equation (28), is not part of the RTC data pipeline and is performed off-line on a separate interface computer, in our case a standard desktop device.For a state vector with 278 components, solving the Riccati equation took only about 1 to 2 seconds using the efficient doubling algorithm [28].
Both on-and off-line computation times are not a limiting factor for routine operations, and the very flexible general formulation given in Eqs. ( 26)-( 29) makes it advantageous.Of course when moving to more complex systems for VLTs, or even ELTs, both off-line and on-line calculations may become an issue.However various approaches have been recently proposed to overcome this problem by exploiting the particular structure of the operations (linear shift invariance, sparsity, structure...) [38][39][40][41], and RTCs' computing power keeps increasing.
Another important practical and theoretical issue is the influence of loop delay on closed-loop performance.In the system model presented in Section 3, we make the standard assumption that the total loop delay between DM control u and WFS output y in the measurement equation ( 8) is equal to two frames.However, during on-sky operations the actual system loop delay d was measured at around d = 1.6 frames for a loop frequency of 150 Hz (the value usually selected in order to achieve adequate signal-to-noise ratio for a 10-12 magnitude star at visible wavelengths).In all the simulations and replays presented in this paper, the system loop delay was therefore set to d = 1.6.With integrator control, it is easily seen that lower delay means better disturbance rejection (in addition, as we shall see, it also renders integrator performance considerably less sensitive to changes in the controller gain).The LQG controller is based on an explicit model of the AO system, including delays.In our case, it has been designed using a model with a loop delay of two frames, which simplified its implementation in DARC RTC.This creates a system/model mismatch which may lead to a slight degradation of LQG performance.However, in order to be consistent with on-sky implementation, the value of two frames for the model loop delay in the LQG controller has been kept in all what follows.

On-sky results
We present in this Section the first results obtained on-sky using an LQG control on all modes.Taking advantage of the CANARY RTC's flexibility, LQG control (with and without vibration mitigation) and standard integral control were alternatively put on line over successive periods, during which sequences of residual slopes, open-loop slopes and PSF images were recorded.This enabled to compare the performance of these three controllers (LQG with and without vibration compensation and integrator) in the most possible similar on-sky conditions.However, sudden and important variations of the Fried parameter r 0 regularly occur, which makes fair comparisons of performance series often difficult.We have thus preferred to compare performance as a function of r 0 , estimated as further explained.These SCAO results have been obtained during the nights of August 28 th and 29 th 2012 (time is indicated hereafter in GMT+1, La Palma local time) on the star BD + 10 3697 of magnitude m V = 9.7.
The WFS slope noise estimated on the slope data is of the order of 80 nm rms, expressed as an optical path difference at the edge of a WFS subaperture.This noise is mainly a consequence of photon noise at the level of the WFS subaperture images detection, since we use a low detector noise EMCCD ANDOR camera.Its impact on performance is difficult to estimate precisely and depends on the choice of the controller.We still estimate in our case that an upper bound of its contribution is 10% of the global residual phase variance in H band.

Integrator
The control computation for the integrator is given by: where g int is the global feedback gain (i.e., the same gain was used for all modes).This gain has been optimized on-sky for each night of the CANARY runs.In Fig. 1, we give an example of on-sky determination of the gain.We can see that performance is not very sensitive to the gain value.For this particular night, the gain was set to 0.6.This flatness of the on-sky performance sensitivity seems counterintuitive, since one would normally expect a parabolic-like performance curve.However, such parabolic performance would correspond to a total loop delay of two frames, whereas for d = 1.6 the expected performance becomes essentially flat (more on this in Section 5.3).

LQG
A convenient tuning parameter in LQG control is the measurement noise variance Σ w taken in the model.Performance is generally improved when setting a noise level higher than the physical noise values so as to absorb various model uncertainties inherent to experimental operation [9,20,33].Prior to starting on-sky operation, this relationship was confirmed through investigation using the CANARY telescope simulator in similar noise conditions to that expected on sky.An amplification of the noise variance by a factor of 2 improved performance and thus this tuning value was used for subsequent on-sky tests.

Performance results
We first present overall performance in terms of Infra-Red (IR) Strehl Ratio (SR), and then focus more specifically on the TT correction.Figure 2 presents the SRs obtained on August 28 th 2012 in H-band and on August 29 th 2012 in H and K-band.We compare in each case integrator control (gain set to 0.6), LQG control with no specific handling of vibrations and LQG control with vibration filtering.Despite strong dispersion in seeing conditions, we clearly see that LQG with vibration mitigation indeed brings a gain in performance compared to LQG without vibration mitigation.Moreover, the LQG controller behaves as expected, globally better than integrator, in terms of performance and stability.
The turbulent priors used in the model building for LQG control (V , r 0 , Σ w ) are identified thanks to the LEARN algorithm [31,37] which uses an open-loop WFS data sequence of 5000 samples.At a sampling frequency of 150 Hz, this corresponds to 33 seconds of acquisition.These data sets are also used for the identification of the TT models.We record a new set of such open-loop data every time we change configuration or asterism on the sky to re-estimate all priors and parameters.LQG matrices were updated in this fashion about every hour.
The fact that performance is not very sensitive to r 0 may seem surprising at first glance.However a gross estimate of error contributors for the integrator shows that more than half the residual phase variance is related to non-seeing dependent terms (Fig. 3).These terms are mainly due to uncorrected bench static errors and to vibrations.Figures 4 and 5 show some examples of Point Spread Functions (PSFs) obtained with the 3 control strategies (integrator, LQG with and without vibration mitigation).As the controller becomes more efficient, SR improves, the PSF becomes more symmetric and first Airy ring more contrasted.

Performance analysis in replay mode
These on-sky comparisons between LQG and integrator controller raise two important questions.First, would LQG retain its edge when compared with a properly tuned Optimized Modal Gain Integrator (OMGI), where different feedback gains are used for different turbulent modes in lieu of the global gain g int ?A second relevant question is whether the gap in performance between LQG and integrator is mainly due to vibration peaks and their particular localization, in which case the said performance gap might just as well vanish altogether for slightly different on-sky conditions.In order to address these issues, we exploited the fact that during CANARY experiments a sequence of open-loop WFS slopes is systematically recorded every time a control is tested.To evaluate the performance that any other controller would have achieved on-sky at this particular moment, it therefore suffices to feed this open-loop WFS data sequence in a filter that emulates its closed-loop rejection transfer function.This produces sequences of residual slopes, from which we extracted residual phase trajectories by applying a minimum-variance slopes-to-Zernike modes reconstructor.This so-called 'replay mode' enabled to replicate with excellent accuracy the actual on-sky performance of both LQG and integral controllers.Furthermore, replay mode simulations reveal that integrator performance is highly sensitive to changes in the system loop delay d. Figure 6 presents the performance variation as a function of the gain, obtained by replaying the data in two cases: d = 1.6 frame and d = 2 frame delay.The same behaviour has been obtained on all data sets: the 1.6 frame delay curve is essentially flat, which is consistent with on-sky results in Fig. 1.

Optimized modal gain integrator performance in replay
As already mentioned, the same gain has been used for all modes in on-sky integrator tests.To explore in replay the benefit of choosing a different integrator gain for each mode, Zernike modes were reconstructed from open-loop on-sky WFS data; we then applied to each sequence the integrator rejection transfer functions corresponding to values of g int ranging from 0.2 to 0.8 and draw performance curves in terms of residual phase variance (in rad 2 ) for each reconstructed mode.This detailed modal analysis showed that, due to the presence of vibrations affecting only tip and tilt and to the moderate WFS noise level, best OMGI performance could be achieved by using only two different gains: a lower value g TT for tip and tilt modes and a higher value g HO set identical for all higher-order modes.
Figure 7 shows representative performance levels (in rad 2 ) obtained on the same open-loop data-set (August 29th at 22:07) for this gain tuning strategy.The blue curve is the performance of the integral controller when tip and tilt gains are set to g TT = 0.3 while the higher orders gain varies from g HO = 0.2 to g HO = 0.8.The red curve corresponds to another tuning, where the higher orders gain is set to g HO = 0.6 while the tip and tilt gain g TT varies from 0.2 to 0.8.Note that for g TT = 0.3 the performance curve exhibits a minimum around g HO = 0.6 and is almost flat for 0.5 ≤ g HO ≤ 0.8.Thus, a finely tuned OMGI controller would give roughly the same performance as g TT = 0.4 and g HO = 0.6 -a quasi-optimal setting which corresponds to the lowest point on the red curve.The on-sky setting for the global integrator gain during this night was g int = 0.6, which corresponds to a quasi-optimal choice of g HO but a markedly sub-optimal choice of g TT .It also appears that using this setting instead of the optimal tuning would increase the residual variance by about 5%, whereas LQG residual phase variance (in magenta) is 25% lower than this best integrator performance.This trend was confirmed on the many data sets recorded during different nights.

Vibration mitigation/amplification and performance
Comparing on-sky performance of different controllers is difficult because disturbance conditions vary over very short time periods.Replay mode offers a neat way to circumvent this difficulty, since it enables to compare the performance that different controllers would achieve when confronted to exactly the same disturbance trajectory.To assess the impact of vibration on closed-loop integrator and LQG performance, we thus present here results obtained by replaying an open-loop sequence of WFS measurements acquired on the night of August 29th with the 3 controllers (integrator, LQG with and without vibration compensation) that had been tuned and tested on-sky at that point in time.For LQG control, this means that we have used here the matrices that have been computed during on-sky operations.
From this sequence of on-sky open-loop slopes and these 3 sequences of residual slopes obtained in replay (one for each controller), we extracted tip and tilt modes trajectories using a modal reconstructor (pseudo-inverse of D) and computed the corresponding 'tip plus tilt' cumulated PSDs (i.e. the sum of tip and tilt cumulated PSDs at each frequency).These cumulated PSDs, which thus provide for each controller a frequency-by-frequency total tip and tilt residual variance, are presented in Fig. 8.The open loop curve presents several steps indicating vibration peaks at different frequencies.At frequencies up to a few hertz they are strong but have very limited impact on performance since all controllers have a good rejection (see Section 6.2).The last vibration peak, located at 22 Hz, is more critical since it is amplified by the integrator control (see Section 6.2 for details) and not corrected by standard LQG control.However, LQG with dedicated vibration identification and filtering mitigates this contribution.On these data we estimate that the 22 Hz vibration contribution amounts, for the integrator control, to a loss of SR in K band by a factor 0.9 .
For results such as the ones in Fig. 8, it is clear that the vibration peak affecting tip and tilt modes around 22 Hz is amplified by the integrator and mitigated by the LQG controller.To make a quantitative assessment of the impact on overall performance for this particular on-sky dataset, we have evaluated the corresponding cumulated open-loop and residual tip plus tilt DSPs over the frequency range [10 Hz − 75 Hz].The step for the open-loop PSD (tip+tilt) is 5 10 −4 arcsec 2 , to be compared with 3.4 10 −4 for LQG and 10.9 10 −4 for integrator, which clearly amplifies the vibrations.On this example (fairly representative of our on-sky results), the tip+tilt residual variance of the integrator is 107% higher than LQG.We could expect an OMGI controller to give better results, and we have thus evaluated its performance in replay.However, even with a finely tuned OMGI controller (in our simulation, g TT = 0.3 and g HO = 0.7), the residual variance would still remain 43% higher than with LQG.
While many real AO systems exhibits vibration spectra spanning a wide frequency range (see for example [12]), we may also consider ideal 'vibration-free' situations, i.e.where perturbations exhibiting no vibration peak in any 'critical' frequency range.Would high-performance SCAO controllers still improve performance?We have thus conducted tests to give an idea of what could be gained with an LQG controller on such data for which the vibrations within the integrator overshoot zone have been removed.This was achieved by modifying on-sky data sets to obtain filtered data without vibrations in the range 10 − 40 Hz.An example of such spectra, corresponding to the data of August 29 th 2012 at 22:02, is given in Fig. 9.We have evaluated the performance obtained by replaying the loop on these data for the LQG controller -with TT model based on Eqs. ( 30)-( 31) -and for the integrator controller with a fixed gain of 0.7 (the best tuning in replay).The LQG residual variance is around 6% better, thus keeping slightly ahead.Depending on the processed data set, this percentage varies between 5 and 10 %.Note that some improvement could be gained by designing the LQG controller with a system loop delay of 1.6.

LQG vs. integrator: frequency-domain and robustness analysis
This Section presents a comprehensive comparative analysis of LQG and integrator performance from a frequency-domain perspective.We seek to provide answers to the following questions: 1/ Do the disturbance models used for LQG control accurately match the turbulence and vibration modes' PSDs reconstructed from on-sky data?2/ Is the quality of LQG correction for these modes, as measured by their modal rejection transfer functions, significantly better than with an integrator?3/ Is LQG control robust to modeling uncertainties?
These three issues are in fact closely intertwined, because of the fundamental performance/stability trade-off inherent to any control loop: on the one hand, efficient disturbance rejection requires high feedback gains; on the other, (too) high feedback gains render the loop's behavior unstable or at least too sensitive to modeling uncertainties.In order to negotiate a sensible trade-off between these inherently conflicting requirements, the controller gain thus needs to be finely adjusted according to the disturbance's PSD and also to the measurement signal-to-noise ratio.

Model relevance for turbulence and vibration
We discuss here the relevance of the models of high-order (Section 6.1.1)and tip-tilt modes (Section 6.1.2) used for LQG control by comparing their temporal PSDs to those of on-sky data.
6.1.1.Model matching for high orders Figure 10 compares the temporal PSDs of the AR2 models used as priors for turbulence with the experimental PSDs computed on reconstructed Zernike coefficients (Z4 and Z28).We recall that the high order AR2 models are based on the physical description in Eqs. ( 30)- (31), with a cut-off frequency proportional to radial order and wind norm, and an energy related to Kolmogorov spectrum (see Section 4.1).The prior temporal PSDs of these models (with wind norm and r 0 estimated on the data) are plotted in Fig. 10.The experimental PSDs have large dynamics, a high level at low frequency (of course larger for low orders), a sharp decrease after the cut-off frequency and a noise floor at high frequency.We see that the cut-off frequency increases as expected with the radial order of the Zernike polynomial [36].The good agreement between modeled and measured PSDs (even if the cut-off frequencies on left column could be better adjusted) validates the model structure and the parameter estimation.On these curves, the turbulent AR2 models of course do not include the high-frequency noise plateau due to the measurement noise.We now compare the PSDs of on-sky x and y mean slopes with the model identified from this data set as a sum of AR2 processes.Recall that in the turbulent model, the vibration components and the noise plateau due to measurement noise are identified in the frequency domain with the method presented in [17].Figure 11 shows two examples of fits provided by the identification procedure.A good agreement is found both for the smooth turbulence PSDs and for the vibration peaks thanks to the flexibility provided by the adaptation of a damping coefficient for each component.On Fig. 10 (right), a different tuning of the parameter identification algorithm could have led to an even better fit around 15 Hz.However, during all on-sky tests, we retained the same tuning of the identification algorithm.

Rejection transfer functions
In the sequel, we denote by E tur (z) the AO loop's disturbance rejection transfer function (a.k.a. the discrete-time feedback loop's sensitivity function), which is the transfer function between the disturbance φ and the residual phase φ res .To compute E tur (z), one needs to take into account the value of the system loop delay, that is d = 1.6 frame.In practice, a total loop delay 1 < d ≤ 2 means that the control u k is not applied between WFS sampling instants t = (k + 1)T and t = (k + 2)kT , but instead between t = (k + d − 1)T and t = (k + d)T .Thus, the correction phase φ cor k seen by the WFS during the sampling interval so that the measurement equation that needs to be accounted for in the computation of E tur (z) Using this formula, it is immediately checked that the rejection transfer function is given by where G is the controller transfer function (i.e. the open-loop transfer function from residual slopes to control inputs).To characterize the effect of the rejection transfer function on a particular Zernike mode, it suffices to consider the corresponding diagonal term of the sensitivity function.

Rejection transfer functions for high orders modes
In this Section, we compare the rejection transfer functions of integrator and of LQG for the two Zernike modes Z4 and Z28.Because of a fundamental result in control theory known as Bode's integral theorem, which in its discrete-time version applies to all feedback loops with delay d ≥ 1 (see e.g.[42], page 79), increasing rejection for some particular frequency will be repaid in kind by disturbance amplification elsewhere (the so-called 'water-bed effect').As a consequence, it is easily checked that for feedback loops with reasonably good signal-tonoise ratio, the gain of an optimal minimum-variance controller should at low-and-medium frequency be roughly inversely proportional to the square root of the PSD of the perturbation while remaining bounded at high frequencies so as not to amplify measurement noise.This rule of thumb provides a convenient way to gauge, by simple visual inspection, whether or not a particular controller provides optimal or near-optimal rejection [11].
Figure 12 shows the AR2 model temporal PSDs and the rejection transfer functions.As expected, the PSD cut-off frequency increases according to the radial order [36].Note also that LQG's rejection exhibits a low-frequency 'plateau' which closely matches the turbulence's PSD, resulting in more efficient correction.The red curve corresponds to the normalized square root of the PSD of the perturbation's model.

Rejection transfer functions for tip and tilt
Figures 13 and 14 show the rejection transfer functions of the TT modes for LQG control and integrator together with the theoretical PSD of the global AR2 TT model (in red), which is the sum of 11 AR2 processes, one for the TT turbulent component and 10 AR2 for identified vibration peaks.As explained above, the integrator has a high rejection at low frequency, but the rejection degrades rapidly and presents even an overshoot centered at about 20 Hz.On CANARY, as we have shown in Section 6.1.2,there are 2 energetic vibration peaks near this overshoot frequency which significantly degrades the performance.LQG controller with vibration filtering reject those vibrations peaks, leading to a much more efficient correction.be retrieved from E tur i (z) using the identity: Under the assumption that L i (z) is a so-called 'regular' loop transfer function (i.e., without zeros or poles outside the unit circle), the Nyquist stability criterion requires that in order for E tur i (z) to be stable, the Nyquist plot of L i (z) must remain at the right of the critical point -1 for all positive frequencies.The next step is to evaluate the associated gain, phase and delay margins.Thus, denoting as θ p,1 , . . ., θ p,m i the normalized frequencies at which arg(L i (e jθ )) = π, the largest extra gain which can be added without destabilizing the loop is given by the gain margin M G,i defined through Likewise, denoting as θ c,1 , . . ., θ c,n i the normalized frequencies at which |L i (e jθ )| = 1, the largest negative phase shift and the largest extra delay which can be added without destabilizing the loop are given respectively by the phase margin M φ ,i and the delay margin M D,i defined through   show the Nyquist plots of L i (z) for 3 representative Zernike modes (tip, defocus Z3 and Z28) for a system loop delay of 1.6 frames.These Nyquist plots keep clear of the critical point −1, leading to good robustness margins.Table 1 presents the gain, phase and delay margins for 7 different Zernike modes (low, medium and high-order).In order to provide a point of comparison, the table also gives the gain, phase and delay margins for the integrator controller with the standard default setting g = 0.5.This suggests excellent robustness of LQG, especially since in this particular case the Kalman filter and the control gain have been computed for an incorrect total loop delay of two frames.

Conclusion
We have presented in this paper the first SCAO results obtained on-sky with LQG control on all modes.These results, that include vibration mitigation on tip and tilt, have been obtained on the CANARY demonstrator installed at the William Herschel Telescope on the Roques de Los Muchachos in La Palma (Canaries Islands, Spain).A detailed performance analysis is performed and comparison is conducted with the integral action controller.Note that preliminary LQG results obtained in MOAO with CANARY have been presented in [23].LQG control gives better on-sky performance than a single gain integrator controller, and proved to be reliable with perfect loop stability and robust performance.Like most AO systems, the CANARY demonstrator is submitted to several vibrations affecting the tip and tilt modes.Using the frequency domain identification method presented in [17], LQG with vibration filtering leads to significant additional improvement in performance.This validates the tip-tilt control strategy adopted on the European extrasolar planet finder SPHERE at the VLT.Extension to vibration mitigation on other modes is straightforward, as indicated in the formulation presentation.
Replays on on-sky data have allowed to test an Optimized Modal Gain Integrator (OMGI).We show that LQG control still leads to better performance.Besides, even in the absence of vibrations in the overshoot frequency range of the integrator, LQG control keeps a slight advantage.
The designed LQG controller relies on a simple and physical state model adapted to the disturbance's PSD which leads to a better compromise between perturbation rejection and noise propagation in the loop.The selected second order autoregressive perturbation models are in good agreement with CANARY on-sky data.Performance is also analysed through evaluation of rejection transfer functions (RTF) and stability margins.Important attenuations in the RTFs confirm the ability to efficiently damp vibratory components.As for stability, excellent gain, phase and delay margins have been evaluated on the on-sky LQG controllers.
This allows to confidently envision more informative models.We could for instance account for wind direction since estimation of wind profiles in norm and direction is achievable with ever better accuracy [43][44][45][46].This may be of particular interest in the context of WFAO [47].Robustness in stability does not mean robustness in performance, and sensitivity to errors in wind direction estimation thus also needs to be investigated.

Fig. 1 .
Fig. 1.On-sky optimization of the integrator gain during the observation done the 28th of August 2012.Different colors represent three gain tuning sequences at different times of the night.

Fig. 2 .
Fig. 2. Strehl ratios as a function of the Fried parameter for integrator and LQG controllers.Top left: SR in H-band recorded on Aug. 28 th 2012.Top right: SR in H-band recorded on Aug. 29 th 2012.Bottom: SR in K-band recorded on Aug. 29 th 2012.Colors: red stars = integrator control, green triangle = LQG control without vibration mitigation, and magenta square = LQG control with vibration mitigation.

Fig. 3 .
Fig. 3. Example of integrator control relative error budget in H band for a data sample on August 29th 2012 with r 0 = 21 cm.

Fig. 6 .
Fig. 6.Residual phase variance in rad 2 for integrator as a function of global feedback gain g int obtained by replaying the open-loop on-sky data of the 28th August at 23:10.Two different loop delays are considered: d = 1.6 (blue) and d = 2 (magenta).

2 )Fig. 8 .
Fig. 8. Cumulated PSDs for open-loop and residual tip plus tilt disturbance and residuals computed from WFS data recorded on-sky on Aug. 29th 2012.Tip and tilt are evaluated in equivalent on-sky angle, expressed in arcsec, hence PSDs in arcsec 2 .(Open-loop data on bottom plot has been lowered by 0.013 arcsec 2 to increase vertical resolution.)

Fig. 9 .
Fig. 9. On-sky original x-and y-slopes CANARY data (red) and a filtered version (blue) where vibrations have been removed in the range 10-40 Hz.Data set of August 29 th 2012 at 22:02.

Fig. 10 .
Fig. 10.PSDs of high orders Zernike Z4 and Z28.Black = Regularized reconstruction from open-loop data, red = AR2 model.The theoretical cut-off frequency is indicated by a magenta vertical dotted line.The data have been recorded the August 29 th 2012.Left: slopes data at 22h07m32s (V = 10.3m/s).Right: slopes data at 23h07m07s (V = 7.8m/s).Top: Zernike number 4. Bottom: Zernike number 28.The AR2 models have been identified on the sequence of open-loop data at: left: 21h53m26s and right: 22h54m21s.

Fig. 11 .
Fig. 11.Identification of vibrations on Tip (left) and Tilt (right) and fit of the AR2 model for vibration filtering in the perturbation modeling.

Fig. 12 .
Fig. 12. Left: Rejection transfer functions on Z4.Right: Rejection transfer functions on Z28 (observing night 29 August 2012).Blue for integrator control and Magenta for LQG.The red curve corresponds to the normalized square root of the PSD of the perturbation's model.

Fig. 15 .− 1 Fig. 16 .Fig. 17 .
Fig. 15.Nyquist plot of LQG loop transfer function, tip Because NM com D ≈ I, the matrix E tur (z) is approximately diagonal, with the same rejection transfer function for all Zernike modes.