© ESO, 2012

1. Introduction

Ever since the early days of optical interferometry, the motion of fringes due to atmospheric turbulence has been a major limitation to this observation technique. One possible way to avoid the blurring of fringes on the detector is to use short fringe exposures, since short integration times temporally “freeze” atmospheric motion. In combination with the low optical throughput of interferometric systems, however, the sensitivity associated with short integrations is low, even on the largest-telescope systems.

Tests with prototype interferometers (Shao & Staelin 1980; for a further overview, see Monnier 2003) demonstrated the use of the control technique fringe tracking to increase the sensitivity of interferometric observations. Correcting for the atmospherically introduced optical path difference (OPD) between different interferometer arms indeed allows the continuous observation of the fringes. In a fringe tracker, a phase sensor estimates the OPD corresponding to the current fringe exposure. On the basis of this value, a command is calculated by the OPD controller, which is then sent to an OPD-correcting actuator system. The result is that fringes in the scientific beam combiner are stabilized to atmospheric motion, and can be integrated much longer than the typical atmospheric coherence time. In addition, the synchronous tracking of the fringes on a bright source allows the interferometric observation of much fainter sources.

Apart from turbulence, an important limitation of current real-time wavefront correction or tracking is instrumental vibrations. Initial on-sky tests of NAOS, the first adaptive-optics (AO) system on the Very Large Telescope (VLT), attributed a loss in Strehl ratio of up to 15% to mechanical vibrations (Rousset et al. 2003). In the case of the VLT Interferometer (VLTI), longitudinal vibrations on the level of the baselines cause fringes to move onto the detector. Power spectra of sequences of OPD values estimated during a run on the VLTI/PRIMA instrument (Delplancke 2008) indicate the presence of a discrete number of vibration components during fringe observations (Sahlmann et al. 2009). In an effort to reduce the vibration level of the VLTI system, a campaign to identify the contribution of different VLT subsystems has been initiated (Haguenauer et al. 2008; Poupar et al. 2010). The efforts to isolate noisy instruments has led to significant improvements.

Active correction techniques can also be applied to reduce vibration contribution. On the one hand, a vibration-sensing accelerometer system operates on the first three mirrors of the four 8.2-m unit telescopes (UTs) and effectively corrects vibrations in the range 10–30 Hz (Haguenauer et al. 2008). On the other hand, the control algorithm of the OPD controller can be extended with additional vibration blocks that correct for a defined set of vibration components. The Vibration TracKing algorithm (VTK; Di Lieto et al. 2008) is such an example that operates as a narrow-band filter. For the Keck Interferometer, a similar narrow-band control system has been designed (Colavita et al. 2010).

In recent years, there has been an effort in AO to design new control schemes. Of special importance are the techniques based on Kalman filtering (Kalman 1960), which is a technique designed for general data processing. It is used to track and correct observed systems in a statistically optimal way, in that it minimizes the norm of the residual signal (for linear systems with Gaussian noise statistics). Le Roux et al. (2004) designed an AO control scheme based on Kalman filtering. One of the important characteristics of this scheme is that it can be extended (Petit 2006) to coherently treat turbulence and vibrations, rather than add extra vibration control blocks to an unsuitable classical controller. As usual for Kalman filters, it is based on an a priori model for the controlled system, which here is referred to as turbulence and vibration perturbations.

The purpose of this work is the design of Kalman-filtering control schemes for fringe tracking, starting from the work that has been done for AO control. In particular, we consider the case of four-telescope interferometry, and develop control schemes that will be applicable to the fringe-tracking subsystem of GRAVITY (Eisenhauer et al. 2011, 2005), which is a four-telescope beam combiner to be installed on the VLTI. The main purpose of this second-generation instrument is high-precision narrow-angle astrometry and phase-referenced interferometric imaging in the near-infrared K band. The main scientific target of GRAVITY is Sgr A^∗, which is currently the most well-studied and likely supermassive black hole, located at the Galactic center. On-sky operation is expected in 2014.

Conceptually, the idea of using Kalman filtering to correct for atmospheric fringe motion dates back to Reasenberg (1990), but has never been implemented under that form. A simple Kalman-filter based strategy for single-baseline fringe tracking was effectively used at the Palomar Testbed Interferometer (PTI, Colavita et al. 1999), which to our knowledge, is the only example of on-sky fringe tracking using Kalman-filter techniques. Lozi et al. (2010) demonstrated the operation of Kalman filtering for vibration correction on a laboratory prototype for single-baseline space-based nulling interferometry.

The use of Kalman filtering for (ground-based) fringe tracking will have several advantages compared to classical control strategies (e.g. integrator control):

• Kalman filtering allows us to control turbulence and longitudinalvibrations in a coherent and equivalent way: the contribution ofdifferent perturbation components is decomposed;

• all information on disturbances is used via an a priori model. This makes Kalman filtering a statistically optimal control method when considering linear systems with Gaussian white noise statistics;

• the use of Kalman filters may solve issues with existing active vibration control at VLTI (e.g. artifacts in frequency spectra of residual VTK-corrected data; see Haguenauer et al. 2008);

• the predictive nature of the filter allows us to deal with the loss of observables due to measurement noise.

In the case of the latter point, fringe-tracking control schemes need to be robust against flux losses owing to imperfect beam combination. Failures in tip-tilt correction, for instance, can lead to the short-time absence of fringe signal. The prediction capability of the Kalman filter controller is leveraged for this purpose.

The outline of this work is the following. In Sect. 2, we present the AO Kalman-filter control scheme, translated into a form applied to two-telescope fringe tracking. The Kalman filter is based on a parametric model and we present a technique for parameter identification in Sect. 3. Since an aim of this work is to develop fringe-tracker control schemes applicable to GRAVITY, we extend the control in Sect. 2 to four-telescope systems (Sect. 4). After considering some specific improvements in Sect. 5, we introduce the complete control schemes in Sect. 6. Before concluding this work, the results of the first simulations of four-telescope Kalman-filter fringe tracking for GRAVITY are presented (Sect. 7).

2. Basic case: two-telescope Kalman-filter control

The input information provided for fringe-tracking control are the OPD values estimated by the phase sensor. The aim of the OPD controller is to give an appropriate output command to the piston actuators (piezo actuators and delay lines). This section introduces the control law that is appropriate for a two-telescope interferometer. The approach used is based on a formalism developed for AO modeling, hence we can directly apply it here to fringe tracking. For a more complete introduction to this model, we refer to Le Roux et al. (2004) and Petit (2006), which are the papers on which the current section is based. The notation we use is (mainly) that of Meimon et al. (2010).

2.1. Modeling the fringe-tracking control

We have already mentioned before that Kalman filters need a model for the controlled system. We successively introduce the model for the fringe tracking and the evolution of the disturbances (turbulence and vibrations).

2.1.1. Description of the fringe-tracking process

Turbulence and longitudinal vibrations add up to introduce OPD perturbations. The OPD in a two-telescope interferometer can be related to a phase value ϕ as $\mathit{\varphi }\mathrm{\equiv }{\mathit{\varphi }}^{\mathrm{tur}}\mathrm{+}\sum _{\mathit{i}}{\mathit{\varphi }}^{\mathrm{vib}\hspace{0.17em}\mathit{i}}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathit{\lambda }}\mathrm{OPD}\mathit{.}$(1)In reality, the quantity y observed by a fringe tracker is not exactly ϕ, but contains two additional components. The first of these is an additive measurement noise w. Secondly, a fringe tracker in operation will apply correction phases, or commands, u using the piston actuators. We can therefore express the measured value of the fringe tracker as $\mathit{y}\mathrm{=}\mathit{\varphi }\mathrm{-}\mathit{u}\mathrm{+}\mathit{w}\mathit{.}$(2)In what follows, we refer to y as the residual OPD. For simplicity, we ignore conversion factors, e.g. 2π/λ, and assume that all quantities are expressed in the same units (say, microns)¹.

To get into the problem of discretization, we denote by T the elementary integration time of the fringe tracker. The discretized variables are then defined in the following way:

• ${\mathit{\varphi }}_{\mathit{n}}^{\mathit{i}}$ is defined as the average phase of the disturbance component i during the time interval  [(n − 1)T,nT] : ${\mathit{\varphi }}_{\mathit{n}}^{\mathit{i}}\mathrm{\equiv }\frac{\mathrm{1}}{\mathit{T}}{\mathrm{\int }}_{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{T}}^{\mathit{nT}}\mathrm{d}\mathit{t}\hspace{0.17em}{\mathit{\varphi }}^{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}$(3)The sum of all these phases is referred to as \begin{lxirformule}$\varphi_n$\end{lxirformule}.

• The command is a fixed value during one integration (assuming that the response time of the actuators is much shorter than T). We therefore define u_n as the command applied at time step n, i.e. applied during the time interval  [nT,(n + 1)T] .

• The final parameter, the noisy component w_n at time step n, could in principle be defined in a similar way as ${\mathit{\varphi }}_{\mathit{n}}^{\mathit{i}}$. In this model, however, w will be considered as simple zero-mean white Gaussian noise that adds to y. The standard deviation is fixed to σ_w, which makes the values w_n independent of other variables.

The assumption that the data processing after integration amounts to one time step T leads to the following discretized version of Eq. (2): ${\mathit{y}}_{\mathit{n}}\mathrm{=}{\mathit{\varphi }}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{-}{\mathit{u}}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{+}{\mathit{w}}_{\mathit{n}}\mathit{.}$(4)The quantity y_n is thus the measurement that is available at time step n. We include a schematic representation of the fringe measurement model in Fig. 1.

A final note to this model for the fringe-tracking process concerns the notion of a fixed measurement error σ_w. Further on in this work (Sect. 5), we introduce an adaptive weighting according to the (variable) signal-to-noise ratio (S/N) conditions. The current assumption of a fixed σ_w is only for clarity reasons and not for the actual implementation.

2.1.2. Description of the OPD evolution

The formalism of Kalman filtering that is used to estimate the commands, requires an a priori model not only for the measurement process, but also for the evolution of the observed system, that is, of the turbulence and the longitudinal vibrations. In essence, this requires the description of new states of the system in terms of previous realizations.

The OPD spectra obtained at the VLTI typically contain a number of peaks (e.g. Sahlmann et al. 2009), which are associated with vibration modes excited by a variety of processes. In the following, we can assume each mode to be independent. The simplest way is then to describe the dynamics of such a vibration mode i is to consider it as a discrete damped oscillator, excited by a Gaussian white noise v. It can be shown (Petit 2006) that the one-dimensional equation of motion for a vibration mode i of natural frequency ${\mathit{f}}_{\mathrm{0}}^{\mathit{i}}$ and damping coefficient kⁱ can be discretized into the recursive equation ${\mathit{\varphi }}_{\mathit{n}\mathrm{+}\mathrm{1}}^{\mathrm{vib}\hspace{0.17em}\mathit{i}}\mathrm{=}{\mathit{a}}_{\mathrm{1}}^{\mathrm{vib}\hspace{0.17em}\mathit{i}}{\mathit{\varphi }}_{\mathit{n}}^{\mathrm{vib}\hspace{0.17em}\mathit{i}}\mathrm{+}{\mathit{a}}_{\mathrm{2}}^{\mathrm{vib}\hspace{0.17em}\mathit{i}}{\mathit{\varphi }}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{vib}\hspace{0.17em}\mathit{i}}\mathrm{+}{\mathit{v}}_{\mathit{n}}^{\mathrm{vib}\hspace{0.17em}\mathit{i}}\mathit{,}$(5)where $\begin{array}{ccc}{\mathit{a}}_{\mathrm{1}}^{\mathrm{vib}\hspace{0.17em}\mathit{i}}& \mathrm{=}& \mathrm{2}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{\pi }\hspace{0.17em}{\mathit{k}}^{\mathit{i}}{\mathit{f}}_{\mathrm{0}}^{\mathit{i}}\mathit{T}}\hspace{0.17em}\cos (\mathrm{2}\mathit{\pi }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathit{i}}\mathit{T}\sqrt{\mathrm{1}\mathrm{-}{{\mathit{k}}^{\mathit{i}}}^{\mathrm{2}}})\mathit{,}\\ {\mathit{a}}_{\mathrm{2}}^{\mathrm{vib}\hspace{0.17em}\mathit{i}}& \mathrm{=}& \mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{k}}^{\mathit{i}}{\mathit{f}}_{\mathrm{0}}^{\mathit{i}}\mathit{T}}\mathit{.}\end{array}$Equation (5) has a form that is better known as an autoregressive model of order two, AR(2). We note that the conversion from continuous to discrete involves a Taylor approximation.

The modeling of longitudinal vibrations as randomly-excited damped oscillators leads in a natural way to an AR(2) description. One may ask whether the contribution of turbulence can also be described as an AR(p) process. Le Roux et al. (2004) use an AR(1) turbulence model, based on temporal evolution characteristics of the phenomenon and basic parameters such as wind speed and the telescope diameter. A slightly more advanced first-order temporal model is discussed in Looze (2007).

When taking a damping coefficient k > 1, the previous AR(2) vibration model can also be used to describe non-resonating signals, as in the case of turbulence². This property led Meimon et al. (2010) to use an AR(2) model for the turbulence description. We follow their argumentation and write ${\mathit{\varphi }}_{\mathit{n}\mathrm{+}\mathrm{1}}^{\mathrm{tur}}\mathrm{=}{\mathit{a}}_{\mathrm{1}}^{\mathrm{tur}}{\mathit{\varphi }}_{\mathit{n}}^{\mathrm{tur}}\mathrm{+}{\mathit{a}}_{\mathrm{2}}^{\mathrm{tur}}{\mathit{\varphi }}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{tur}}\mathrm{+}{\mathit{v}}_{\mathit{n}}^{\mathrm{tur}}\mathit{.}$(8)

