Comparisons of Split-Linear Fitting of Wind Curves

The detection of slope change points in wind curves depends on linear curve-fitting. Hall and Titterington's algorithm based on smooth- ing is adapted and compared to a Bayesian method of curve-fitting. After prior spline smoothing of the data, the algorithms are tested and the er- rors between the split-linear fitted wind and the real one are estimated. In our case, the adaptation of the edge-preserving smoothing algorithm gives the same good performance as automatic Bayesian curve-fitting based on a Monte Carlo Markov chain algorithm yet saves computation time. This study is aimed at the improvement of the aircraft autopilot concep- tion process. The autopilot allows landings in bad weather conditions and must guarantee passengers safety, touchdown comfort, and precision. We study in par- ticular the influence of wind during automatic landing. We focus on the effect of the linear wind components in the last 30 seconds to show that they are a decisive factor in touchdown precision. This is achieved by comparing simulated landings with either a real wind or its piecewise linear approximation. This has led us to develop a method of split-linear fitting based on slope change detection adapted to our data. This method is similar to those proposed by Jones (1998) which aim at predicting the influence of discrete gusts on linear systems. Slope change detection is associated with edge and peak detection. Two kinds of method can be adapted : a classic one, based on smoothing and a Bayesian one based on Monte Carlo Markov Chains. Hall and Titterington (1992) pro- posed an edge-preserving smoothing algorithm. It is based on edge detection by comparisons of three smoothings. The aim of this method is to compute, for each given point, three smoothed estimates of the function, based on the data to the right, to the left and on both sides of the point. Each discontinuity is asso- ciated with a local maximum of the difference between the three fits. Wu and Chu (1993) took this algorithm and modified it using kernel smoothing instead