2.2. State-space description of the system control

The state-space representation is a general framework to describe problems involving an evolving system on which measurements are performed (for a general introduction, see Durban 2004). Under the assumption of linearity and time invariance (LTI state-space models), the model has the form of two recursive (constant) matrix equations: one describes the evolving system and the other the observational process on that system.

At least qualitatively, we have good indications that the considered problem can be described as an LTI problem. On the one hand, the states of turbulence and vibrations can be considered as realizations of a time-invariant system (for reasonable timescales). On the other hand, the OPD measurements are clearly a form of observation of that system. Using standard state-space notation, we can cast Eqs. (4), (5), and (8) into the form (Meimon et al. 2010) $\begin{array}{ccc}\begin{array}{c}\\ & & & \end{array}& & \\ \begin{array}{c}\\ & & & \end{array}& & \end{array}$where (for the sake of the example, we assume two vibration components, vib   1 and vib   2) $\begin{array}{ccc}{{x}}_{\mathit{n}}^{\mathrm{\top }}& \mathrm{=}& ({\mathit{\varphi }}_{\mathit{n}}^{\mathrm{tur}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\varphi }}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{tur}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\varphi }}_{\mathit{n}}^{\mathrm{vib}\hspace{0.17em}\mathrm{1}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\varphi }}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{vib}\hspace{0.17em}\mathrm{1}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\varphi }}_{\mathit{n}}^{\mathrm{vib}\hspace{0.17em}\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\varphi }}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{vib}\hspace{0.17em}\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em})\mathit{,}\\ {{v}}_{\mathit{n}}^{\mathrm{\top }}& \mathrm{=}& ({\mathit{v}}_{\mathit{n}}^{\mathrm{tur}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{v}}_{\mathit{n}}^{\mathrm{vib}\hspace{0.17em}\mathrm{1}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{v}}_{\mathit{n}}^{\mathrm{vib}\hspace{0.17em}\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0})\mathit{,}\\ \mathrm{C}& \mathrm{=}& (\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1})\mathit{,}\end{array}$and $\mathrm{A}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\left(\begin{array}{c}\\ {\mathit{a}}_{\mathrm{1}}^{\mathrm{tur}}& {\mathit{a}}_{\mathrm{2}}^{\mathrm{tur}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\mathit{a}}_{\mathrm{1}}^{\mathrm{vib}\hspace{0.17em}\mathrm{1}}& {\mathit{a}}_{\mathrm{2}}^{\mathrm{vib}\hspace{0.17em}\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{a}}_{\mathrm{1}}^{\mathrm{vib}\hspace{0.17em}\mathrm{2}}& {\mathit{a}}_{\mathrm{2}}^{\mathrm{vib}\hspace{0.17em}\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right)\mathit{.}$(14)The vector x is called the state vector of the system. In addition, we denote by Σ_v and ${\mathrm{\Sigma }}_{\mathit{w}}\mathrm{=}{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}$ the covariance matrices of the system noise (v) and measurement noise (w), respectively. We note that for the current example, Σ_v is of the form ${\mathrm{\Sigma }}_{\mathit{v}}\mathrm{\equiv }\mathrm{E}\mathrm{\left\{}{v}{{v}}^{\mathrm{\top }}\mathrm{\right\}}\mathrm{=}\mathrm{diag}({\mathit{\sigma }}_{\mathit{v}}^{\mathrm{tur}}\mathit{,}\mathrm{0}\mathit{,}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{vib}\hspace{0.17em}\mathrm{1}}\mathit{,}\mathrm{0}\mathit{,}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{vib}\hspace{0.17em}\mathrm{2}}\mathit{,}\mathrm{0}{)}^{\mathrm{2}}\mathit{.}$(15)As usual in state-space representation, the state variable x_n is a hidden quantity: only y_n is observed. It is important to highlight the block-diagonal form of the matrix A: every perturbation component (turbulence and vibrations) contributes a single AR(2) block.

2.3. Control: the (asymptotic) Kalman filter

Kalman filtering (Kalman 1960) is a technique designed for general data processing. It is used to track and correct observed systems in a statistically optimal way, in that it minimizes the residual signal (Chui & Chen 2009). Kalman filters start from a linear state-space description of the system and provide estimations of the system’s (hidden) state vector derived from the observations.

We denote by $\hat{{x}}{}_{\mathit{n}\mathrm{|}{\mathit{n}}^{\mathrm{\prime }}}$ the estimation of the state vector at time step n, based on all observations up to time step n′ (i.e., y₀,...,y_n′). Given the state-space model in Eqs. (9) and (10), the Kalman filter equations³ take the form $\begin{array}{ccc}\hat{{x}}\mathit{n}\mathrm{|}\mathit{n}& \mathrm{=}& \hat{{x}}\mathit{n}\mathrm{|}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{+}{\mathrm{G}}_{\mathrm{\infty }}\mathrm{(}{\mathit{y}}_{\mathit{n}}\mathrm{-}\mathrm{C}\hat{{x}}\mathit{n}\mathrm{|}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{+}{\mathit{u}}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{)}\mathit{,}\\ \hat{{x}}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{|}\mathit{n}& \mathrm{=}& \mathrm{A}\hat{{x}}\mathit{n}\mathrm{|}\mathit{n}\mathit{.}\end{array}$The former of these equations updates the information about the disturbance state x_n, taking into account the measurement y_n that has become available at time step n. The latter then predicts the state at the subsequent moment in time.

The matrix G_∞ in Eq. (16) is the (asymptotic) Kalman gain matrix, calculated as ${\mathrm{G}}_{\mathrm{\infty }}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{\infty }}{\mathrm{C}}^{\mathrm{\top }}\mathrm{(}\mathrm{C}{\mathrm{\Sigma }}_{\mathrm{\infty }}{\mathrm{C}}^{\mathrm{\top }}\mathrm{+}{\mathrm{\Sigma }}_{\mathit{w}}{\mathrm{)}}^{-1}\mathit{,}$(18)where Σ_∞ is the solution of the matrix (Riccati) equation ${\mathrm{\Sigma }}_{\mathrm{\infty }}\mathrm{=}\mathrm{A}{\mathrm{\Sigma }}_{\mathrm{\infty }}{\mathrm{A}}^{\mathrm{\top }}\mathrm{-}\mathrm{A}{\mathrm{\Sigma }}_{\mathrm{\infty }}{\mathrm{C}}^{\mathrm{\top }}\mathrm{(}\mathrm{C}{\mathrm{\Sigma }}_{\mathrm{\infty }}{\mathrm{C}}^{\mathrm{\top }}\mathrm{+}{\mathrm{\Sigma }}_{\mathit{w}}{\mathrm{)}}^{-1}\mathrm{C}{\mathrm{\Sigma }}_{\mathrm{\infty }}{\mathrm{A}}^{\mathrm{\top }}\mathrm{+}{\mathrm{\Sigma }}_{\mathit{v}}\mathit{.}$(19)The gain matrix determines the weight by which the new measurement y_n is taken into account for updating our estimate of the state vector x_n. The mathematical form of G_∞ can be derived by requiring that the expectation value $\mathrm{E}\{\mathrm{(}{{x}}_{\mathit{n}}\mathrm{-}\hat{{x}}{}_{\mathit{n}\mathrm{|}\mathit{n}}\mathrm{)}^{\mathrm{2}}\}$ is minimal (e.g. Chui & Chen 2009). We note that the gain matrix contains the noise characteristics (Σ_v,Σ_w), making it optimally adapted to the considered system.

The final piece of information that is needed to design a two-telescope control scheme is an equation to determine commands. An optimal command will minimize the observed residual OPD y_n, which was given by Eq. (10). The command u needs to compensate the x-part (w is the noisy part, which cannot be reduced by u). It is not difficult to see that the most appropriate choice of u is ${\mathit{u}}_{\mathit{n}}\mathrm{=}\mathrm{K}\hspace{0.17em}\hat{{x}}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{|}\mathit{n}\mathit{,}$(20)where (still for the two-vibration example) $\mathrm{K}\mathrm{=}(\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0})\mathit{.}$(21)This result forms the final equation needed for control.

3. Parameter identification

The state-space description that we need to model the fringe-tracking process contains a possibly large number of parameters. Reconsidering the state-space model in Eqs. (9) and (10), we see that the parameters to identify are $\{\mathrm{\left(}{\mathit{a}}_{\mathrm{1}}^{\mathit{i}}\mathit{,}{\mathit{a}}_{\mathrm{2}}^{\mathit{i}}\mathit{,}{\mathit{\sigma }}_{\mathit{v}}^{\mathit{i}}\mathrm{\right)}\hspace{0.17em}\mathrm{|}\hspace{0.17em}\mathrm{\forall }\mathit{i}\mathrm{\in }\mathrm{\{}\mathrm{disturbances}\mathrm{\}}\} \mathrm{and} {\mathit{\sigma }}_{\mathit{w}}\mathit{.}$(22)Once the parameters are obtained, we can “fill” the system matrices, and apply the Kalman filter.

We assume that OPD measurements from the phase sensor are the only available piece of information about the disturbances. In other words, it is assumed that no other a priori information about the disturbance components is available (e.g. from other observation techniques).

3.1. The spectral method

The method we use to calculate model parameters is based on the one proposed by Meimon et al. (2010, hereafter the MPFK method). The MPFK method is a power-spectrum-based routine designed for the parameter identification of the Kalman-filter AO control scheme of Le Roux et al. (2004) and Petit (2006), the scheme that lies at the basis of Sect. 2. We refer to the original MPFK paper for details of the method, and give here a brief summary and some information on our modifications.

The operation of the MPFK method is based on modeling the (estimated) power spectrum of an observed perturbation-phase sequence. In a time interval before identification, a sequence of closed-loop data is obtained and corrected for the applied commands (pseudo-open loop data, POL). In the notation of the state-space model, a sequence of data $\mathrm{\{}{\mathit{y}}_{\mathit{n}}^{\mathrm{POL}}\mathrm{\equiv }{\mathit{y}}_{\mathit{n}}\mathrm{+}{\mathit{u}}_{\mathit{n}\mathrm{-}\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\mathrm{|}\hspace{0.17em}\hspace{0.17em}\mathit{n}\mathrm{=}\mathrm{0}\mathit{,}\mathit{...}\mathit{,N}\mathrm{-}\mathrm{1}\mathrm{\}}$(23)is collected. The periodogram of this sequence contains the imprint of the different perturbation components that contribute to the phase signal (Fig. 2). In the first fitting step, the measurement noise is estimated from the flat tail of the periodogram. By stepwise adding AR(2)-model spectra to the estimated noise spectrum, first for the turbulence and then for the vibration components, the global power spectrum is assembled until a stopping criterion is satisfied. The iterative fitting itself is based on the maximization of the likelihood function associated with each model periodogram. We note that the MPFK method is designed to operate in non-real time, i.e. outside the real-time tracking loop.

In the standard MPFK method, a vibration-peak fitting step makes use of a large parameter grid for the three parameters associated with a single vibration component (${\mathit{a}}_{\mathrm{1}}^{\mathit{i}}\mathit{,}{\mathit{a}}_{\mathrm{2}}^{\mathit{i}}\mathit{,}{\mathit{\sigma }}_{\mathit{v}}^{\mathit{i}}$). The method uses no a priori information about the vibration components, but simply maximizes in each iteration the likelihood criterion corresponding to the given grid. In other words, the fitting of vibration peaks is done in a “blind” way on a large grid. To limit the computation time, we decided to add a simple subroutine that selects peaks from the power spectrum. In every iteration, we flag the points that differ significantly from the estimated power spectrum resulting from the previous iteration. Picking each time the highest peak allows us then to use a much smaller grid around a first guess for the peak parameters, which is made by estimating the maximum, width, and central frequency of the peak. In addition to the gain in computation time, peak picking limits the risk that multiple peaks are being fit with a single large peak, degrading the quality of the fit. We note that our peak-picking method selects peaks based on their height, rather than their energy.

3.2. Practical application

In the MPFK paper, it is shown that N = 2000 is a proper choice for the length of the POL-data acquisition. Since POL data is constructed from closed-loop data, the acquisition of data sequences does not require the loop to be (re)opened. This property is advantageous for the operation of the fringe tracker and opens perspectives to perform updates of the model parameters in closed loop. There is however a fundamental difference between the case of AO, for which the MPFK method was designed originally, and interferometry. In the former case, the initialization of the Kalman filter can be based on a sequence of open-loop data (no POL data is available at initialization). In the case of fringe tracking, however, it is impossible to obtain a sequence of open-loop data without taking the risk that fringes are lost. Unlike the case of AO, where wavefront observables remain well-defined at the level of the wave-front sensor (e.g. microlens images for a Shack-Hartmann wavefront sensor), uncorrected fringes have a typical root mean square (rms) motion that largely surpasses typical detection ranges.

The initialization of the Kalman-filter controller for fringe tracking thus inevitably requires a different strategy than open-loop measurements. A straightforward solution is to close the loop temporally using a classical control algorithm. POL data obtained in this way can then be fitted, and consequently we switch from the classical to the Kalman-filter controller. We note that this technique assumes that we take into account the difference in measurement error: the measurement noise σ_w for a classical controller, extracted from the POL-data fit, will be higher than for the Kalman filter. The measurement error for the latter can be estimated based on Eq. (59), which is introduced in Sect. 7.

3.3. Other methods

The spectral identification method is based on a maximum-likelihood method to model periodograms of measurement sequences. Alternative methods exist to identify parameters in state-space models (see also Meimon et al. 2010).

The advantage of the spectral method is that it is rather transparent, and gives satisfying results (see Sect. 7). The reason is that the imprint of the different perturbation components can be largely decoupled in frequency space, which makes it possible to separate the different disturbances. Future tests on real data will have to justify the use of the method for fringe-tracking purpose, and allow us to optimize the parameter-grid construction. A disadvantage of the method, however, is that it is difficult in general to make an iterative fitting routine robust against outliers and non-standard situations. A good alternative to the spectral-fit method could be a method of a purely algebraic nature, e.g. similar to the one proposed by Hinnen et al. (2005). We are currently checking the possibility of using methods of this nature for parameter identification, as an alternative to the spectral method introduced here. Tests of identification methods based on the reconstruction of the measured time sequences (Shumway & Stoffer 1982; Ziskand & Hertz 1993) were not successful in the current case, basically owing to the large number of components that need to be characterized.

4. Four-telescope fringe-tracking control

In the two previous sections, we have presented a general and complete control scheme for two-telescope fringe tracking. It is generally known that two-telescope interferometers provide a rather limited amount of information per observation (one point in the uv-plane per wavelength bin, no direct phase information). Optical interferometry needs to be extended to include more apertures, especially when considering image reconstruction.

The purpose of the current section is to find how two-telescope control can be extended to more telescopes. We chose to consider the case of four-telescope fringe tracking, in the framework of the second-generation near-infrared beam combiner GRAVITY for VLTI. GRAVITY will combine the beams of the four UTs (or four 1.8-m Auxiliary Telescopes, ATs) in a completely redundant way, i.e. on six baselines. Per elementary fringe exposure, the instrument therefore measures six OPDs, rather than one OPD for a two-telescope beam combiner. The OPDs are corrected by a system of four actuators, i.e. one per telescope. We now consider how the six observables can be transformed into proper commands for the actuator system. We note that this problem is absent in the case of a two-telescope beam combiner: one OPD measurement corresponds to one command. Our four-telescope approach can easily be generalized to n-telescope fringe tracking.

4.1. General modal-based control

Piston disturbances acting on a four-telescope interferometer can be considered as vectors in a four-dimensional space. In this space, an arbitrary disturbance P can be expressed as ${P}\mathrm{=}({\mathit{P}}^{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}^{\mathrm{1}}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}^{\mathrm{2}}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}^{\mathrm{3}}{)}^{\mathrm{\top }}\mathit{.}$(24)The actual piston disturbances P are not directly observable: only piston differences or OPDs are measured. In a similar way to above, we can define OPD vectors as ${O}{P}{D}\mathrm{=}({\mathrm{OPD}}^{\mathrm{01}}\hspace{0.17em}\hspace{0.17em}{\mathrm{OPD}}^{\mathrm{02}}\hspace{0.17em}\hspace{0.17em}{\mathrm{OPD}}^{\mathrm{03}}\hspace{0.17em}\hspace{0.17em}{\mathrm{OPD}}^{\mathrm{12}}\hspace{0.17em}\hspace{0.17em}{\mathrm{OPD}}^{\mathrm{13}}\hspace{0.17em}\hspace{0.17em}{\mathrm{OPD}}^{\mathrm{23}}{)}^{\mathrm{\top }}\mathit{,}$(25)where OPD^ij ≡ P^j − Pⁱ. For simplicity, we have assumed that pistons and OPDs are expressed in the same units, say microns. The conversion P → OPD is then given by the matrix equation ${O}{P}{D}\mathrm{=}\mathrm{M}\hspace{0.17em}{P}\mathit{,}$(26)where ${\mathrm{M}}^{\mathrm{\top }}\mathrm{=}\left(\begin{array}{c}\\ -1& -1& -1& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& -1& -1& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}& -1\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{1}\end{array}\right)\mathit{.}$(27)It is easy to show that the system in Eq. (26) is not invertible. Mathematically, the matrix M is of rank three, rather than four. Physically, this corresponds to it being impossible to recover the global piston P^tot ≡ P⁰ + P¹ + P² + P³ from OPD observations (in a similar way to the piston not being observable with a single-dish telescope). A four-telescope fringe tracker therefore only needs to correct for three net observables.

Rather than using a description that is based on pistons, we perform a matrix transformation, based on a singular value decomposition of the matrix M, to construct a set of variables that directly isolates the invisible global piston. We define the mode four-vector Q (Qⁱ for i = 0,1,2,3) as ${Q}\mathrm{\equiv }{\mathrm{V}}^{\mathrm{\top }}\hspace{0.17em}{P}\mathit{,}$(28)where $\mathrm{V}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\begin{array}{c}\\ -1& -1& -1& \mathrm{1}\\ -1& \mathrm{1}& \mathrm{1}& \mathrm{1}\\ \mathrm{1}& -1& \mathrm{1}& \mathrm{1}\\ \mathrm{1}& \mathrm{1}& -1& \mathrm{1}\end{array}\right)\mathit{.}$(29)This orthogonal transformation isolates (a multiple of) P^tot into Q³. The Q-modes form an equivalent to the Zernike modes in AO, and their relation to the pistons is graphically represented in Fig. 3. The conversion between vectors OPD and Q is then given by ${Q}\mathrm{=}\mathrm{H}\hspace{0.17em}\hspace{0.17em}{O}{P}{D}\mathit{,}$(30)where $\mathrm{H}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\left(\begin{array}{c}\\ \mathrm{0}& \mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{1}& -1& \mathrm{0}& \mathrm{1}\\ \mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& -1& -1\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{.}$(31)We note that the matrix H automatically puts Q³ to 0, indicating again that we have no information about P^tot from vectors OPD. The three components Q⁰, Q¹, and Q² now form a proper description of the observable system with three degrees of freedom. A general control scheme based on modes consists then of three steps:

• 1.

  we calculate the residual mode vector y_Q from the measured residual OPD 6-vector y using the transformation in Eq. (30) ${{y}}_{\mathit{Q}}\mathrm{=}\mathrm{H}\hspace{0.17em}{y}\mathrm{;}$(32)

• 2.

  using three controllers (e.g. integrator controllers), we calculate a proper Q-command u_Q to correct the Q-modes. We fix ${\mathit{u}}_{\mathit{Q}}^{\mathrm{3}}\mathrm{\equiv }\mathrm{0}$, that is, the global piston is fixed (by convention);

• 3.

  transform u_Q into the four-component piston command u using the inverse transformation of Eq. (28): u = V   u_Q. These four commands are sent to the actuators for piston correction.

4.2. Modal Kalman-filter control

Above, we introduced a properly defined control scheme based on a modal approach. We now apply a similar philosophy to a Kalman-filter based control model, and find the first appropriate control law for four-telescope fringe tracking.

4.2.1. Modal state-space model

Starting from the state-space description in Sect. 2.3, we derived (see Appendix A) a full modal state-space model for a four-telescope fringe-tracking system, which can be easily generalized to n telescopes. Individual OPD measurements are assumed to be uncoupled in this model. Additionally, all piston perturbations are considered as being independent. We note that this condition is imposed on two scales. First, on the scale of one telescope, we ignore the coupling between different vibration components. Secondly, on the scale of all telescopes, we assume neither that coupling exists between different telescope vibrations, nor that the turbulence perturbation on the pistons is coupled⁴.

The result of this analysis is the modal state-space model (compare with Eqs. (9) and (10)) given by $\begin{array}{ccc}{{x}}_{\mathit{Q,n}\mathrm{+}\mathrm{1}}& \mathrm{=}& \begin{array}{}\mathrm{􏽥}\\ \mathrm{A}\end{array}\hspace{0.17em}{{x}}_{\mathit{Q,n}}\mathrm{+}{{v}}_{\mathit{Q,n}}\mathit{,}\\ {{y}}_{\mathit{Q,n}}& \mathrm{=}& \begin{array}{}\mathrm{􏽥}\\ \mathrm{C}\end{array}\hspace{0.17em}{{x}}_{\mathit{Q,n}}\mathrm{-}{{u}}_{\mathit{Q,n}\mathrm{-}\mathrm{2}}\mathrm{+}{{w}}_{\mathit{Q,n}}\mathit{.}\end{array}$The vector quantities have the following definitions:

• x_Q:

  4-block modal state vector;

• y_Q:

  modal measurement 4-vector (${\mathit{y}}_{\mathit{Q}}^{\mathrm{3}}\mathrm{\equiv }\mathrm{0}$);

• u_Q:

  modal command 4-vector (${\mathit{u}}_{\mathit{Q}}^{\mathrm{3}}\mathrm{\equiv }\mathrm{0}$);

• v_Q:

  4-block modal system noise vector, of covariance matrix Σ_Q,v;

• w_Q:

  measurement noise 4-vector, of covariance; Σ_Q,w note that ${\mathit{w}}_{\mathit{Q}}^{\mathrm{3}}\mathrm{\equiv }\mathrm{0}$.

On the other hand, the matrices $\begin{array}{}{\mathrm{􏽥}}^{}\\ \mathrm{A}\end{array}$ and $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathrm{C}\end{array}$ are $\begin{array}{}\mathrm{􏽥}\\ \mathrm{A}\end{array}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{A}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{A}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{A}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{A}\end{array}\right) \mathrm{and} \begin{array}{}\mathrm{􏽥}\\ \mathrm{C}\end{array}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{C}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{C}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{C}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{,}$(35)where the blocks A and C have the same structure as the corresponding system matrices in the two-telescope state-space model. All scalar quantities (y_n, u_n, w_n) have now become four-vectors, and the vector/matrix quantities (x_n, A, etc.) have obtained a four-block structure (x_Q,n, $\begin{array}{}{\mathrm{􏽥}}^{}\\ \mathrm{A}\end{array}$, etc.). This is all because there are four modes. There are several interesting properties about the above result, which we discuss now. We note that we first neglect the noise terms.

First, all four terms in Eq. (34) have a last component that is immediately defined to be 0 (see also matrix $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathrm{C}\end{array}$). In other words, the global piston is automatically defined as not being observable. Thus, although it is still a real evolving quantity described by the last block component of the vectors in Eq. (33), we can neglect it for the analysis.

Secondly, when looking at the structure of the state-space model and of Eq. (35), we see that all modes have a decoupled evolution (still neglecting the noise). What is more is that the sub-matrices of $\begin{array}{}{\mathrm{􏽥}}^{}\\ \mathrm{A}\end{array}$ are the same, and the same observation holds for $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathrm{C}\end{array}$. Neglecting thus the noise, we see that the evolution is described by four identical copies of a two-telescope state-space model (with the modification that for the last copy the observation equation is 0). Modes therefore essentially behave in a similar way to the quantity ϕ of the previous sections⁵.

Considering the noise terms.

A few things change when we also take into account the noise terms. Since modes are calculated as combinations of pistons, Eq. (28), or combinations of OPDs, Eq. (32), the covariance matrices are in general no longer diagonal (given our assumption that piston perturbations as well as OPD measurements are independent).

Since we have individual OPD measurements at our disposal, it is possible to estimate the measurement errors ${\mathit{\sigma }}_{\mathit{w}}^{\mathit{ij}}$ associated with the measurements of OPD^ij. This allows us to calculate the mode measurement error covariance matrix Σ_Q,w as ${\mathrm{\Sigma }}_{\mathit{Q,w}}\mathrm{=}\mathrm{H}\hspace{0.17em}{\mathrm{\Sigma }}_{\mathit{w}}\hspace{0.17em}{\mathrm{H}}^{\mathrm{\top }}\mathit{,}$(36)where ${\mathrm{\Sigma }}_{\mathit{w}}\mathrm{=}\mathrm{diag}({\mathit{\sigma }}_{\mathit{w}}^{\mathrm{01}}\mathit{,}{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{02}}\mathit{,}{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{03}}\mathit{,}{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{12}}\mathit{,}{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{13}}\mathit{,}{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{23}}{)}^{\mathrm{2}}$(37)(we refer to Appendix A). The evaluation of the covariance matrix Σ_Q,v, however, is a more fundamental problem. The matrix Σ_Q,v depends on quantities that require direct piston information, which is inaccessible (owing to the unobservability of the absolute pistons). We neglect the off-diagonal entries; under this approximation, it can be shown that Σ_Q,v has a diagonal built of four equal diagonal parts (e.g. like the structure of $\begin{array}{}{\mathrm{􏽥}}^{}\\ \mathrm{A}\end{array}$). In this model, the evolution of the four modes by Eq. (33) is thus genuinely independent and equal.

Convention.

Up to now, we have carried the global piston mode as a superfluous burden in our matrix equations. Since we assume that the modes are fully decoupled, we can drop all matrix and vector blocks corresponding to Q³. Under this convention, the fundamental matrices V and H become $\mathrm{V}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\begin{array}{c}\\ -1& -1& -1\\ -1& \mathrm{1}& \mathrm{1}\\ \mathrm{1}& -1& \mathrm{1}\\ \mathrm{1}& \mathrm{1}& -1\end{array}\right) \mathrm{and} \mathrm{H}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\left(\begin{array}{c}\\ \mathrm{0}& \mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{1}& -1& \mathrm{0}& \mathrm{1}\\ \mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& -1& -1\end{array}\right)\mathit{.}$(38)It is not difficult to generalize this operation to the other quantities. For example, the measurement-error covariance matrix Σ_Q,w loses its last column and line and becomes a 3 × 3 matrix⁶.