Introduction
This study is aimed at the improvement of the aircraft autopilot conception process.The autopilot allows landings in bad weather conditions and must guarantee passengers safety, touchdown comfort, and precision.We study in particular the influence of wind during automatic landing.We focus on the effect of the linear wind components in the last 30 seconds to show that they are a decisive factor in touchdown precision.This is achieved by comparing simulated landings with either a real wind or its piecewise linear approximation.This has led us to develop a method of split-linear fitting based on slope change detection adapted to our data.This method is similar to those proposed by Jones (1998) which aim at predicting the influence of discrete gusts on linear systems.
Slope change detection is associated with edge and peak detection.Two kinds of method can be adapted : a classic one, based on smoothing and a Bayesian one based on Monte Carlo Markov Chains.Hall and Titterington (1992) proposed an edge-preserving smoothing algorithm.It is based on edge detection by comparisons of three smoothings.The aim of this method is to compute, for each given point, three smoothed estimates of the function, based on the data to the right, to the left and on both sides of the point.Each discontinuity is associated with a local maximum of the difference between the three fits.Wu and Chu (1993) took this algorithm and modified it using kernel smoothing instead of smoothing by linear combinations of observations.Loader (1996) showed convergence properties of this kind of change point estimate adapted to the search for one discontinuity in a function.Müller and Song (1997) improved this kind of algorithm to increase the convergence rate of one discontinuity estimate.Green (1995) and Denison et al. (1998) applied reversible Monte Carlo Markov Chains and Bayesian model determination to curve-fitting, including, in particular, piecewise linear polynomial fitting.The automatic Bayesian curve-fitting algorithm proposed by Denison et al. (inspired by Green's (1995) algorithm) determines the number and the position of the slope change points, knots, by the estimation of their joint posterior distribution.
Our aim here is to derive an efficient split-linear fitting algorithm.Unlike other studies, the quality of fit is not measured by the usual residual squared error.The criteria to be met are estimation of a parsimonious model by detecting only the significant slope changes in the curve and achievement of the same landing properties with the real wind and with the split-linear approximation.This is an important feature of this paper, which especially involves adaptations of the algorithms which consist in prior smoothing of the data and algorithm refinements.
A variation on Hall and Titterington's (1992) method is computed on samples of noise-free wind data measured during flight tests.The results are compared with those obtained with the automatic Bayesian curve-fitting algorithm.Denison et al. (1998) algorithm is taken here as benchmark but considered as too computational expensive for real application.The Hall and Titterington's method consists in computing and comparing three smoothings.The local extremes of the difference between the three smoothings are used as estimates of slope change points.Moreover kernel smoothing is here replaced by spline smoothing to improve the quality of the fitting as further detailed.The proposed variation is based on a single slope sign change detection by slope comparisons.These slope comparisons are combined with smoothing comparisons to localize other slope changes between the previously detected ones.
Tests and comparisons show that the performances of the two methods are similar.Differences between touchdown with real wind and touchdown with split-linear fit are estimated with a complex simulation program.This allows us to perform free-choice but time-costly landing simulations.Implementation is carried out using S-plus (1997).
The specific problem and the data transformations are explained in Section 2. The algorithms will be presented in Section 3. The numerical results and the comparison between the methods are explored further in Section 4. We will conclude by a brief discussion in Section 5.

Aeronautical Context and Data Transformations
We aim to decompose wind speed evolution into linear components to study later their own influence on the behaviour of an aircraft autopilot.The stage of the flight we are interested in is the final approach, specifically the last 30 seconds before touchdown when wind perturbations are critical because of the imminence of the landing.Split-linear fit of wind curves will help to confirm the idea that the linear components of the wind are the main factor towards unusual landings.Then il will be possible to discriminate between severe landing due to particular external conditions and bad tuning of the autopilot.Moreover, the split-linear fit of real flight winds will add to the wind databases usually used for landing simulations.
Only the longitudinal wind component W x (t) along the runway axis is considered.Landing quality is estimated by the distance in meters between the beginning of the runway and the touchdown on the runway axis (X T P ).Our sample of data consists of 59 winds measured during flight tests.As we focus on the latest 30 seconds before landing, each wind recording is a time serie of N =240 points sampled every 1 8 th of a second.Simulation software allows us to perform landing simulations.This is convienent since the criterion to be met is expressed in term of error between touchdown performance (X T P ) instead of the usual residual square error.This software is used to set the different values of the parameters and to evaluate the algorithms.The difficulty we have to cope with, is that it cannot be included in an optimisation or a cross validation program.
Direct applications of the following algorithms on raw data are not satisfactory.Successful computations have been carried out using smoothed wind curves.We aim to determine only the significant slope change points and avoid finding too many knots.Experiments with raw data lead to non-significant point detection when using Hall and Titterington's algorithm and to difficulties in setting up the prior distribution of the number of knots in the Bayesian algorithm.Moreover, the estimated errors between touchdowns with the split-linear fits of raw data are too large.
So a necessary first step consists in clearing winds of non significant information.That is, components which do not change the touchdown.It may be measurement noise as well as some frequency components which do not disturb the aircraft during the flight.To be more specific, high frequencies have little influence on our aircraft trajectory because of its inertia.Smoothing splines are used to get rid of the high frequencies and allows the slope change detection to be efficient.Note that Fourier methods are not applicable here since the data are neither stationary nor periodic.
Let v(t i ) be the cubic spline fit of the velocity values The smoothing parameter λ governs the trade-off between smoothness and goodness of fit and is usually estimated by cross validation.In our particular case, we aim at keeping only the significant components of the wind towards the landing.Therefore goodness of fit is estimated here by the error between the real X T P touchdown with raw data and the one with smoothed data.This criterion replaces the usual residual error minimization.Empirical cross validation has led us to choose λ = 10 −5 as the largest smoothing parameter minimizing error mean and scatter (4% error maximum).

Split-linear Fitting Algorithms
The aim is to detect the slope change times t * j among the times {t 1 , t 2 , . . ., t n }.The split-linear algorithms are detailed in this section.They are applied on smoothed wind data as previously explained in Section 2. Observations are denoted (t i , w i ) whereas the values of the smoothed wind curves are denoted by v i .
Since one matching condition is to detect only significant slope changes, each method implies the search for the optimum values of the different parameters which govern the precision.

Variation on Hall and Titterington algorithm
The features of this algorithm inspired by Hall and Titterington (1992) are as follows: 1.The smoothed wind curve is split into intervals I k , k = 1, . . ., K.
2. In each interval, the slope a k of the best straight line is estimated by ordinary least-squares. 5. Since every t * j , t * * j has been detected, the split-linear fitting of the smoothed wind curves is carried out by linking the v * j = v(t * j ) or v(t * * j ).
Linear least-squares model (step 1. and 2.) The size of the K intervals governs the trade-off between accuracy of the method and calculation time.In our case, tests lead us to select 4-points intervals.In each interval, the piecewise linear least-squares model is fit to the data :

Smoothing comparison (step 4.)
The bandwidth parameter h has to be chosen for each smoothing.It determines how large is the neighbourhood of each point used to calculate the three smoothings.
For each point in [1 + h,N−h], left, right and centre smoothing are computed.Hall and Titterington propose smoothers such as particular linear combinations of the observations around each point.This smoothing has been tested and gave worse results on our data than spline smoothing.So spline smoothing has been chosen for left and right smoothing whereas centre smoothing is obtained by averaging.
Let f r (i) be the value of the right smoothing at the point i computed on For the two spline smoothings, the smoothing parameter λ has to be chosen.Moreover, the bandwidth h has to be large enough to avoid side-effects of the smoothing.But since 2h points are lost for the slope change detection, a trade-off has to be sought.
These two parameters λ and h can be optimized by cross validation.In our case, we have to choose prior satisfactory empirical values and perform the choice further as regards the numerical results in Section 4.
The diagnosis is made by estimating for each point If the point i is a slope change point then this difference will be a local maximum.
Note that since the data are smoothed before the algorithm is applied, f r (i) The detection of the local maxima is made by the comparison of each point with its two neighbours.A point i is a local maximum if : In brief, the parameters to set are : the number of intervals K, the optional T min , the smoothing parameter λ and the bandwidth h.

Bayesian automatic curve-fitting
Recall here that we aim at identifying the slope change times t * j among the t i .Therefore, the algorithm of Denison et al. is applied on our data so as to obtain a split-linear fit and only modified to obtain the final knots.The Bayesian basis of the algorithm is briefly described here according to Denison's et al. (1998) and Green's (1995) articles.
The piecewise polynomial model to fit is as follows: The dimension K of the true model is unknown as well as θ (K) = (t * 1 , . . ., t * K , σ 2 ) representing knots positions and variance of the model error.Inference about K and θ (K) is based on the joint distribution of (K, θ (K) , v).The best (K, θ (K) ) are assumed to maximize the posterior distribution p(K, θ (K) ).The coefficients β n,j in (3.1) are then estimated by least-squares minimization.
Since an analytic or numerical expression of this complicated distribution is impossible, samples of the joint posterior distribution have to be simulated.Denison et al. (1998) propose the design of three types of move (addition,deletion and movement of a knot) between the subspace Θ K for θ K which allows us to explore all the subspaces, searching for the (K, θ K ) which minimizes Each move is chosen with a specific probability which depends on the prior distribution of K.
The prior probability of the dimension K of the model is assumed to be a Poisson distribution with parameter λ p : These plots confirm that the Denison et al. (1998) algorithm gives better performance on "smooth" curves and show that the value of λ p is impossible to set for raw data.
The algorithm of Denison et al. can be summarised in three steps as follows: 1. Initialisation : choose λ p knots uniformly among the t i and set k equal to the number of interior knots.
2. Choice of the type of move : generate u uniformly and compare it to b k and d k to choose a birth, death or move transition.
3. Iteration : repeat step 2 until the MSE no longer decreases.
In brief, the parameters to be set for this algorithm are the smoothing parameter λ and the parameter λ p of the prior distribution of K.
Note that the quality of fit is measured by the MSE since, as already explained, it is impossible to replace the criterion by the one we focus on, that is difference between touchdown with the fitted wind curves and the real ones.However, verifications are carried out (see estimation of error in Section 4.) to ensure that the fits were satisfactory according to our criterion.

Practical Comparisons
In this section, both algorithms are applied on the wind data described in Section 2. The test sample is composed of all 59 available winds.Each algorithm is tested separately at first and then a numerical comparison is carried out.The variation on Hall and Titterington's algorithm is compared to the Bayesian approach.The Bayesian approach is here considered as a benchmark whose quality is well-known but computational costly.
Graphic comparisons are made between raw data and split-linear fitting to check the quality of the methods.Then touchdown are simulated with the two kinds of wind, the errors are estimated and compared so as to choose the best algorithm.
The tests of the algorithm inspired by Hall and Titterington's one are carried out with empirically best values of the different parameters to be set.To carry out a fair comparison with the Denison's one, we have to focus on the final number of knots K (slope change point) as we are searching for a parsimonious linear-fits.The comparison protocol is as follows: 1. Choose the maximum number of knots allowed K max (in our case, K max = 20 since too close slope changes do not have any physical meaning).
2. For each wind, check that the number K class < K max (K class is the number of knots obtained by the classic algorithm).
3. For each wind, set λ p (parameter of the Poisson distribution of the number of knots) such as |K class − K Bayesian | ≤ 3 with K Bayesian number of knots given by the Bayesian linear curve-fitting algorithm.
4. Carry out the graphical comparisons for each wind.
5. Carry out the numerical comparisons for each wind by simulating the X T P impact obtained with the two split-linear fits.

Test of the variation on Hall and Titterington's algorithm
The parameters are set to the best empirical ones (λ = 10 −5 and h = 15 points).Further results in Section 4.2 will show that they also are the best as regards to error minimization.
Figure 2 et 3 illustrate the interest of the adaptation performed on the Hall and Titterington's algorithm.

Test of the Bayesian automatic curve-fitting algorithm
The graphic example shown in Section 4.1 (Figure 3) is then compared to the results given by the Denison et al. algorithm on the same wind.The comparison protocol is followed.K class is estimated then λ p is chosen such that |K class − K Bayesian | ≤ 3.
For the wind shown here as an example, K class = 15, K Bayesian = 13 for λ p = 0.15.The graphical result seems to be less satisfactory than the one given by the first algorithm.But as our criterion to be met is mainly to obtain the same X T P touchdown with the split linear fit and the raw wind as explained in Section 4.2, we carry out numerical comparisons.

Estimation of error
As previously said, the main criterion of satisfaction consists in the estimation of error.The usual residual error minimization is then replaced by the minimization of the difference between touchdown with raw data and the split-linear fit.
A subset of 10 winds is extracted.A split-linear fit is estimated with each algorithm using the same values of parameters as in Section 4.1.Then touchdown simulations are computed comparing the value of X T P obtained with each raw wind and its three fits.The errors are estimated in percent of the X T P calculated with the raw wind.
The numerical comparisons are carried out in the same way as the graphic comparisons.The comparison protocol is applied here.For each wind, K class is estimated and λ p is chosen such as |K class − K Bayesian | ≤ 3.
Figure 5 depicts the boxplot of errors of the two algorithms.The performance of these two methods are similar.

Discussion
The purpose was the split-linear fitting of time series of wind velocity before touchdown.Two algorithms have been discussed in this article and compared on real-life data issued from wind measurement during flight tests.We aimed at determining all the significant slope changes in the wind curves.The criterion to meet was the achievement of the same simulated touchdown with the raw wind and its linear fit.Two kinds of approach have been considered.The algorithm based on slope or smoothing comparisons (Hall and Titterington) has been adapted and compared to a Bayesian curve-fitting algorithm used as a benchmark.This involved particularly the choice of best values for the parameters of each method.For the first algorithm, the smoothing bandwidth and parameter had to be set while the prior distribution of the number of slope changes was the main parameter to consider for the Bayesian algorithm.The main features of this article is the interest of a prior smoothing of the data and the setting protocol for all parameters.These settings have to be carried out empirically in our case because of the specificity of the satisfactory criterion to meet.The performance of these two methods are anyway satisfactory and similar.They both proove that only a piecewise approximation of the smoothed wind must be taken into account in further landing simulation processes.The choice between the two methods can only be made according to the difficulty to arrange the parameters and time calculation.
3. The consecutive slopes a k and a k+1 are compared.If a k .ak+1 < 0, a point of slope sign change t * j is detected.4. If the interval [t * j , t * j+1 ] is larger than T min , points of single slope change t * * j are sought by the comparison of three smoothings further detailed.

Figure 1 :
Figure 1: This Figure depicts the influence of λ p on the final number of knots K obtained for two raw and smoothed winds.The difficulty to set λ p for raw data leads us to use prior spline smoothing on our data.

Figure 2 :Figure 3 :
Figure 2: Example of an unsatisfactory graphic comparison between the raw data and the split-linear fitting obtained with the Hall and Titterington's algorithm.At least five slope changes are forgotten at times 15, 190 and between times 210 and 250.

Figure 4 :
Figure 4: The raw data and the split-linear fit obtained on smoothed data with the Bayesian algorithm are plotted.This Figure is associated with Figure 3.

Figure 5 :
Figure5: Boxplots of errors of the two algorithms estimated in percents of the X T P calculated with the raw wind.The results are similar.In the two cases, the maximum error is less than 4.5%.