4.2.2. Modal Kalman filter

The non-diagonal nature of the covariance matrix Σ_Q,w shows that the modes are not independently measured. It follows that the control scheme has to involve one global Kalman filter for the three modes, rather than one per mode. The full Kalman-filter equations then take the form $\begin{array}{ccc}\hat{{x}}\mathit{Q,n}\mathrm{|}\mathit{n}& \mathrm{=}& \hat{{x}}\mathit{Q,n}\mathrm{|}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{+}{\mathrm{G}}_{\mathrm{\infty }}\hspace{0.17em}\mathrm{(}{{y}}_{\mathit{Q,n}}\mathrm{-}\begin{array}{}\mathrm{􏽥}\\ \mathrm{C}\end{array}\hspace{0.17em}\hat{{x}}\mathit{Q,n}\mathrm{|}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{+}{{u}}_{\mathit{Q,n}\mathrm{-}\mathrm{2}}\mathrm{)}\mathit{,}\\ \hat{{x}}\mathit{Q,n}\mathrm{+}\mathrm{1}\mathrm{|}\mathit{n}& \mathrm{=}& \begin{array}{}\mathrm{􏽥}\\ \mathrm{A}\end{array}\hspace{0.17em}\hat{{x}}\mathit{Q,n}\mathrm{|}\mathit{n}\mathit{.}\end{array}$In this equation, the Kalman gain G_∞ is calculated as in Eqs. (18) and (19), with the changes $\mathrm{A}\mathit{,}\mathrm{C}\mathit{,}{\mathrm{\Sigma }}_{\mathit{v}}\mathit{,}{\mathrm{\Sigma }}_{\mathit{w}}\mathrm{\to }\begin{array}{}{}^{\mathrm{􏽥}}\\ \mathrm{A}\end{array}\mathit{,}\begin{array}{}{}_{\mathrm{􏽥}}\\ \mathrm{C}\end{array}\mathit{,}{\mathrm{\Sigma }}_{\mathit{Q,v}}\mathit{,}{\mathrm{\Sigma }}_{\mathit{Q,w}}$. Finally, the optimal mode command is calculated as ${{u}}_{\mathit{Q,n}}\mathrm{=}\begin{array}{}\mathrm{􏽥}\\ \mathrm{K}\end{array}\hspace{0.17em}\hat{{x}}\mathit{Q,n}\mathrm{+}\mathrm{1}\mathrm{|}\mathit{n}\mathit{,}$(41)where the matrix $\begin{array}{}{\mathrm{􏽥}}^{}\\ \mathrm{K}\end{array}$ is defined in the same way as $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathrm{C}\end{array}$ in Eq. (35), with blocks K defined as in Sect. 2.3.

4.2.3. Practical considerations

The modal Kalman-filter control is an interesting control scheme: it is completely symmetric, its conversion matrices are simple, its input and output have the same dimensions (three), and (under the above approximations) each mode has exactly the same state-space model. An additional advantage of the modal scheme is presented in Sect. 4.4, where we consider the identification of the parameters in the state-space model described in Eqs. (33) and (34).

The symmetry involved in the definition of the modes implies that the control will work most effectively when all other involved quantities are symmetric. It is easy to see, however, that this is not the case. The design and implementation of the GRAVITY system are largely symmetric, whereas actual observation conditions (e.g. wind shake, vibrations) and baseline configurations are not. Whenever a single component in the system fails, for example one telescope or one baseline, at least two of the three observables are spoiled in the current control scheme (each OPD contributes to two modes). A priori, the loss of observables might be expected to have a negative impact on the control. We therefore propose a second four-telescope control scheme, which is based on the tracking of the individual OPDs.

4.3. OPD Kalman-filter control

An alternative to a global Kalman filter that tracks all modes is to use six Kalman filters to control each individual OPD. In this (redundant) scheme, each telescope pair is considered as a single two-telescope interferometer.

One might wonder why it would be more interesting to track twice the number of variables (three modes versus six OPDs). The main interest of the modal control is that it was defined in a perfectly symmetric way to transform the six observables and four pistons into three quantities that are tracked. Here, the conversion to and from modes is done using fixed matrices. In the OPD control, however, the redundancy of the conversion to and from OPDs allows us to choose different ways of defining commands. In particular, one can consider non-symmetric combinations of observables or – even better – adaptively weight observables with respect to the physical conditions of observations. The latter possibility is our main point of interest: whenever observing conditions (low flux, low-signal baselines, etc.) make certain OPDs unfavorable, we can lower their weight in the command calculation. Additionally, it even allows us to let the interferometer work with only three or two telescopes.

We postpone a deeper discussion of the above problem to Sect. 5, where we also modify the modal scheme to a more robust version. For completeness, we already give an appropriate method for command-vector calculation. By x_OPD, we denote the six-block OPD state vector, i.e. the state vector corresponding to each OPD organized into a column vector. The four-component piston command u_n is then calculated as $\begin{array}{ccc}{{u}}_{\mathit{n}}\mathrm{=}{\mathrm{M}}_{\mathit{W}}^{†}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ \mathrm{K}\end{array}\hspace{0.17em}\hat{{x}}\mathrm{OPD}\mathit{,n}\mathrm{+}\mathrm{1}\mathrm{|}\mathit{n}\mathit{,}& & \end{array}$(42)where $\begin{array}{}{\mathrm{􏽥}}^{}\\ \mathrm{K}\end{array}$ is the six-block version of the previously introduced matrix K, in order to operate on the six substate vectors in $\hat{{x}}{}_{\mathrm{OPD}}$. The matrix , on the other hand, is a weighted generalized inverse to the matrix M (Eq. (27)), calculated as $\begin{array}{ccc}{\mathrm{M}}_{\mathit{W}}^{†}\mathrm{=}\mathrm{(}{\mathrm{M}}^{\mathrm{\top }}\mathrm{W}\hspace{0.17em}\mathrm{M}{\mathrm{)}}^{†}\hspace{0.17em}{\mathrm{M}}^{\mathrm{\top }}\mathrm{W}\mathit{.}& & \end{array}$(43)In this equation, the 6 × 6 weighting matrix W, defined further on in Eq. (47), distributes the weight among the different OPD substates for the command calculation.

4.4. Four-telescope parameter identification

The MPFK method (including the peak-pick extension) can be directly applied to two-telescope fringe tracking, for the simple reason that it is exactly designed for the corresponding state-space model. In the four-telescope case, the sequences of ${\mathit{y}}_{\mathit{n}}^{\mathrm{POL}}$ (Eq. (23)) are replaced by equivalent sequences of six-dimensional vectors ${{y}}_{\mathit{n}}^{\mathrm{POL}}$, i.e. the POL OPD measurement six-vectors.

In the modal state-space model, the three observed modes are decomposed into the elementary disturbances. Therefore, the most logical option is to apply the spectral fit to the time sequences $\{{y}\begin{array}{}\mathrm{POL}\\ \mathit{Q,n}\end{array}\mathrm{\equiv }\mathrm{H}{{y}}_{\mathit{n}}^{\mathrm{POL}}\hspace{0.17em}\mathrm{|}\hspace{0.17em}\mathit{n}\mathrm{=}\mathrm{0}\mathit{,}\mathit{...}\mathit{,N}\mathrm{-}\mathrm{1}\}\mathit{,}$(44)where ${{y}}_{\mathit{Q,n}}^{\mathrm{POL}}$ are the three-component POL mode measurements (we use Eq. (32)). Under the convention of neglecting the off-diagonal entries of the covariance matrix Σ_Q,v, we have shown in Sect. 4.2.1 that the state-space model is the same for the three modes. The result of this approximation is not only that we can apply the spectral fit to each of the three non-trivial modes, but also that we can average the spectra of the three modes and perform one global fit. Indeed, assuming three independent modes, the three periodograms are independent realizations of the same state-space model. Averaging periodograms for identification translates into a modification of the periodogram statistics, and thus a new likelihood function for the identification method (see Appendix B for some details). The important advantage is that the periodogram is a more accurate estimation of the power spectrum, and that even small perturbation components become significantly discernible from statistical noise.

The OPD Kalman-filter control scheme, on the other hand (Sect. 4.3), is based on the redundant tracking of all OPDs by six two-telescope Kalman filters. In this case, it thus suffices to apply the spectral fit to each of the six OPD measurement sequences. The disadvantage is, of course, that we have to apply the identification method six times, instead of once in the modal scheme. Additionally, since we have only one sequence per OPD, we cannot apply the advantageous averaging of periodograms.

5. Fine-tuning the fringe-tracking control

Our aim to find appropriate Kalman-filter control schemes led us to consider two options: tracking three observable modes in a three-component Kalman filter and tracking six OPDs in six individual Kalman filters. We now consider some specific cases encountered in real observations, and present a way of dealing with them in our control schemes.

5.1. Low S/N baselines

Different processes can lead to lower-quality OPD measurements on a given baseline. Examples are non-central fringe measurements in the fringe envelope and variable splitting ratios of the beams. Whenever the fringe amplitude on a certain baseline is low with respect to the noise level, the control strategy should be designed to maintain the fringe tracking at as optimal a level as possible.

An advantage of the interferometer concept considered above is that its redundant architecture enables it to take into account low fringe-amplitude baselines, in which different OPD combinations can compensate for other OPD measurements (e.g. OPD⁰¹ = OPD⁰² − OPD¹²). The key step will be to take a weighted recombination of the six baselines, in order to mutually improve the raw measurements.

In essence, we can write the process of weighting the residual OPD values as $\begin{array}{ccc}{y}\mathit{W}\mathrm{=}{\mathrm{I}}_{\mathit{W}}\hspace{0.17em}{y}\mathit{,} \mathrm{where} {\mathrm{I}}_{\mathit{W}}\mathrm{=}\mathrm{M}\hspace{0.17em}{\mathrm{M}}_{\mathit{W}}^{†}\mathit{,}& & \end{array}$(45)with an associated measurement-error covariance matrix ${\mathrm{\Sigma }}_{\mathit{w,W}}\mathrm{=}{\mathrm{I}}_{\mathit{W}}\hspace{0.17em}{\mathrm{\Sigma }}_{\mathit{w}}\hspace{0.17em}{\mathrm{I}}_{\mathit{W}}^{\mathrm{\top }}\mathit{.}$(46)The 6 × 6 matrix I_W combines in a weighted way the OPD measurements, which allows us to mutually improve each OPD value. The weights W^ij attributed to each OPD^ij are organized in a diagonal weighting matrix W, which is used to calculate the previously introduced matrix (Eq. (43)). In accordance with the purpose of weighting, a proper definition of W is $\mathrm{W}\mathrm{=}{\mathrm{\Sigma }}_{\mathit{w}}^{-1}\mathit{,}$(47)which is the inverse of the (diagonal) measurement-error covariance matrix. In other words, we take ${\mathit{W}}^{\mathit{ij}}\mathrm{=}\mathrm{(}{\mathit{\sigma }}_{\mathit{w}}^{\mathit{ij}}{\mathrm{)}}^{-2}$. We note that either global weights W, i.e. fixed weights during operation, or instantaneous weights W_n, i.e. defined for each OPD measurement, can be considered. The latter is possible owing to an instantaneous measurement error that is provided by the phase sensor (we refer to Eq. (59) in Sect. 7).

In the modal Kalman-filter control scheme, we now simply apply the standard OPD-to-mode transformation in Eq. (32) to the weighted modes y_W, and accordingly use a weighted measurement covariance matrix ${\mathrm{\Sigma }}_{\mathit{Q,w,W}}\mathrm{=}\mathrm{H}\hspace{0.17em}{\mathrm{\Sigma }}_{\mathit{w,W}}\hspace{0.17em}{\mathrm{H}}^{\mathrm{\top }}\mathit{.}$(48)For the OPD scheme, the six Kalman filters are applied to the six weighted OPDs. We note that in this scheme, the weighting is assumed not to introduce coupling into the OPDs. In other words, we ignore off-diagonal terms in the covariance matrix Σ_w,W in the OPD scheme. This allows us to keep using the 6 independent Kalman filters, and will prove useful for the adaptive-gain definition below.

A final note on this OPD weighting is that it is only useful for non-resolved fringe-tracking targets. When the target is (significantly) resolved, the intrinsic visibility phases may be non-zero and hence phase closures cannot be assumed to be 0 (e.g. OPD⁰¹ + OPD¹² − OPD⁰² ≠ 0, in general). In any case, fringe tracking of significantly resolved objects is disadvantageous owing to the lower fringe contrasts.

5.2. Flux dropouts

OPD weighting resolves the issue of low S/N’s on individual baselines. A second problem occurs when the flux acquired at a telescope is low, i.e. the so-called flux dropouts. The notion of flux dropouts contains the different processes by which the number of photons coming from a certain telescope drops. These include:

• atmosphere-related effects, where scintillation and otherturbulence effects can lead to a failure of the AO system. Theresulting Strehl ratio can be insufficient for proper fringeformation and detection. In addition, cirrus may lower the fluxarriving at the focal plane of a telescope;

• errors in the fiber injection, where mechanical problems (e.g. vibrations) and other propagating tip-tilt errors in the optical train can lead to problems in the injection of the beams into the fibers, coupled to the beam combiner.

As far as the second point is concerned, GRAVITY will have a fiber coupler subsystem to perform an internal tip-tilt correction (Pfuhl et al. 2010). Owing to the residual errors, however, a specially developed fringe-tracking strategy is still required.

It should be clear that the problem of flux dropouts is very different from the problem introduced in the previous section. When the flux drops, it is impossible to compensate for the low signals on the baselines associated with that telescope: all these baselines are affected in a similar way. For non-predictive control schemes, no other solution exists apart from waiting until the lost signal re-emerges, while keeping the command for that telescope fixed. The advantage of using a predictive method, such as a Kalman filter, is that it allows us to continue operating based on the state-space model.

In Appendix C, we argue for the necessity to modify the gain to take into account flux problems at the level of a telescope. In essence, we modify the constant Kalman gain and allow it to decrease when observing conditions on associated telescopes degrade. The result is that the estimation of the next state is based more on the model, i.e. the observed (noisy) residual OPD is not taken into account as effectively.

To introduce the instantaneous-gain definitions, we first define the instantaneous versions Σ_Q,w,W,n and Σ_w,W,n of the previously introduced noise-covariance matrices (Eqs. (46) and (48)). The instantaneous covariances are recalculated at each time step using instantaneous weighting matrices W_n (rather than a constant W). The discussion in Appendix C then leads us to the instantaneous-gain definitions $\begin{array}{ccc}{\mathrm{G}}_{\mathrm{\infty }\mathit{,n}}\mathrm{=}\frac{\mathrm{tr}\{{\mathrm{\Sigma }}_{\mathit{Q,w,W}}\}}{\mathrm{tr}\{{\mathrm{\Sigma }}_{\mathit{Q,w,W,n}}\}}\hspace{0.17em}{\mathrm{G}}_{\mathrm{\infty }} \mathrm{\left[}\mathrm{modalscheme}\mathrm{\right]}& & \end{array}$(49)(“tr” denotes taking the trace) and ${\mathrm{G}}_{\mathrm{\infty }\mathit{,n}}^{\mathit{ij}}\mathrm{=}\frac{({\mathrm{\Sigma }}_{\mathit{w,W}}{)}^{\mathit{ij}}}{({\mathrm{\Sigma }}_{\mathit{w,W,n}}{)}^{\mathit{ij}}}{\mathrm{G}}_{\mathrm{\infty }}^{\mathit{ij}} \mathrm{\left[}\mathrm{OPDscheme}\mathrm{\right]}$(50)(the index ij to the parentheses indicates that we take the matrix diagonal element corresponding to the considered OPD^ij; for instance, for ij = 01 and ij = 12, we take the zeroth and third diagonal element, respectively).

The advantage of the OPD scheme is now clear. Here we consider the extreme case when we lose all information provided by one telescope. In the modal scheme, the loss of information immediately means that we lose all mode observables, which implies that we cannot take into account the observational information at that step (the instantaneous gain matrix becomes zero). On the other hand, losing a telescope in the OPD scheme means we still have three OPDs that are properly measured, and the gain for those Kalman filters can stay fixed. This implies that the command is still based on half the number of usual observables (versus none for the modal scheme). The cost we have to pay is that we need to track six variables instead of three.

5.3. Piston isolation

For some specific reasons, it might be interesting to completely decouple one or more telescopes from the command estimate. This can occur when a telescope (or its AO system) temporally becomes very unreliable during an observing run. The flexibility of the OPD scheme allows us to track blindly on one/two telescope(s) without affecting the tracking of the others, after which we can switch back to the usual scheme.

To temporally decouple telescope i, we attribute a weight zero to all associated baselines to telescope i. Doing so, however, results in putting the command on telescope i to 0. This can be disadvantageous to the stability of the system, as it generally involves a large jump in all actuator positions (the average actuator position is always 0, hence all actuator positions change when one command is put to 0). The basic idea is then to add an extra component to the command vector, $\begin{array}{ccc}{{u}}_{\mathit{n}}\mathrm{=}{\mathrm{M}}_{\mathit{W}}^{†}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ \mathrm{K}\end{array}\hspace{0.17em}\hat{{x}}\mathrm{OPD}\mathit{,n}\mathrm{+}\mathrm{1}\mathrm{|}\mathit{n}\mathrm{+}{\mathrm{L}}^{\mathit{i}}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ \mathrm{K}\end{array}\hspace{0.17em}\hat{{x}}\mathrm{OPD}\mathit{,n}\mathrm{+}\mathrm{1}\mathrm{|}\mathit{n}\mathit{,}& & \end{array}$(51)where the components of $\hat{{x}}{}_{\mathrm{OPD}\mathit{,n}\mathrm{+}\mathrm{1}\mathrm{|}\mathit{n}}$ that correspond to i are blindly estimated (gain 0). The newly-introduced matrix Lⁱ, on the other hand, has a specific structure that adds a blindly estimated command to i, and subtracts an equal value from all other commands such that the sum of the commands is zero. For instance, for i = 0 (i.e. to isolate P⁰) and i =  { 1,2 }  (i.e. to isolate both P¹ and P²), we have ${\mathrm{L}}^{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\left(\begin{array}{c}\\ \mathrm{3}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -1& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -1& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -1& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right) \mathrm{and} {\mathrm{L}}^{\mathrm{1}\mathit{,}\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\begin{array}{c}\\ \mathrm{0}& -1& -1& \mathrm{0}\\ \mathrm{0}& \mathrm{2}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{2}& \mathrm{0}\\ \mathrm{0}& -1& -1& \mathrm{0}\end{array}\right)\mathit{,}$(52)respectively. As we subtract an equal command from each non-decoupled actuator, these commands on the decoupled telescope(s) do not affect the non-decoupled telescopes.

6. Full Kalman-filter operation

We have introduced the basic theory, developed several control schemes, and considered strategies that allow us to minimize the problems associated with imperfections of the system/control schemes. It is a good point now to define the full Kalman-filter operation by comparing the two schemes. The operation is split into an identification phase and a (real-time) tracking phase, the former being needed to find the parameters.

Identification phase

• 1.

  Pseudo-open loop (POL) data acquisition: a set of data pointsand estimated measurement errors is obtained in closed loop andcorrected for the applied commands to get a POL sequence.

• 2.

  Weighting: we estimate a global measurement noise for each POL OPD sequence by taking the median of the measurement errors. We apply the weighting in Eq. (45) to get the weighted OPD sequence. In addition, for the modes, apply the OPD-to-mode transformation in Eq. (32) to calculate the weighted mode sequence. We also calculate the weighted covariance matrices using Eq. (46) or (48).

• 3.

  Spectral fit: the spectral-fit method is applied to either the averaged mode periodograms (modal scheme) or the individual weighted OPD periodograms (OPD scheme) to get the fit parameters.

• 4.

  Filling of the Kalman matrices: the fitted parameters are used to define the system matrices, and to calculate the gain G_∞.

Tracking phase

For each elementary fringe exposure:

• 1.

  Weighting: the residual OPD estimates for the six baselines,obtained from the phase sensor, are combined in a weighted wayto give optimal residual mode estimates (modal scheme) or OPDestimates (OPD scheme). The weighting factors are calculatedfrom instantaneous estimates of the measurement errors. Wenote that this step requires a singular-value decomposition tocalculate the weighted generalized inverse of M. We calculate the associated instantaneous measurement covariances Σ_Q,w,W,n or Σ_w,W,n (Eq. (46) or (48) with instantaneous weights).

• 2.

  Kalman step: we apply the Kalman filter, using the instantaneous gain definitions (Eq. (49) or (50)).

• 3.

  Command calculation: we calculate the actuator commands.

7. Simulations for VLTI/GRAVITY

To test the performance of our Kalman-filter control, we construct a simple simulator of the control process in a fringe tracker, which also allows us to distinguish between the different control schemes. The simulations are performed in function of the key science case of the four-telescope beam combiner GRAVITY: the astrometric observation of the Galactic-center neighborhood. As described in Gillessen et al. (2010), stable fringe tracking of a magnitude K = 10 star is required, with a specification of 300   nm for the residual OPD.

More complete simulations of the fringe tracking process for GRAVITY, including frame-by-frame acquisition, are postponed to a future paper (Choquet et al., in prep.). In the same paper, a comparative study with classical control algorithms will be presented and loop parameters will be optimized.

7.1. Disturbance simulation

Data is simulated based on model power spectra for the associated phenomena. Since turbulence and vibration components are assumed to be mutually independent, we can sum individual perturbation sequences to assemble complete data sequences.

For each data sequence, we take a number of vibration components under the form of randomly excited damped oscillators. The frequencies chosen are inspired by OPD power spectra obtained with the VLTI/PRIMA fringe tracker, shown in Sahlmann et al. (2009). The total rms deviation per piston is chosen between 200 and 240   nm, which is a realistic estimate of the current vibration level at VLTI (Poupar et al. 2010). The vibration profiles used for each telescope are shown in Fig. 4 (superimposed on the turbulence profile), with some specifications in Table 1.

Considering the turbulence contribution, Conan et al. (1995) derive the asymptotic laws for single turbulence-layer OPD spectra $\mathit{S}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{\propto }\{\begin{array}{c}\\ {\mathit{f}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}& \mathrm{for} \mathrm{0.2}\hspace{0.17em}\mathit{v}\mathit{/}\mathit{B}\mathit{>}\mathit{f,}\\ {\mathit{f}}^{\mathrm{-}\mathrm{8}\mathit{/}\mathrm{3}}& \mathrm{for} \mathrm{0.2}\hspace{0.17em}\mathit{v}\mathit{/}\mathit{B}\mathit{<}\mathit{f}\mathit{<}\mathrm{0.3}\hspace{0.17em}\mathit{v}\mathit{/}\mathit{D,}\\ {\mathit{f}}^{\mathrm{-}\mathrm{17}\mathit{/}\mathrm{3}}& \mathrm{for} \mathrm{0.3}\hspace{0.17em}\mathit{v}\mathit{/}\mathit{D}\mathit{<}\mathit{f}\end{array}$(53)(pure Kolmogorov turbulence, i.e. with infinite outer scale), where v is the wind speed parallel to the baseline, B is the baseline, and D is the diameter of the telescopes. The first two power laws in Eq. (53) are clearly observable in OPD spectra obtained with the Mark III interferometer (Colavita et al. 1987). In this case, the steep f^−17/3 cut-off is not observed because of to the small telescope diameters (high cut-off frequency). Nevertheless, the f^−17/3 law, which is also expected for mono-layer tip-tilt spectra, is also in practice not observed on large telescopes. The tip-tilt power spectra of the AO-system VLT/NAOS, used by Meimon et al. (2010) to verify the spectral identification method (MPFK method, Sect. 4.3), are well-modeled by an AR(2)-model spectrum (lacking a very steep cut-off). In addition, OPD spectra obtained with the VLTI/PRIMA fringe tracker using the ATs (D = 1.8   m) are compatible with an f ^− 8/3 power law up to 200 Hz, where a flat noise tail takes over (Sahlmann et al. 2009). We thus conclude that, in real-life situations (multi-layer turbulence, integrated ${\mathit{C}}_{\mathit{N}}^{\mathrm{2}}$-profiles, multi-directional wind), the f^−17/3 law is not observable/observed.

We therefore decide to use the OPD-spectrum model in Eq. (53), but discard the f^−17/3 cut-off. A very similar model is used by Gitton (2010) in the context of VLTI instrumentation. In Table 1, we present the values assigned to the parameters, which are representative of the observing conditions at the VLTI. Each piston sequence is scaled to an rms deviation of 10 μm, such that the resulting OPD sequence has an rms value of  ~14   μm ($\mathrm{10}\sqrt{\mathrm{2}}\hspace{0.17em}\mathit{\mu }$m). The theoretical turbulence profile is shown in Fig. 4.

7.2. Flux dropouts and measurement noise

In addition, we need to introduce the phenomenon of varying flux in the simulations, which manifests itself as variable measurement noise. This means that, in parallel to the simulated OPD sequences, we need to simulate a time sequence of instantaneous noise characteristics, and add a random noise-realization sequence to the noise-free OPD sequence. The measurement error ${\mathit{\sigma }}_{\mathit{w}}^{\mathit{ij}}$ on a baseline ij is inversely proportional to the mean observed visibility V^ij and the S/N ρ^ij on that baseline, where ${\mathit{V}}^{\mathit{ij}}\mathrm{=}\frac{\mathrm{2}\sqrt{{\mathit{N}}^{\mathit{i}}{\mathit{N}}^{\mathit{j}}}}{{\mathit{N}}^{\mathit{i}}\mathrm{+}{\mathit{N}}^{\mathit{j}}}$(54)and${\mathit{\rho }}^{\mathit{ij}}\mathrm{=}\frac{{\mathit{N}}^{\mathit{i}}\mathrm{+}{\mathit{N}}^{\mathit{j}}}{\sqrt{{\mathit{N}}^{\mathit{i}}\mathrm{+}{\mathit{N}}^{\mathit{j}}\mathrm{+}\mathrm{4}\mathrm{·}{\mathrm{RON}}^{\mathrm{2}}}}\mathrm{·}$(55)Here, RON is the read-out noise per pixel and Nⁱ is the number of coherent photons arriving from telescope i for a single phase measurement. The factor 4, on the other hand, is due to the phase sensing for GRAVITY being based on 4-pixel measurements (ABCD algorithm, see Shao & Staelin 1977). For simplicity, we only consider OPD estimation using phase-delay estimation. Phase-unwrapping techniques for resolving the 2π ambiguity in measured fringe phases, for example using a combination of phase-delay and group-delay algorithms (e.g. Blind et al. 2011), are postponed to future work

In terms of the number of coherent photons N arriving per exposure on one telescope, we have ${\mathit{N}}^{\mathit{i}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{·}\frac{\mathrm{1}}{\mathrm{5}}\hspace{0.17em}\mathrm{𝒯}\hspace{0.17em}\mathit{N,}$(56)where is the total coherent throughput (AO-corrected atmospheric turbulence, delay-line system and GRAVITY). The factor $\frac{\mathrm{1}}{\mathrm{3}}$ comes from the splitting of arriving photons on the three baselines corresponding to the telescope, while $\frac{\mathrm{1}}{\mathrm{5}}$ refers to the dispersion of fringes on five spectroscopic channels (needed for the group-delay estimation, not considered in the current simulations).

The variability in the coherent throughput is simulated based on a model for the propagating tip-tilt errors in the AO system, thus simulating improper fiber injection. Basically, we assume that $\mathrm{𝒯}\mathrm{=}{\mathrm{𝒯}}_{\mathrm{0}}\hspace{0.17em}\exp (\mathrm{-}{\mathit{\phi }}^{\mathrm{2}})\mathit{.}$(57)The quantity φ is the relative offset of proper fringe injection (the ratio of tip-tilt offset to the fiber mode field radius). On the basis of tip-tilt measurements at the VLTI, the temporal power spectrum for φ is derived to be ${\mathit{S}}^{\mathit{\phi }}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{\propto }\{\begin{array}{c}\\ \log \mathit{f}& \mathrm{for} \mathrm{2}\hspace{0.17em}\mathrm{Hz}\mathit{<}\mathit{f}\mathit{<}\mathrm{8}\hspace{0.17em}\mathrm{Hz}\mathit{,}\\ \mathrm{-}\log \mathit{f}& \mathrm{for} \mathrm{8}\hspace{0.17em}\mathrm{Hz}\mathit{<}\mathit{f}\mathit{<}\mathrm{50}\hspace{0.17em}\mathrm{Hz}\mathit{,}\hspace{0.17em}\mathit{f}\ne \mathrm{18.1}\hspace{0.17em}\mathrm{Hz}\mathit{,}\\ \mathit{\delta }\mathrm{\left(}\mathit{f}\mathrm{\right)}& \mathrm{for} \mathit{f}\mathrm{=}\mathrm{18.1}\hspace{0.17em}\mathrm{Hz}\end{array}$(58)(Gitton, priv. comm.). The total rms tip-tilt deviation is scaled to 14.6 mas, of which 5 mas is present in the vibration component. A deeper discussion of this power spectrum would lead us too far for current purposes, and will be presented in future work (Choquet et al., in prep.). An extract of a simulated normalized throughput sequence is shown in Fig. 5.

Since phase delays estimated on the five spectroscopic channels can be combined, the estimator for the measurement error is derived to be⁷${\mathit{\sigma }}_{\mathit{w}}^{\mathit{ij}}\mathrm{=}\frac{\mathit{\lambda }}{\mathrm{2}\mathit{\pi }}\sqrt{\frac{\mathrm{2}}{\mathrm{5}}}\frac{\sqrt{{\mathit{N}}^{\mathit{i}}\mathrm{+}{\mathit{N}}^{\mathit{j}}\mathrm{+}\mathrm{4}\mathrm{·}{\mathrm{RON}}^{\mathrm{2}}}}{\mathrm{2}\sqrt{{\mathit{N}}^{\mathit{i}}{\mathit{N}}^{\mathit{j}}}}$(59)(Houairi et al. 2008). The parameters in this and the above equations are specified in Table 1, and are representative of the astrometric observation of the Galactic center, which is the principal science case of GRAVITY. At maximal throughput (), we derive a coherent photon number per exposure of Nⁱ = 20 and a corresponding measurement error of ${\mathit{\sigma }}_{\mathit{w}}^{\mathit{ij}}\mathrm{=}\mathrm{68}\hspace{0.17em}$nm.

7.3. Turbulence-only simulations

Before considering simulations with the complete set of expected OPD perturbations (turbulence and longitudinal vibrations), we perform a number of fringe-tracking runs including only the effect of turbulence (but with flux dropouts). The simulation process is started by generating a sequence of 2000 POL-data points. To keep the focus of our simulations on investigating the principle and performance of Kalman-filter fringe tracking, we chose to construct POL data in a simple way by adding measurement errors (Eq. (59)) to a turbulence OPD sequence. In complete and fully realistic simulations of the fringe tracking, the simulated POL data should be acquired in a frame-by-frame simulation with the phase-sensing algorithm. As mentioned before, this problem will be addressed in our future work (Choquet et al., in prep.).

Following the operation scheme defined in Sect. 6, we simulate 200 full fringe-tracking runs of 100 s for both the modal and the OPD Kalman filters. To test our definition of the instantaneous Kalman gains when handling flux dropouts, we distinguish between simulations without modifying the gains (F, fixed) and with instantaneous gains (I).

Once the fit parameters are obtained and the Kalman matrices are generated, both recombination schemes can be tested. The fringe tracking is initiated with an empty state vector ($\hat{{x}}{}_{\mathrm{0}\mathrm{|}\mathrm{-}\mathrm{1}}\mathrm{=}\mathrm{0}$). The convergence of the Kalman filter to the actual state of the system typically takes a few to a few tens of tracking points. We reject this transition region when calculating the residual OPD.

The results of the turbulence-only simulations are shown in Fig. 6. In this plot, all residual rms OPD values (6 per run) are grouped in blocks of 2 nm, giving us a measure of the distribution of residual OPD values. We found that the performance of both control schemes (modal and OPD) in both gain settings (fixed and instantaneous) is similar. The residuals are distributed with some degree of asymmetry around a nominal residual value in the range of 125–145 nm. The total tracking residuals of 100-s simulations are thus only about twice as large as the measurement error ${\mathit{\sigma }}_{\mathit{w}}^{\mathit{ij}}$ at maximum throughput for a single OPD observation (68   nm, see previous section). This indicates that the tracking performance of the Kalman-filter schemes is good. For both filters (modal versus OPD), the gain modification slightly reduces the residual OPDs, by about 3–4 nm (about 5–6% in energy).

7.4. Full perturbation simulations

We now turn to simulations based on the full power spectrum in Fig. 4. We put a maximum of 40 on the number of vibrations that can be identified during the mode identification, and 20 per baseline for the OPD identification. These numbers correspond to a maximum matrix size of 240 × 240 for Σ_∞.

As in the case of the turbulence-only simulations, we used 2000-point data sequences for the identification. The result of a typical spectral fit is shown in Fig. 7. By definition, the average mode spectrum contains all vibration components. This is clearly apparent in the upper plot, where the number of components to identify is about twice as high (a single baseline has the vibrations of only two telescopes). Another thing to note is that the statistical spread of the average modal periodogram is indeed lower than that of the OPD periodogram.

A total of 200 full simulations of identification and 100-s fringe-tracking run were performed. The corresponding distributions of OPD residuals are shown in Fig. 8, organized into blocks of 10 nm. The addition of vibration perturbations has introduced a considerable spread in the residual OPD values. In addition, unlike in the turbulence-only simulations, a clear difference in performance between the modal and the OPD scheme is found. The residual rms value for the modal scheme is about 290   nm, on average, where about 30% of the measurements are above 300   nm. The OPD scheme, on the other hand, has an average of about 240   nm for the rms OPD value, and only 6% of the residuals surpass 300   nm. When choosing between fixed or instantaneous gains, there is no apparent significant difference. We describe these observations in greater detail in Sect. 7.5.

Since the vibration perturbations considered are on the level of 200–240 nm, the corresponding OPD vibration perturbations are roughly 300 nm (e.g. $\mathrm{220}\sqrt{\mathrm{2}}\hspace{0.17em}\mathrm{nm}\mathrm{=}\mathrm{311}\hspace{0.17em}$nm). Comparing the smallest residuals of turbulence-only simulations (~125 nm) to the best residuals of the current simulations (~200 nm), our simulations indicate that as much as 75% of the vibration energy can be filtered out by our Kalman filters (we calculate that (200² − 125²)/311² = 25% of the vibration energy is left).

7.5. Discussion

The above simulations give a first quantitative result for the Kalman-filter schemes developed in this work. We have illustrated the operation of the spectral fit on realistic perturbation profiles, and have found the performance one would expect from the modal and OPD control schemes.

When considering our definition of instantaneous gains, it is clear that no significant differences between fixed and instantaneous gains are seen, for the considered flux perturbation sequences. One conclusion that might be drawn is that the instantaneous gains do not improve the fringe-tracker performance. To verify that this would be the case, we performed simulations without the flux dropouts, for which we added the results to Fig. 8. Surprisingly, although there is some apparent difference, the global pattern in the distribution of OPD residuals seems to be rather unaffected by the flux dropouts. At this point, it might be expected that the Kalman-filter based control scheme is not limited by the expected fluctuations in flux. Future comparison to other control strategies would be required to determine whether this property is exclusive to Kalman filtering. Although the instantaneous-gain definitions do not change the Kalman-filter performance for the considered flux perturbations, we assume that the gains are still useful. When poorer observing conditions prevail, causing longer flux dropouts for instance, the gain can still be expected to limit the risk of improper tracking. In this respect, using instantaneous gains might be very advantageous to the fringe tracking on the ATs, which lack a higher-order AO correcting system. Additionally, when lower fluxes are assumed for the observed object, the gain modification becomes useful.

The main observation of our simulations is that the OPD Kalman filter provides the best results. This is somewhat unexpected, since the modal scheme has the most reliably estimated power spectrum to model, and has the optimal number of variables to track. One of the main disadvantages of the modal scheme, however, is that all vibrations are fitted at once. As a result, the periodogram fit is more difficult:

• A higher fraction of the periodogram points are part of a vibrationpeak, thus fewer points are available to estimate the turbulencecomponent. Since the latter component is the most importantperturbation, a good estimation of its associated power spectrumis essential.

• Vibration peaks place a much stronger bias on the estimation of the turbulence component. The reason is that the statistics of an average mode periodogram are “narrower” (i.e. have a smaller spread, see the discussion in Appendix B), hence vibration peaks have a stronger influence on the turbulence estimation than in the case of an OPD-periodogram fit.

• The more peaks that are present in the periodogram, the higher the probability that the peaks will overlap. Owing to the insufficient resolution of the periodograms, overlapping peaks will be fitted as a single broader peak, which is more difficult to correct (more transient behavior).

All of these properties might lower the quality of the power-spectrum estimation for the modal scheme, hence lower the Kalman filter performance.

Another property that might have been incorporated to explain the performance differences between modal and OPD control is that flux dropouts have different implications for the schemes. We recall that a low level of flux at one telescope can significantly affect the quality of all modes, whereas three of the six OPDs are well-measured in the OPD scheme. However, since simulations with flux dropouts give similar results to simulations without, the performance difference cannot be attributed to flux dropouts. Temporally losing information about all mode observables thus does not significantly influence the tracking performance of the Kalman filter. In contrast, for non-predictive control algorithms applied to the modes, losing observables can be expected to be a critical limitation.

When considering the actual shape of the distributions in Fig. 8, we find a quite large spread in the residual OPD values. Not taking into account local statistical jitter, the distributions have a clear global maximum with some “leaking” to higher residuals (~10–15%). We ascribe these outliers to improper estimation of the power spectra of the perturbation components. Improvements to the identification method could indeed be made to make it more robust (e.g. optimal grid construction, reconsidering the number of identification points).

As a final note, it cannot be ruled out that the performance difference is related to the conceptual difference between the Kalman filters. For a vibration acting on one telescope, the three baselines that are not associated remain unaffected from this vibration. In the modal scheme, however, the vibration needs to be filtered in all modes. It might be this property of indirectly affecting the other telescopes (via the modal filtering) that causes the lower average performance of the modal scheme. More tests are needed to verify this reasoning.

7.6. Performance with respect to GRAVITY

The main conclusion of these simulations is that the average performance of the Kalman-filter control schemes is compatible with the performance required for GRAVITY, specified at 300   nm (Gillessen et al. 2010). For the OPD scheme, at least 90% of the performed tracks should attain this stability level, and a nominal performance of λ/10 is within reach. This makes the latter scheme clearly the most suitable of the two.

The most important working point to reach this stability level for our Kalman-filter schemes is optimizing the operation of the identification method. In this respect, a decent implementation should limit both the risk of low-quality identification and the calculation time. For the latter point, an optimal adaptation of the parameter grid to the typical prevailing disturbances at the VLTI should be defined (Chouquet et al., in prep.).

Compared to the simulations done in Choquet et al. (2010), based on a classical integrator controller, our simulations show that the nominal Kalman-filter performance might be compatible or better. In any case, we should immediately nuance this comparison, since the assumed parameters and perturbation profiles are rather different than the ones used here. For example, the simulated turbulence profiles used here exclude the unobserved f^−17/3 part and a higher and more realistic vibration contribution is assumed. Our future work will present a firmer base for comparing Kalman-filter fringe tracking to classical control, something that is beyond the reach of the current work.

8. Summary

In the framework of new-generation interferometric instruments, we have designed and analyzed a Kalman control scheme for fringe tracking. The starting point of this work was the work done in AO control, in particular a control scheme proposed by Le Roux et al. (2004) based on Kalman filtering, which is an optimal data processing algorithm. Apart from the usual correction for turbulence disturbance, a Kalman filter allows us to correct for vibration disturbances, present in VLTI observations under the form of longitudinal vibrations.

We started by formulating a control scheme for two-telescope fringe tracking, based on a parametric model for the fringe-tracking disturbances. To estimate the involved parameters, we included an identification method based on the characterization of perturbation power spectra (Meimon et al. 2010). As a direct application to multi-telescope fringe tracking, we extended the fringe-tracking formalism to four-telescope control schemes. We note that our approach can easily be generalized to n-telescope fringe tracking.

Two Kalman-filter control schemes are considered in the four-telescope case: a modal-based scheme, where three modes (linear combinations of pistons) are tracked in a single three-component Kalman filter, and an OPD-based scheme, where six OPDs are tracked using six Kalman filters. The advantage of the former is that only one spectral fit needs to be performed: the three modal periodograms can be averaged. To ensure that the control schemes are unaffected by low S/N baselines and low fluxes at the level of the telescopes (flux dropouts), we included a weighting step and a modification to the gain.

Simulations performed in the context of the four-telescope beam combiner VLTI/GRAVITY allowed to test the two control schemes and the concept of instantaneous gains. Two important observations that follow are:

• 1.

  Despite being operationally heavier, the OPD control schemeprovides the best results. We (partly) attribute the lowerperformance of the modal control to the more difficultidentification (all vibration parameters need to be extracted atonce). There might also be a conceptual performance limitationrelated to Kalman filtering of the modes, in that all perturbationsneed to be corrected in all modes.

• 2.

  For the considered flux perturbations, the instantaneous gains do not significantly improve the fringe-tracking performance.

When considering the latter point, comparison to simulations without the flux dropouts have indicated that Kalman-filter control seems to be rather unaffected by the considered flux perturbations. Yet, the instantaneous gain might prove useful when worse observing conditions prevail (e.g. fainter tracking source, longer flux dropouts).

As far as the science case of VLTI/GRAVITY is concerned, both schemes show nominal performances that meet the specified performances. The results open perspectives for attaining even higher quality residuals (λ/10 performance). The most important working point for the future is optimizing the identification method.

---

Acknowledgments

We would like to thank Pierre Fedou for useful interactions on the design of VLTI/GRAVITY. We are grateful to Frank Eisenhauer, the editor (Thierry Forveille) and the anonymous referee for comments and suggestions that helped improving this manuscript.

References

Appendix A: Derivation of the modal state-space model

The derivation below is based on the state-space model in Sect. 2. We use a natural intermediate state-space model to develop the modal model.

A.1. Canonical state-space model

Before going to the complete set of four telescopes, we recast the problem of a two-telescope interferometer (Sect. 2) in a more useful form. The result allows a direct generalization to an n-telescope interferometer.

A.1.1. Two-telescope interferometer in a piston scheme

In two dimensions, piston disturbances and associated OPD measurements can be written like P = (P⁰      P¹)^⊤ and (trivially) OPD = (OPD⁰¹). Finally, the matrix M that governs the conversion P → OPD is $\mathrm{M}\mathrm{=}(\mathrm{-}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1})\mathit{.}$(A.1)We define a new state-space model in terms of fundamental piston disturbances P (rather than the phase disturbances ϕ in Sect. 2). For instance, a vibration on the level of piston P^k (k = 0,1) can be written as ${\mathit{P}}_{\mathit{n}\mathrm{+}\mathrm{1}}^{\mathit{k,}\mathrm{vib}}\mathrm{=}{\mathit{a}}_{\mathrm{1}}^{\mathit{k,}\mathrm{vib}}{\mathit{P}}_{\mathit{n}}^{\mathit{k,}\mathrm{vib}}\mathrm{+}{\mathit{a}}_{\mathrm{2}}^{\mathit{k,}\mathrm{vib}}{\mathit{P}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathit{k,}\mathrm{vib}}\mathrm{+}{\mathit{v}}_{\mathit{n}}^{\mathit{k,}\mathrm{vib}}\mathit{,}$(A.2)where we note that we simply include one extra index, k, to refer to the considered piston. In this state-space model, new state vectors x are defined as ${{x}}_{\mathit{n}}\mathrm{=}({{x}}_{\mathit{n}}^{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{{x}}_{\mathit{n}}^{\mathrm{1}}{)}^{\mathrm{\top }}\mathit{,}$(A.3)where, taking e.g. two vibrations on piston 0 and one on piston 1, $\begin{array}{ccc}{{x}}_{\mathit{n}}^{\mathrm{0}}& \mathrm{=}& ({\mathit{P}}_{\mathit{n}}^{\mathrm{0}\mathit{,}\mathrm{tur}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{0}\mathit{,}\mathrm{tur}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}_{\mathit{n}}^{\mathrm{0}\mathit{,}\mathrm{vib}\hspace{0.17em}\mathrm{1}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{0}\mathit{,}\mathrm{vib}\hspace{0.17em}\mathrm{1}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}_{\mathit{n}}^{\mathrm{0}\mathit{,}\mathrm{vib}\hspace{0.17em}\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{0}\mathit{,}\mathrm{vib}\hspace{0.17em}\mathrm{2}}{)}^{\mathrm{\top }}\mathit{,}\\ {{x}}_{\mathit{n}}^{\mathrm{1}}& \mathrm{=}& ({\mathit{P}}_{\mathit{n}}^{\mathrm{1}\mathit{,}\mathrm{tur}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{1}\mathit{,}\mathrm{tur}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}_{\mathit{n}}^{\mathrm{1}\mathit{,}\mathrm{vib}\hspace{0.17em}\mathrm{1}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{P}}_{\mathit{n}\mathrm{-}\mathrm{1}}^{\mathrm{1}\mathit{,}\mathrm{vib}\hspace{0.17em}\mathrm{1}}{)}^{\mathrm{\top }}\mathit{.}\end{array}$The full equation of state can then be written in block form as $\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left(\begin{array}{c}\\ \\ \end{array}\right)}\\ {{x}}_{\mathit{n}\mathrm{+}\mathrm{1}}\end{array}\mathrm{=}\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left(\begin{array}{c}\\ & \\ & \end{array}\right)}\\ \mathrm{A}\end{array}\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left(\begin{array}{c}\\ \\ \end{array}\right)}\\ {{x}}_{\mathit{n}}\end{array}\mathrm{+}\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left(\begin{array}{c}\\ \\ \end{array}\right)}\\ {{v}}_{\mathit{n}}\end{array}\mathit{,}$(A.4)where block diagonal matrices Aⁱ are defined exactly as before (for each fundamental piston). We assume that the system noise v is composed of independent white Gaussian noise components, hence the covariance matrix Σ_v is still diagonal.

Given that we consider the same two-telescope interferometer as before, the observable is still the same: the residual OPD ${\mathit{y}}_{\mathit{n}}^{\mathrm{01}}$ (note that we give the index “01” to refer to the baseline on which we work). The new state vector requires an updated definition of the system matrix C, which we denote $\overline{\mathrm{C}}$. An appropriate definition is $\overline{\mathrm{C}}\mathrm{=}\mathrm{(}\mathrm{-}{\mathrm{C}}^{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathrm{C}}^{\mathrm{1}}\mathrm{)}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathrm{)}\left(\begin{array}{c}\\ {\mathrm{C}}^{\mathrm{0}}& \mathrm{0}\\ \mathrm{0}& {\mathrm{C}}^{\mathrm{1}}\end{array}\right)\mathrm{=}\mathrm{M}\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left(\begin{array}{c}\\ & \\ & \end{array}\right)}\\ \mathrm{C}\end{array}\mathit{,}$(A.5)where C⁰ and C¹ are row matrices defined as before. For example, the corresponding C⁰ and C¹ to Eq. (A.3) are ${\mathrm{C}}^{\mathrm{0}}\mathrm{=}\mathrm{(}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathrm{)} \mathrm{and} {\mathrm{C}}^{\mathrm{1}}\mathrm{=}\mathrm{(}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\mathrm{)}\mathit{.}$(A.6)In accordance with the definition of the matrix A in Eq. (A.4), we have defined C in Eq. (A.5) as the block-diagonal matrix of matrices Cⁱ. Considering the piston control, we wish to calculate the commands to be applied directly to the pistons. Command vectors are hence of length two (we consider two actuators). The full measurement equation in its new form is then ${\mathit{y}}_{\mathit{n}}^{\mathrm{01}}\mathrm{=}\mathrm{M}\hspace{0.17em}\mathrm{C}{{x}}_{\mathit{n}}\mathrm{-}\mathrm{M}\hspace{0.17em}{{u}}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{+}{\mathit{w}}_{\mathit{n}}^{\mathrm{01}}\mathit{.}$(A.7)

A.1.2. Four-telescope interferometer in piston scheme

The equation of state in Eq. (A.4) in now trivially generalized to the four-double equivalent $\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left(\begin{array}{c}\\ \\ \\ \\ \end{array}\right)}\\ {{x}}_{\mathit{n}\mathrm{+}\mathrm{1}}\end{array}\mathrm{=}\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)}\\ \mathrm{A}\end{array}\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left(\begin{array}{c}\\ \\ \\ \\ \end{array}\right)}\\ {{x}}_{\mathit{n}}\end{array}\mathrm{+}\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left(\begin{array}{c}\\ \\ \\ \\ \end{array}\right)}\\ {{v}}_{\mathit{n}}\end{array}\mathit{.}$(A.8)We still assume that the covariance matrix Σ_v is diagonal, i.e. the vectors ${{v}}_{\mathit{n}}^{\mathit{k}}$ (k = 0,1,2,3) have independent white Gaussian noise components. The fact that fringes are tracked on six baselines leads to the use of residual OPD vectors y and measurement noise vectors w, where ${y}\mathrm{=}({\mathit{y}}^{\mathrm{01}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{y}}^{\mathrm{02}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{y}}^{\mathrm{03}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{y}}^{\mathrm{12}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{y}}^{\mathrm{13}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{y}}^{\mathrm{23}}{)}^{\mathrm{\top }}\mathit{,}$(A.9)and analogously for w. The system matrix C and the vectors u_n also have their four-dimensional equivalent, and the full measurement equation finally reads ${{y}}_{\mathit{n}}\mathrm{=}\mathrm{M}\hspace{0.17em}\mathrm{C}\hspace{0.17em}{{x}}_{\mathit{n}}\mathrm{-}\hspace{0.17em}\mathrm{M}\hspace{0.17em}{{u}}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{+}{{w}}_{\mathit{n}}\mathit{.}$(A.10)As we did with Σ_v, the measurement error covariance matrix Σ_w is assumed to be diagonal.

A.1.3. Problems with the state-space model

The combination of Eqs. (A.8) and (A.10) forms a very natural full state-space model of a four-telescope six-baseline fringe tracking system. On the one hand, it describes in a natural way the disturbances on the level of the four apertures, under the reasonable assumption that these are decoupled. On the other hand, the measurement equation describes how the six measured residual OPDs are the result of the contributions of the four fundamental pistons, and of the correction that is applied by the four actuator systems (piezos + delay lines). Owing to the very natural description given by this model, we can refer to it as the canonical state-space model. It seems that we have found an appropriate input model to build a Kalman filter.

However, there is a serious difficulty of the canonical model: it is impossible to identify its parameters. In Sect. 3, the parameter identification method that is introduced is based on applying a spectral fit to sequences of (pseudo-)open loop measurements. Here, however, the matrix M couples all information about the evolution of pistons, in a non-trivial way. In essence, the problem is that absolute piston information cannot be unambiguously recovered from observations. Difficulties with the system identification for a piston-based state-space model could hence have been expected. However, the state-space model defined above now forms a perfect starting point for the construction of a modal-based state-space model.

A.2. Modal state-space model

We now consider how the canonical state-space model can be put into a modal form, in the spirit of what was done in Sect. 4.1. This direct derivation of a modal state-space model is more reliable than simply applying one Kalman-filter controller per mode. In particular, although the modes are orthogonally defined, they have coupling terms that are related to there being linear combinations of pistons/OPDs. The power of the Kalman filter is that it is constructed to deal with vector quantities and covariance matrices which describe coupling between different terms.

Modal equation of state.

For the new equation of state, we introduce a piece of new notation. We use the tilde sign (~) to denote extended vector and matrix quantities with respect to the previously introduced versions. First we introduce the extended state-space vectors ${\begin{array}{}{\mathrm{􏽥}}^{}\\ {x}\end{array}}_{\mathit{n}}$ as ${\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}}_{\mathit{n}}\mathrm{=}\mathrm{(}\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}\begin{array}{}\mathrm{0}\\ \mathit{n}\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}\begin{array}{}\mathrm{1}\\ \mathit{n}\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}\begin{array}{}\mathrm{2}\\ \mathit{n}\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}\begin{array}{}\mathrm{3}\\ \mathit{n}\end{array}{\mathrm{)}}^{\mathrm{\top }}\mathit{,}$(A.11)where e.g. $\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}\begin{array}{}\mathrm{2}\\ \mathit{n}\end{array}\mathrm{=}({{0}}^{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{{0}}^{\mathrm{1}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{{x}}_{\mathit{n}}^{\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{{0}}^{\mathrm{3}}{)}^{\mathrm{\top }}\mathit{,} \mathrm{(}{{0}}^{\mathit{j}}\mathrm{=}\mathrm{0}\hspace{0.17em}{{x}}_{\mathit{n}}^{\mathit{j}}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathit{j}\mathrm{=}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{3}\mathrm{)}\mathit{.}$(A.12)In words, the components $\begin{array}{}{}^{\mathrm{􏽥}}\\ {x}\end{array}\begin{array}{}\mathit{i}\\ \mathit{n}\end{array}$ consist of four zero-columns of which the kth zero-column has the same length as the vector ${{x}}_{\mathit{n}}^{\mathit{k}}$ (k = 0,1,2,3), but the ith zero-column is replaced by ${{x}}_{\mathit{n}}^{\mathit{i}}$ itself. It follows that, by definition, all $\begin{array}{}{}^{\mathrm{􏽥}}\\ {x}\end{array}\begin{array}{}\mathit{i}\\ \mathit{n}\end{array}$ have the same length (unlike the vectors ${{x}}_{\mathit{n}}^{\mathit{i}}$), and we define $\mathit{r}\mathrm{\equiv }\mathrm{lengthofvectors}\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}\begin{array}{}\mathit{i}\\ \mathit{n}\end{array}\mathit{.}$(A.13)In exactly the same way as in Eq. (A.11), we define the vectors ${\begin{array}{}{\mathrm{􏽥}}^{}\\ {v}\end{array}}_{\mathit{n}}$ as extensions to v_n, and the corresponding covariance matrix is written as ${\mathrm{\Sigma }}_{\begin{array}{}{\mathrm{􏽥}}^{}\\ \mathit{v}\end{array}}$. Finally, we define the matrix $\begin{array}{}{\mathrm{􏽥}}^{}\\ \mathrm{A}\end{array}$ as $\begin{array}{}\mathrm{􏽥}\\ \mathrm{A}\end{array}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{A}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{A}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{A}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{A}\end{array}\right)\mathit{,}$(A.14)i.e. a block diagonal matrix with the square system matrix A on the diagonal. A completely equivalent equation of state to the one in Eq. (A.8) is then given by ${\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}}_{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{=}\begin{array}{}\mathrm{􏽥}\\ \mathrm{A}\end{array}\hspace{0.17em}{\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}}_{\mathit{n}}\mathrm{+}{\begin{array}{}\mathrm{􏽥}\\ {v}\end{array}}_{\mathit{n}}\mathit{.}$(A.15)To pass to a modal equation of state, we follow an approach similar to that used in the Q-mode definition Q = V^⊤   P. We define a (still unitary) block-matrix version of the matrix V $\begin{array}{ccc}{\mathrm{V}}_{\square }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\begin{array}{c}\\ & & & \\ & & & \\ & & & \\ & & & \end{array}\right)\mathit{,}& & \end{array}$(A.16)where “1” denotes a unit r × r matrix. The modal state vector is then defined as $\begin{array}{ccc}{{x}}_{\mathit{Q,n}}\mathrm{\equiv }{\mathrm{V}}_{\square }^{\mathrm{\top }}\hspace{0.17em}{\begin{array}{}\mathrm{􏽥}\\ {x}\end{array}}_{\mathit{n}}\mathit{.}& & \end{array}$(A.17)Multiplying Eq. (A.15) by and using the equality , the modal equation of state finally reads ${{x}}_{\mathit{Q,n}\mathrm{+}\mathrm{1}}\mathrm{=}\begin{array}{}\mathrm{􏽥}\\ \mathrm{A}\end{array}\hspace{0.17em}{{x}}_{\mathit{Q,n}}\mathrm{+}{{v}}_{\mathit{Q,n}}\mathit{,}$(A.18)where we have defined the mode system noise $\begin{array}{ccc}{{v}}_{\mathit{Q,n}}\mathrm{\equiv }{\mathrm{V}}_{\square }^{\mathrm{\top }}\hspace{0.17em}{\begin{array}{}\mathrm{􏽥}\\ {v}\end{array}}_{\mathit{n}}\mathit{.}& & \end{array}$(A.19)The covariance matrix of v_Q is . Using some numerical examples, it can be shown that Σ_Q,v has a quadruple diagonal structure, i.e. four equal diagonal blocks, but also possesses off-diagonal terms.

Modal observation equation.

We multiply the observation equation by H (Eq. (31)), which leads to ${{y}}_{\mathit{Q,n}}\mathrm{\equiv }\mathrm{H}\hspace{0.17em}{{y}}_{\mathit{n}}\mathrm{=}\mathrm{R}\hspace{0.17em}{\mathrm{V}}^{\mathrm{\top }}\mathrm{C}{{x}}_{\mathit{n}}\mathrm{-}\mathrm{R}\hspace{0.17em}{\mathrm{V}}^{\mathrm{\top }}{{u}}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathrm{+}\mathrm{H}\hspace{0.17em}{{w}}_{\mathit{n}}\mathit{,}$(A.20)where R = diag(1,1,1,0). We note that we used the property H   M = R   V^⊤. We verify that this equation indeed defines the last component of y_Q,n to be zero, in accordance with our convention that Q³ ≡ 0. The defined transformation gives the observing equation in terms of the residual Q-mode observables y_Q. Next, we need to pass from the piston state vector x_n to the equivalent x_Q,n for the modes. For this, we define an extended version of the matrix C as⁸$\begin{array}{}\mathrm{􏽥}\\ \mathrm{C}\end{array}\mathrm{=}\left(\begin{array}{c}\\ {\mathrm{C}}_{\mathrm{-}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathrm{C}}_{\mathrm{-}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\mathrm{C}}_{\mathrm{-}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{,} \mathrm{where} {\mathrm{C}}_{\mathrm{-}}\mathrm{=}\mathrm{(}{\mathrm{C}}^{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathrm{C}}^{\mathrm{1}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathrm{C}}^{\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathrm{C}}^{\mathrm{3}}\mathrm{)}\mathit{.}$(A.21)It takes some pen and paper work to show that one can rewrite the first term in the right hand side of Eq.~(\ref{eq:oemodified}) as$\mathrm{R}\hspace{0.17em}{\mathrm{V}}^{\mathrm{\top }}\hspace{0.17em}\mathrm{C}\hspace{0.17em}{{x}}_{\mathit{n}}\mathrm{=}\begin{array}{}\mathrm{􏽥}\\ \mathrm{C}\end{array}\hspace{0.17em}{{x}}_{\mathit{Q,n}}\mathit{.}$(A.22)Two final definitions are the modal command and the mode measurement noise ${{u}}_{\mathit{Q,n}}\mathrm{\equiv }\mathrm{R}\hspace{0.17em}{\mathrm{V}}^{\mathrm{\top }}\hspace{0.17em}{{u}}_{\mathit{n}} \mathrm{and} {{w}}_{\mathit{Q,n}}\mathrm{\equiv }\mathrm{H}\hspace{0.17em}{{w}}_{\mathit{n}}$(A.23)(for both, the last component is 0). The full modal equation is finally given by ${{y}}_{\mathit{Q,n}}\mathrm{=}\begin{array}{}\mathrm{􏽥}\\ \mathrm{C}\end{array}\hspace{0.17em}{{x}}_{\mathit{Q,n}}\mathrm{-}{{u}}_{\mathit{Q,n}\mathrm{-}\mathrm{2}}\mathrm{+}{{w}}_{\mathit{Q,n}}\mathit{.}$(A.24)We note that the covariance matrix of w_Q is Σ_Q,w = H   Σ_w   H^⊤.

Appendix B: Details of four-telescope identification

We denote by P^y the periodogram of a POL-data sequence on a single baseline, i.e. as in Eq. (23). In addition, we define S as the power spectrum associated with the periodogram P^y. Under the assumption of Gaussian white noise excitation, it can be shown that, for fixed frequencies f, the distribution of P^y(f) asymptotically approaches an exponential distribution of expectation value S(f). In other words, for long identification sequences, the likelihood function L of P^y is close to $\mathrm{L}\mathrm{\left(}{\mathit{P}}^{\mathit{y}}\mathit{,}{\lambda }\mathrm{\right)}\mathrm{=}\underset{\mathit{f}}{\mathrm{􏽙}}\frac{\mathrm{1}}{\mathit{S}\mathrm{\left(}\mathit{f}\mathrm{\right)}}\exp \left(\mathrm{-}\frac{{\mathit{P}}^{\mathit{y}}\mathrm{\left(}\mathit{f}\mathrm{\right)}}{\mathit{S}\mathrm{\left(}\mathit{f}\mathrm{\right)}}\right)\mathit{,}$(B.1)where λ denotes all model parameters of S.

Averaging three exponential distributions leads to a Gamma-Erlang distribution of order three. Hence, denoting by $\overline{{\mathit{P}}^{\mathit{y}}}$ the averaged periodogram, the periodogram likelihood function in Eq. (B.1) is modified into $\mathrm{L}\mathrm{\left(}\overline{{\mathit{P}}^{\mathit{y}}}\mathit{,}{\lambda }\mathrm{\right)}\mathrm{=}\underset{\mathit{f}}{\mathrm{􏽙}}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathrm{3}}{\mathit{S}\mathrm{\left(}\mathit{f}\mathrm{\right)}}\right)}^{\mathrm{3}}(\overline{{\mathit{P}}^{\mathit{y}}}\mathrm{\left(}\mathit{f}\mathrm{\right)}{)}^{\mathrm{2}}\exp \left(\mathrm{-}\frac{\mathrm{3}\overline{{\mathit{P}}^{\mathit{y}}}\mathrm{\left(}\mathit{f}\mathrm{\right)}}{\mathit{S}\mathrm{\left(}\mathit{f}\mathrm{\right)}}\right)\mathrm{·}$(B.2)At a fixed frequency f, the expectation value of the averaged periodogram remains the same, whereas the standard deviation is reduced by a factor $\sqrt{\mathrm{3}}$. The smaller standard deviation is the major advantage of periodogram averaging: the power spectrum is more accurately estimated, and even low vibration peaks become significant. The spectral-fit method then requires the following three modifications:

• 1.

  The maximum-likelihood criterion needs to be adapted to thelikelihood function in Eq. (B.2).

• 2.

  The threshold for vibration-peak detection (i.e. peak-picking criterion) can be lowered, since lower peaks become significantly detectable. More concretely, we choose the detection criterion P^y(f) > 4   S(f) to detect vibration peaks (probability of 0.05% per point for false detection), rather than the detection criterion P^y(f) > 7   S(f) for non-averaged periodograms (probability of 0.1% per point for false detection).

• 3.

  A more subtle point is the characterization of the measurement-noise covariance matrix Σ_Q,w (see Eq. (36)). In the standard MPFK method, applicable to single-baseline data, the measurement error is determined from the horizontal tail of the POL-data periodogram. In modal periodograms, the information of the different baselines is mixed. As a solution, we can first extract the measurement errors of the 6 OPD periodograms. Alternatively, we can simply combine the measurement errors estimated by the phase sensor (see Eq. (59)) into a global measurement error (one for each baseline). We note that the measurement noise of the average modal periodogram is calculated by summing the diagonal of Σ_Q,w and dividing by three (i.e. averaging over the three modes).

Appendix C: Instantaneous-gain definition

The calculation of an arbitrary Kalman gain matrix G_∞ (e.g. as in Eq. (18)) is of the form ${\mathrm{G}}_{\mathrm{\infty }}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{\infty }}{\mathrm{C}}^{\mathrm{\top }}\mathrm{(}\mathrm{C}{\mathrm{\Sigma }}_{\mathrm{\infty }}{\mathrm{C}}^{\mathrm{\top }}\mathrm{+}{\mathrm{\Sigma }}_{\mathit{w}}{\mathrm{)}}^{-1}\mathit{,}$(C.1)where Σ_∞ is calculated by solving a Riccati equation. It is clear that the gain depends on a measurement covariance Σ_w. Every Kalman step should ideally involve the calculation of a new gain, to compensate for the variable measurement noise (flux dropouts). As this would be a heavy computational burden, we propose adding an artificial term to the gain.

It is clear from Eq. (C.1) that the gain is inversely related to the Σ_w (neglecting that it is also used to calculate Σ_∞). It then makes sense to define an instantaneous gain G_∞,n as a quantity of the form ${\mathrm{G}}_{\mathrm{\infty }\mathit{,n}}\mathrm{\equiv }\frac{\mathit{f}({\mathrm{\Sigma }}_{\mathit{w}})}{\mathit{f}({\mathrm{\Sigma }}_{\mathit{w,n}})}\hspace{0.17em}{\mathrm{G}}_{\mathrm{\infty }}\mathit{,}$(C.2)where f is a real function. The matrix Σ_w,n is defined as the instantaneous measurement covariance at time step n. For flux of the usual level of noise, the definition of G_∞,n implies that we are tracking using the optimal gain G_∞ (the multiplicative term is 1). Whenever the instantaneous measurement noise increases, the multiplicative factor must ensure that the gain is decreased. In the extreme case when the measurements are very noisy and hence unusable, f can be chosen such that the gain approaches 0. In that case, the Kalman filter is based completely on the deterministic part of the evolution.

For the modal and the OPD Kalman filter, we take the instantaneous-gain definitions as in Sect. 5.2. It is important to understand the difference between the modal and the OPD scheme that led to Eqs. (49) and (50). In the modal scheme, each mode estimate is limited by the noise on the worst OPD triplet (by triplet, we mean three baselines that are connected to the same telescope). This is because these modes are defined as mixes of OPDs associated with all telescopes. A first order multiplicative correction of the gain thus preferably contains a measure for all error terms, hence we consider the trace of the covariance matrix, neglecting the off-diagonal terms. For the OPD scheme, however, each Kalman filter tracks a variable in an independent way. The gain of the different Kalman filters can therefore be adapted purely as a function of the noise on the corresponding variable. Again neglecting the coupling terms, the gain associated with a baseline ij is modified by the corresponding term in the error covariance matrix. In this way, we avoid that the flux problems for one telescope propagating to the tracking on baselines not associated with that telescope.

We remark that the discussion that led us to Eqs. (49) and (50) is probably not the end of the story. It is possible that more elegant “instantaneous gain” definitions exist. Yet, these definitions would likely have a similar form to our choice.

All Tables

All Figures