domain: phys stat spi

text:
Alireza Keramat
email: keramat@ust.hk

Mohamed S Ghidaoui

Xun Wang

INVERSE TRANSIENT ANALYSIS FOR PIPELINE LEAK DETECTION IN A NOISY ENVIRONMENT

Keywords: Leak detection, inverse problem, optimization, adjoint equations, waterhammer. 1

In the transient-based leak detection methods in the time domain, the solution parameter set is determined by minimizing an objective function that fits the numerically modeled pressures to measured pressures. If the measurement is contaminated with noise, then the optimum decision parameters may be some fictitious leaks whose computed pressure head is best matched to the noise not the reflections from leaks. This issue is explained in detail and then a different merit function to decline it, is suggested. It is based on maximizing the signal-to-noise ratio (SNR) which in literature usually termed as Matched Field Processing. This function is then compared with the conventional least square objective function and its efficiency subject to noise is assessed. It was found that the efficiency of each method mostly depends on the variance of noise and the system properties which may change from case to case. More specifically, the proposed merit function is more robust when its optimal variance is smaller than that of the least square cost function.

INTRODUCTION

Intentional excitation of fluid flow in pipe networks and collection of transient pressure data at some points is widely established as a tool for leak detection in fluid-filled pipes [START_REF] Colombo | A selective literature review of transient-based leak detection methods[END_REF]. Many methodologies and tools, such as leak reflection method, inverse transient analysis (ITA), impulse response analysis, transient damping method and frequency domain response analysis have been developed for solving various types of leak detection problems in recent years [START_REF] Colombo | A selective literature review of transient-based leak detection methods[END_REF]. ITA aims at finding properties of the pipe system by means of a transient flow data in the system. In the implementation, the method assumes potential leaks with unknown discharge at computational sections and then minimizes the summation of absolute differences between the measured data and the computational results. First introduced by [START_REF] Liggett | Inverse transient analysis in pipe networks[END_REF], a MOC-based transient solver in conjunction with Levenberge-Marquardt optimization algorithm is employed to minimize the misfits. The methodology proposed in [START_REF] Liggett | Inverse transient analysis in pipe networks[END_REF] has been refined by numerous researchers over the years with particular attention focused on methods for optimization (e.g. [START_REF] Vítkovský | Leak detection and calibration using transients and genetic algorithms[END_REF][START_REF] Kapelan | A hybrid inverse transient model for leakage detection and roughness calibration in pipe networks[END_REF]Jung and Karney 2008;[START_REF] Haghighi | Detection of Leakage Freshwater and Friction Factor Calibration in Drinking Networks Using Central Force Optimization[END_REF] and the hydraulic transient model (e.g. [START_REF] Vítkovský | Inverse analysis and modelling of unsteady pipe flow: Theory, applications, and experimental verification[END_REF][START_REF] Brunone | Detecting leaks in pressurized pipes by means of transients[END_REF][START_REF] Covas | Inverse transient analysis for leak detection and calibration of water pipe systems modelling special dynamic effects[END_REF]. Indeed, the so-called inverse transient analysis demands an accurate mathematical and numerical model which can favorably represent the system behavior. Otherwise, the residuals of the fit are wrongly estimated which finally results in fictitious optimum quantities. To this aim, many works have been conducted to improve the water hammer models so as to decrease modelling errors.

Studies on optimization aspects of the ITA often aim at suggesting efficient algorithms to evaluate gradients [START_REF] Nash | Efficient inverse transient analysis in series pipe systems[END_REF][START_REF] Vítkovský | Inverse analysis and modelling of unsteady pipe flow: Theory, applications, and experimental verification[END_REF]. There are also significant records to address optimal measurement site locations [START_REF] Vítkovský | Optimal measurement site, locations for inverse transient analysis in pipe networks[END_REF] or minimization algorithms [START_REF] Kapelan | A hybrid inverse transient model for leakage detection and roughness calibration in pipe networks[END_REF][START_REF] Haghighi | Detection of Leakage Freshwater and Friction Factor Calibration in Drinking Networks Using Central Force Optimization[END_REF]. Despite that numerous researches have been carried out to improve the ITA for leak detection, still there are several drawbacks which demand more investigations. In fact, any inverse solution should be tested for the presence of data error such as random noise in measurements and uncertain input parameters. Studies dedicated to resolve this reason of possible failure of ITA for leak detection are indeed scanty.

In the leak detection based on ITA, usually there are too many decision parameters because all computational points are treated as potential leaks in the model. A large number of design parameters promote high correlation between the simulated pressures as well as the optimum parameters. As a result, in a noisy environment and when the optimization parameters are too many, the outputs of optimization are no longer unique thus leading to wrong localization. An effective solution to this problem is to limit the number of unknown leak candidates, based on which a penalty function is introduced in the merit function. It is known as model parsimony studied by Vítkovský et al. (2007). Another remedy should be applied to the collected signal to limit the noise effects in the localization. This issue is the purpose of this paper. Based on maximizing the energy of the signal to the energy of the noise, signal-to-noise ratio (SNR), a merit function different than the conventional least square function is applied for the optimization in ITA. In the signal processing literature, this approach is usually termed as matched filter or matched field processing [START_REF] Krim | Two decades of array signal processing research[END_REF]. It is widely used in acoustics for sound source localization [START_REF] Wang | Estimation of multiple sound sources with data and model uncertainties using the EM and evidential EM algorithms[END_REF][START_REF] Tolstoy | Detecting pipe changes via acoustic matched field processing[END_REF]), structural vibration (Turek and Kuperman, 1997;Tippmann and Lanza di Scalea, 2015) and ocean engineering [START_REF] Baggeroer | An overview of matched field methods in ocean acoustics[END_REF]) by matching the modeled and measured signals or blockage detection in sewer pipes by correlations between measurements [START_REF] Tolstoy | Waveguide monitoring (such as sewer pipes or ocean zones) via matched field processing[END_REF].

In this research, the merit function based on matched field processing is incorporated for localization in ITA. It is then compared with the least square norm to investigate its efficiency in detections when the measured signal is contaminated with noise. It is found that the two merit functions arrive at different results when the noise level is high (very low SNR). A criterion regarding the priority of each one for detection based on their optimal variance is provided and discussed. It was concluded that for high noise levels, an improved methodology regarding the choice of the two cost functions can lead to more reliable localizations.

MOTIVATION

In the time-domain based inverse transient analysis, all computational points are treated as potential leaks. An optimization tool is then applied to solve for effective leak areas. Traditionally, the least square fit between measurements and computed heads is utilized as the merit function to find the optimum leak sizes.

A noise-free signal contains a substantial amount of information which are manifested by the wave reflections from pipe system boundaries including anomalies. In the presence of noise, these reflections are contaminated so that the results of the model are not precisely matched to the measured data. It means that the misfit quantity coincides to the estimated variance of noise for the actual leak parameters. However, some issues can lead the detection procedure to wrong leak parameters and overfitting to measurements. In other words, the global minimum is displaced and it is no longer at the actual leak parameters or at least it is not the unique global extremum. These issues are basically driven by the misfit function variables being either the computed head or measure data. Adding too many parameters to the ITA promotes high correlations between computed signals and makes the inverse problem non-unique. This issue which stems from assuming all nodes as potential leaks (too many parameters) in ITA is more pronounced when the measurement contains noise.

The problem is more illustrated by Fig. 1 while a noisy signal (white noise) with SNR = 0 dB in Fig. 1 (a) is applied for leak localization. Typical result from an ITA are presented in Fig. 1 (b) in which the actual and predicted leak size and locations are represented by red and blue bars respectively. Fig. 1 (c) shows the optimization iterations when the starting point (initial guess) is the actual leak parameters. As seen, a set of wrong leak locations and sizes are detected. In this case, the averaged misfit square is less than the noise variance. In fact the method tries to find fictitious leaks to match to the white noise. As a result, the global minimum of the least square fit between measured and computed signal is displaced and it corresponds to another set of leaks which differ with the actual one both in location and size. Note that this is an example with no modeling error with only leak parameters (effective leak area at each node) as decision variables of optimization. When too many decision variables are incorporated (e.g. friction factor, viscoelastic and unsteady friction coefficients, etc.) the noise effect in the localization becomes more problematic and thus the detection gets remarkably inaccurate. The aim of this analysis is to find the properties of the pipe system, particularly leak parameters herein. The leak areas at each node are considered as decision variables of an objective function to be minimized subject to water hammer equations. In fact this leakage detection method is based on several forward and adjoint water hammer analyses which are directed by an optimization scheme. Firstly, the forward leakage problem in case of plain waterhammer waves is explained. Next, the optimization procedure and a new merit function for leak detection in a noisy environment is discussed. Finally, the approach of the adjoint problem is employed in the optimization procedure to efficiently find the gradients of the objective function with respect to the desired parameters.

Mathematical model and numerical solution

Two equations of mass and momentum conservation govern transient flows in pressurized pipes. The classical water hammer theory which assumes that pressure is constant over the flow cross section is used herein. Additionally, an averaged velocity is applied in the momentum and continuity equations, allowing for a one-dimensional representation of water hammer equations. The Method of Characteristics (MOC) is used to solve the governing equations. Using this method, the two equations are transformed to two algebraic equations which are valid on characteristic lines. These equations in conjunction with the leakage orifice relation are solved simultaneously for flow rate at either sides of each node as well as the pressure head. All significant effects including unsteady friction and viscoelasticity can efficiently be incorporated into the resulting characteristics equations [START_REF] Keramat | Waterhammer with column separation, fluid-structure interaction and unsteady friction in a viscoelastic pipe[END_REF]. This allows for the incorporation of a precise water hammer solver which in turn leads to reducing the modelling errors in the detection procedure. During transient flow, leaks are modelled using the orifice formula. The effective leak coefficient (Ae = Cd AL) is assumed to remain as that of steady state, though this is a controversial assumption and could be the source of significant flaw in the simulation. The orifice relation gives

1 1 1 12 ,2 n n n i i i i ei i Q Q Q H A g H         (1)
where Qi is the flow rate of leak i (see Fig. 2) and Q1i and Q2i stand for its upstream and downstream discharge respectively; superscript n stands for time steps and no-superscript indicates steady state.

Figure 2. Schematic of a reservoir-pipe-valve system with multi-leaks to be used for ITA.

Optimization

The nonlinear constrained optimization using a build-in MATLAB function called "fmincon" is carried out. Among several algorithms of this function, the Sequential Quadratic Programming (SQP) is employed. This method was also applied by [START_REF] Shamloo | Leak detection in pipelines by inverse backward transient analysis[END_REF] for the leak detection. It is demonstrated to be efficient in constrained optimization with nonlinear objective and constraint functions. SQP methods solve a sequence of optimization sub-problems, each of which optimizes a quadratic model. The solution involves major and minor iterations. Major iterations generate a sequence of iterations that satisfy linear constraints. Iterations converge to a point satisfying the first-order conditions for optimality. At each iteration, a Quadratic Programming (QP) subsystem is used to generate a search direction toward the next iteration. Solution of the QP subsystem is itself an iterative procedure, with minor iterations of a SQP method being iterations of the QP method [START_REF] Barcelay | SQP methods for large-scale optimization[END_REF]. This procedure is performed by existing MATLAB function. The required gradient vector of the objective functions is calculated using the adjoint solution of waterhammer equations and are fed to the optimizing solver.

Maximum Signal to Noise ratio

In the inverse transient analysis, the pipe parameters (leak area, friction factors, wave speeds, etc.) are fitted in the way that the model outputs h(t) become as close as possible to the measured pressures h m (t). The observed pressures are usually far from the corresponding model simulations due to either data error such as poor measurement instrumentations and calibration, or modelling errors such as inconsistent boundary conditions or incomplete mathematical modelling of the physical phenomena. As a result, the inverse transient analysis

Qi Q1 Q2

Reservoir Valve cannot precisely determine the system specifications unless appropriate mathematical or statistical operations are performed in conjunction with each nature of discrepancies.

The following objective function minimizes the difference between the model results and the measurements:

      2 1 0 T m E h t h t dt    (2)
that is in fact the integration of misfit quantities over all time sections in 0 tT  , where T represents the last collected time. This least square criterion is derived from maximum likelihood estimators. Another merit function can be defined using the measured and computed heads. It can be achieved by finding a filtering function based on computed heads such that it maximizes the energy of the signal. This filter in fact adjusts the measured heads by increasing the reflections from actual leaks and decreasing noise energy. The convolution integral of a unit weight function   Wt with the measured time series   m ht at the last time point T is calculated to this aim:

      * * * m W h T W h T W n T   (3)
where n(t) represents a zero mean white noise and W is assumed to be a unit vector, i.e. 1 W  , to obtain a non-trivial solution. To find the most efficient filtering function, the energy of the signal divided by the expected value of the noise energy, which is defined by signal-to-noise ratio (SNR), should be maximized. Considering Schwarz's inequality, the energy of the signal at each time t is upper bounded by

                2 2 2 2 0 0 0 * .

T T T W h T W T s h s ds W T s ds h s ds

                         (4)
Noise is random and assumed to be independent at different time steps so the denominator of SNR (energy of the noise) is simplified to

            2 2 2 2 0 0 * , T T W n T W T s n s ds W T s ds                                E E (5)
in which  stands for the standard deviation of the noise. As a result, the SNR reduces to

                        2 2 2 2 0 0 2 2 2 0 2 0 * 1 , * t t t t

W t s ds h s ds W h t h s ds W n t W t s ds

                                      SNR = E (6)
In this inequality as well as (4), the function h(s) is deterministic for the specified system properties (leaks specifically), so equality condition corresponds to the maximum value of the energy of the filtered signal.

            2 2 2 2 2 2 0 0 1 1 * T T m m m E W h T W T s h s ds h t h t dt                   h h (7)
This objective function should be maximized to find the best leak parameters which correspond to the highest filtered energy of the measured signal to the noise. It is called matched filter or Matched Field Processing (MFP) in the signal processing literature [START_REF] Krim | Two decades of array signal processing research[END_REF].

Discuss merit functions

Both objective functions (2) and ( 7) in a noise-free condition arrive at a known global extremum being zero and 2 h , respectively. In a noisy environment, the global extremum may be displaced thus leading to wrong localization. To study the robustness of the two merit functions, a statistical analysis is performed. The variance of the two random functions (2) and ( 7) is estimated to this aim. Assuming that the mean value of measured head is equal to the computed head, the variance of Eq. (2) becomes

1 2 2 4 2 m E h N    (8)
in which N is the number of time sections of computed head and m h  is the variance of measured head. For Eq. ( 7) it yields

2 2 2 4 2 2 4 m m E h h     h (9)
Basically, a smaller variance for the cost function corresponds to a better localization because it is then less sensitive to the random noise. As seen, both variances depend on numerical simulations (due to N and h ) and the noise properties. For a fair comparison between the two, these values should first be found and then appropriate cost function is chosen. In Eq. ( 9), the norm of the computed head corresponding to the actual leaks 2 h cannot be found in real experiments because it depends on leak parameters. As a result, general conclusion regarding the two merit functions cannot be made. However, it is demonstrated in this study that when the quantity of ( 9) is smaller than (8), on average, the detection based on matched field processing is more robust, and vice versa.

Adjoint method

It corresponds to define a different problem in which the misfit function has the same gradients with respect to the design parameters as the forward problem. It can be favourably illustrated in its discrete form where the governing equations are discretised and set to a matrix-form representation to solve for the state variables. Therefore, the simultaneous system of equations over the entire spatial and transient period is Mx = R (10) in which M, R and x are respectively coefficient matrix, right-hand-side vector and vector of state variables containing the unknown heads and flow rates at all time and space sections. The procedure of generating the matrix relation ( 10) from the characteristics equations is explained by [START_REF] Liggett | Inverse transient analysis in pipe networks[END_REF]. This matrixrepresentation of the compatibility equations can also be prepared at each time step, but for all computational nodes making analogous simultaneous equations to Eq. ( 10). This approach, termed as implicit method of characteristics, also explained by [START_REF] Afshar | Water hammer simulation by implicit method of characteristic[END_REF], [START_REF] Rohani | Simulation of transient flow caused by pump failure: Point-Implicit Method of Characteristics[END_REF] and most recently, by [START_REF] Wang | Water Hammer Simulation Using Explicit-Implicit Coupling Methods[END_REF]. As in general, the equations are nonlinear, this matrix relation (10) shows the linearized form so matrix M and vector R contain terms of state variables which are updated through iterations. As a result, the transient solver iterates (three to five times) to converge to the solution. The derivative of the objective function with respect to the decision parameters j b (lumped leak coefficients or friction factors) is called gradient vector. Each elements of this vector is represented by

d . d T j j j E E E b b b               x x (11)
The direct calculation of / j b  x by the divided difference method is inefficient because it requires calling the solver for each j b derivative computation. Eq. ( 10) can be used to find each column of the Jacobian matrix   /  xb analytically, as follows: 11) results in:

1 . j j j j j j b b b b b b                         M x R x R M x M M x (12) Substitution for / j b  x in Eq. (
1 d d T j j j j E E E b b b b                         RM M x x (13)
The desired adjoint matrix equation ( 14) is made with regard to the coefficient of the second bracket in Eq. ( 13):

1 T T T E E            M ψ M ψ x x (14)
where ψ contains the adjont variables to state variables x and remains identical when calculating all elements of the gradient vector. It therefore needs to be calculated only once for all bj parameters. Having found the adjoint vector and considering Eq. ( 13), each elements of the gradient vector is given from Eq. ( 15):

d . d T j j j j E E b b b b                RM ψ x ( 15
)
Another approach based on Lagrange multiplier to find this adjoint equations can be found which leads to the same relations. To find the Hessian matrix, the Jacobian matrix is needed. The elements of this matrix can be found using the same procedure as presented above but now, the objective function E is defined as one variable at one specific time point (Vitkovsky 2001;[START_REF] Liggett | Inverse transient analysis in pipe networks[END_REF].

NUMERICAL RESULTS

A reservoir-pipe-valve system is considered to investigate the leak detection method using the explained inverse transient analysis. In order to verify the primary results, the system characteristics of the experimental apparatus in the Robin Hydraulics Laboratory in the School of Civil and Environmental Engineering at the University of Adelaide [START_REF] Bergant | Water hammer and column separation measurements in an experimental apparatus[END_REF]) are adopted to run the numerical tests. Details of this system are provided in Table 1. Considering Eqs. ( 8) and ( 9), two leak cases according to Table 2 are defined to see the robustness of the two merit functions in each case. The two cases are discriminated by their variance given by Eqs. ( 8) and ( 9). In the leak case 1, merit function E1 has lower variance and for the leak case 2, cost function E2 has less variance.

General specification of the reservoir-pipe-valve system Details of leaks used for the ITA and corresponding leak-induced pressure head reflection and variances of each merit function.

Cd AL (m 2 ) 1 2 2 4 4 2 (m ), SNR=-5,0,5 m E h N    2 2 2 4 2 4 2 4 (m ), SNR=-5,0,5 m m E h h     h Leak case 1 7 2 10   2×10 4 , 2×10 3 , 2×10 2 1×10 5 , 4×10 4 , 1×10 4
Leak case 2 6 3 10   5×10 8 , 5×10 7 , 5×10 6 8×10 6 , 2×10 6 , 8×10 5

Based on the data provided in Tables 1 and2, the forward problem is solved to find the pressure head at the valve h(t). A Gaussian white noise with zero-mean unit variance represented by a random variable n(t) was used to hypothetically generate the "measured" data indicated by superscript "m"

      -SNR 20 , 10 mm m L hh h t h t n t h      (16)
In this relation, SNR (dB) is a given number which depends on the system instrumentation and collected experimental data and  is the variance of the noise. One can deduce that this pressure reduction is almost proportional to the leak size using the provided analytical formula by [START_REF] Vítkovský | Inverse analysis and modelling of unsteady pipe flow: Theory, applications, and experimental verification[END_REF]. This way of numerically generating noisy pressure signal allows for having higher noise corresponding to higher leakage size, therefore it roughly makes the noise-level investigations independent from the leakage size.

In order to estimate the efficiency of each merit function subject to noise, an error estimation criteria may be defined:

1 ˆLL XX L    (17)
in which L is the pipe length, L X is the leak location, ˆL X is the estimated leak location which corresponds to the node with maximum effective area obtained from optimization. This criterion does not take the detected leakage size into account because usually it is less critical for clients. Another criterion as follows may be used which takes both leak size and location into account in order to theoretically compare the two merit functions.

 

1 2 ˆN e L ei Li i eL A X A X AX      (18)
In this relation, Ae is the actual leak size and Aei, i =1 ... N, is the detected efficient area for each computational point and N is the number of MOC nodes.

Leak case 1

Considering the problem stated in section 2, the proposed and the conventional merit functions are compared. To this aim, a number of 0-SNR signals are used for localization with the two cost functions in ITA. To make a general conclusion, the detection procedure was repeated for 10 different measurements at the valve contaminated by white noise. The result is then averaged and presented in Fig. 3 for the two merit functions. The results of this figure show that the two cost function behave almost similar. In fact according to Eq. ( 18 To further investigate the two cost functions, the detection is repeated for 60 times and the error of each run is evaluated according to Eq. ( 17). Fig. 5 shows the results for four different SNR. It demonstrates that the two merit functions differ only at low values of SNR and they behave in the same manner for SNR higher than 5. These results are in agreement with the theoretical conclusion which was made based on the variance of the two cost functions. Since the least square function have lower variance, the localization using this cost function is better off. 17) for 60 simulations with different SNR.

Leak case 2

A similar investigation to the previous section but now with another leak is made. Fig. 5 is the averaged results of 10 runs for 0 dB signals and Fig. 6 is the averaged results of for -5 dB noisy signals. Errors evaluated using Eq. ( 18) for the matched field and least square merit functions, respectively are: for SNR=-5 dB (Fig. 6). The comparison reveals that for this leak case, the matched field cost function is definitely preferred since its ptimal variance is smaller than that of the least square function (Table 2). The error analysis based on Eq. ( 17) is carried out for 60 simulations in order to compare the two metrics in leak localization. The mean and 95% confidence interval of error for each SNR is presented in Fig. 7. The conclusion regarding the excellence of the proposed merit function for very low SNR is confirmed. For the reason that the variance of the matched field merit function is lower, the average error in localization is less and it is more robust than the least square function.

Figure 7. The mean and 95% confidence interval of localization error evaluated using Eq. ( 17) for 60 simulations with different SNR.

CONCLUSIONS

The leak detection based on the ITA in the time domain when measured data are contaminated with a white noise is addressed. The problem of non-uniqueness and overfitting which stems from having too many decision variables (MOC nodes as leak candidates) arises when measured signal is noisy. It is more significant under high levels of noise or equivalently in systems with very small leaks which have reflection amplitudes of the same order or smaller than noise standard deviation. The analysis based on maximizing signal-to-noise ratio which is largely stablished as a successful signal processing technique is carried out. The idea is equivalent to the matched filter or Matched Field Processing (MFP) which involve correlations between measured and modeled signals. It was led to a new merit function for the ITA leak detection method instead of the conventional least square criterion. The two possible merit functions are proved to have the same global extremum for localizations using noise-free measurements. The least square function has zero global minimum and the matched field function has global maximum equal to L 2 norm of the signal. These two extremum occur when the modeled and noise-free measured signals coincide. In a noisy environment when the non-uniqueness problem in the ITA becomes a significant issue, the two objective functions behave differently. It was found in this research that the optimal variances of the two cost functions are the distinguishing parameters to decide for the selection of each function for leakage detection.

The merit function with lower variance was more robust for localization. This finding is closely investigated by calculating the variances of the two merit functions for two case problems. The numerical examples are solved and on the basis of a detection-error criterion, the preference of each one to the other is demonstrated.

The incorporation of the matched field merit function in the ITA can lead to an improved localization procedure in a noisy environment. The two merit functions are used for localization with a given measured signal (probably the average of several measurements at a location). If then they arrive at different solution, that of the cost function with lower optimal variance is supposed to be more reliable. Note that a classical solution for the raised non-uniqueness or overfitting problem is to use model parsimony, which favors a lower number of leak candidates. The proposed matched field cost function can still be used with model parsimony. The two cost functions can be applied to enumerate the domain and achieve a more reliable result. Under high noise levels, they may reach to different solutions. Discussions regarding the preferred localization can be valid herein.

Figure 1 .

 1 Figure 1. (a): noisy pressure head with 0-SNR (blue) and best computed fit (red). (b) the actual (red) and detected (blue) leak parameters (location and size). (c): The convergence curve of optimization iterations.

  the first reflection in pressure due to the leak (a step pressure drop) and m h

  equal to -5 but now the average of 30 runs are shown in Fig. 4. Then the evaluated errors according to Eq. (18it demonstrates similar localization though with slightly different sizes.

Figure 3 .

 3 Figure 3. The actual (red) and detected (blue) leak parameters (location and size) for the two merit functions. The results are the average of 10 measurements, each with SNR=0 dB.

Figure 4 .

 4 Figure 4. The actual (red) and averaged detected (blue) leak parameters (location and size) for the two merit functions. The results are the average of 30 measurements, each with SNR=-5 dB.

Figure 4 .

 4 Figure4. The mean and 95% confidence interval of localization error evaluated using Eq. (17) for 60 simulations with different SNR.

Figure 5 .

 5 Figure 5. The actual (red) and detected (blue) leak parameters (location and size) for the two merit functions. The results are the average of 10 measurements, each with SNR=0 dB.

Figure 6 .

 6 Figure6. The actual (red) and averaged detected (blue) leak parameters (location and size) for the two merit functions. The results are the average of 30 measurements, each with SNR=-5 dB.


domain: spi stat nlin

text:
Fedi Zouari
email: zouari.fedi@gmail.com

Xun Wang
email: xunwang00@gmail.com

Moez Louati
email: mzlouati@gmail.com

Mohamed S Ghidaoui
email: ghidaoui@ust.hk

SINGLE EXTENDED BLOCKAGE IDENTIFICATION USING A MODEL-BASED MATCHED-FIELD PROCESSING APPROACH

Keywords: blockage detection, transient flow, frequency response, matched-field processing, inverse problem

A new transient-based approach for single blockage detection is proposed which requires two measurement stations at both ends of a single pipe system. The method uses the frequency response and gives satisfactory results using low frequency bandwidth. The first advantage of the proposed technique is to obtain the location and length of the blockage independently from the blockage's area. The second advantage is that the method uses information not only from the eigen-frequencies, but also from the whole frequency spectrum. This brings more accuracy especially when the identification of the resonant frequencies and their corresponding mode number are not accurate. Moreover, the method has the ability of tolerating noise as long as it is Gaussian of zero-mean and has a well-known structure. Comparison with other blockage detection methods is given and the benefits and drawbacks of the proposed method are discussed.

INTRODUCTION

Extended blockages in water supply systems are ubiquitous and they cause many inefficiencies such as energy loss and increase in potential contamination and leaks. Unlike other defects in pipe, blockages are very difficult to detect as they do not have any external visible trace nor do they induce a detectable sound like leakages under steady state condition. However, researches showed that anomalies such as blockages can be detected if a transient response of a pipe system is measured (e.g. [START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF][START_REF] Meniconi | Experimental investigation of coupled frequency and time-domain transient test-based techniques for partial blockage detection in pipelines[END_REF][START_REF] Louati | In-depth study of plane wave-blockage interaction and analysis of high frequency waves behaviour in water-filled pipe systems[END_REF], Ignarcio et al., 2016).

Similar inverse problems were addressed in different research fields, ranging from the identification of blockage in a sewage line [START_REF] Tolstoy | Waveguide monitoring (such as sewer pipes or ocean zones) via matched field processing[END_REF], and in the cooling system of nuclear reactors (Wu and Fricke, 1990) to the vocal tract shape reconstruction (Shroeder, 1967;[START_REF] Gopinath | Determination of the shape of the human vocal tract from acoustical measurement[END_REF]Gopinath andSondhi, 1970, DeSalis and[START_REF] Desalis | The development of a rapid single spectrum method for determining the blockage characteristics of a finite length duct[END_REF] and the identification of conductivity in electrical lines [START_REF] Bruckstein | Inverse scattering for discrete transmission-line models[END_REF].

In general, the methods for blockage detection or area reconstruction could be divided into two main categories: methods that model the blockage as a step function, and therefore the blockage is simply a pipe section with smaller diameter (see sketch in Fig. [START_REF] Louati | In-depth study of plane wave-blockage interaction and analysis of high frequency waves behaviour in water-filled pipe systems[END_REF]). Other methods model the blockage with a smooth function of the longitudinal coordinate . In both categories, the techniques/tools used for blockage detection are differentiated by whether the time domain or frequency domain signals are used. For example, in time domain wavelet analysis is applied by [START_REF] Meniconi | Experimental investigation of coupled frequency and time-domain transient test-based techniques for partial blockage detection in pipelines[END_REF] and reconstructive Method of Characteristics (or Layer pealing) by [START_REF] Gong | Detection of localized deterioration distributed along single pipelines by reconstructive MOC analysis[END_REF]. In frequency domain, Fourier decomposition and Ehrenfest theorem can be applied to write an explicit relation between the area of the pipe and its acoustic properties (e.g. [START_REF] Schroeder | Determination of the Geometry of the Human Vocal Tract by Acoustic Measurements[END_REF][START_REF] Qunli | Determination of blocking locations and crosssectional area in a duct by eigenfrequency shifts[END_REF][START_REF] Desalis | The development of a rapid single spectrum method for determining the blockage characteristics of a finite length duct[END_REF]; and recently, [START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF] and [START_REF] Louati | In-depth study of plane wave-blockage interaction and analysis of high frequency waves behaviour in water-filled pipe systems[END_REF] combined the frequency response with the dispersion relation to determine the blockage location, length and size. So far, the frequency domain methods [START_REF] Schroeder | Determination of the Geometry of the Human Vocal Tract by Acoustic Measurements[END_REF][START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF] and [START_REF] Louati | In-depth study of plane wave-blockage interaction and analysis of high frequency waves behaviour in water-filled pipe systems[END_REF]) rely on the accurate identification of several resonant frequencies and their corresponding mode number. In many practical situations, a limited number of resonant frequencies are available and their identification may be affected by the presence of noise.

In this work a new model-based matched-field processing approach (MFP) is proposed for the detection of a single blockage in pipes. It relies on measuring the frequency response at both ends of a pipe and makes use of the transfer matrix method. In contrast to other methods in the literature, the proposed technique does not require the identification of the resonant and the anti-resonant frequencies and also tolerates the presence of noise. This approach is recently applied for leak detection in pipes by [START_REF] Wang | Pipeline leak detection using the matched-field processing method[END_REF].

PROPOSED NEW METHOD FOR BLOCKAGE DETECTION

2.1 Theoretical framework In most water supply application the location, the length and the size of the blockage are the most critical information. Therefore an idealized single blockage in a reservoir-pipe-valve system is considered as shown in Fig. [START_REF] Louati | In-depth study of plane wave-blockage interaction and analysis of high frequency waves behaviour in water-filled pipe systems[END_REF]. The blockage has a length and an area and it is located at a distance from the reservoir. is the area of the intact pipe (before the formation of the blockage). Two measurement stations are setup at E-proceedings of the 37 th IAHR World Congress August 13 -18, 2017 -Kuala Lumpur, Malaysia both ends of the pipe at distance and from the reservoir. The pipe flow is assumed frictionless and the wave speed is assumed constant. In this case, for any and , the pressure head variations from the mean in the frequency domain at distance and from the reservoir ( and ) are given from the Transfer Matrix method (Chaudhry, 2013) [START_REF] Louati | In-depth study of plane wave-blockage interaction and analysis of high frequency waves behaviour in water-filled pipe systems[END_REF] [2] where the superscript "c" in denotes computation; is the characteristic impedance of the intact conduit with is the gravitational constant; is the wave number; is a known discharge variation at the valve used to generate the transient. From Eq. ( 1) and (2) it can be shown that [3] [4] With regards to blockage detection, the quantities on the left hand side (LHS) of Eqs (3) and (4) are related to the known "measured" variables; whereas the quantities on the right hand side (RHS) are related to the unknown variables. Hence, the quantities on the LHS of Eqs (3) and (4) will be referred as the "measured" synthetic signals and .

[5]

The main advantage of using these synthetic signals and rather than directly using the measured pressure head and is that these signals ( ) can be written in a "quasi-linear" relation with respect to the change in area (i.e. when and are fixed). This will allow the length and the location of the blockage to be estimated independently from its size. That is the three dimensional problem reduces to a two dimensional one.

In practice, the measured pressure head will not be exactly equal to the predicted pressure head from the model (Eqs. ( 1) and ( 2)), it will rather contain noise such as environmental noise, measurement noise and many other types of noise. Therefore, the "measured" synthetic signal are given by .

Model-based matched-field processing approach (MFP)

Given the measurements and , for a given set of frequencies , Eqs (3) and ( 4) can be written in a compact form as [6] where ; ; in which the superscript "T" denotes transpose operator; and are 2-by-N matrices where the elements of their row are ; ;

;

; and is the added noise vector.

If the noise is zero-mean Gaussian white noise, then the blockage parameters can be found by minimizing the -norm of the difference between the "measured" signal and the model

[7]
If the noise is non-white Gaussian of a well-known structure (known covariance matrix) then a filter can be applied to the signals and to the model to whiten the noise, and thereafter the same method can be used [START_REF] Wang | Pipeline leak detection using the matched-field processing method[END_REF]. Because of the quasi-linearity between the area of the blockage and the "measured" signal, an estimate of which has the minimum mean square error (least square solution) is given by (for any and )

[8]

where the superscript "H" denotes the conjugate transpose. Replacing into Eq. ( 7) leads to [9] where . Using this approach, the number of variables in the optimization problem in Eq. ( 7) has reduced from three variables to two variables in Eq. ( 9). In this case the objective function can be evaluated and plotted for all possible and , and an estimate of the blockage length and location can be easily found by enumeration. This is one of the advantages of the proposed method comparing to other methods in the literature where the blockage identification involves an optimization problem with at least three variables [START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF] and Ignarcio et al., (2016)).

Once the blockage's length and location are estimated, the blockage area can be easily found by substituting the values of the estimated and into Eq. ( 8). Note that the objective function has a symmetry line related to the measurement location given by . Whereas has a fixed symmetry line independent from the measurement location given by . Therefore to avoid multiple solutions of and the measurement location from the valve side should not be at the valve (i.e.

). Also note that the objective functions and are not defined at the trivial solutions , , and . The procedure of the blockage identification using this approach can be summarized in the following algorithm Algorithm 1: Blockage identification using MFP 1: Use and at known locations and and to construct the signals and (Eq. ( 5)) 2: Construct the matrices corresponding to each signal 3: Compute the objective functions at every possible blockage location and length ( ) 4: Find an estimate ( and ) which correspond to the maximum of the objective function 5: Compute two estimates of by substituting ( and ) found into Eq. ( 8) using and 6: Compute the estimates of by and by

RESULTS AND DISCUSSION

To validate the proposed method, and to illustrate its advantages and limitations, a numerical case study is considered. The numerical setup consists of a reservoir pipe valve system containing a single blockage as shown in Fig. [START_REF] Louati | In-depth study of plane wave-blockage interaction and analysis of high frequency waves behaviour in water-filled pipe systems[END_REF]. The lengths of different pipe sections are m, m, m, m, m. The intact pipe area is and the change of area due to the blockage is m 2 . The wave speed is assumed to be constant m/s and the pipe is assumed to be frictionless. The pressure head response and at the measurement locations ( and ) are computed using the transfer matrix method (i.e. using Eq. ( 1) and ( 2)), where the transient is generated by an impulse (i.e.

). The frequency domain is sampled from to with a step size , where is the fundamental frequency of the intact system . Algorithm ( 1) is applied to compute the objective function at every point in the feasible domain discretized with m. Fig.

(2) shows a color plot of the objective function together with the location of the maximum. The results shows that the objective function has a single maximum which corresponds to the blockage's location and length estimates ( m and m). This is how the method provides the blockage length and location. Regarding, the blockage area, two estimates of the change in pipe area are obtained from Eq. ( 8) for each signal and . The estimates of are m 2 . Note that one of the estimated value is negative which may represent a pipe wall thinning. This negative value is ignored since it is assumed that there is only positive change in area (blockage). Hence the estimate of the area is m 2 .

Figure 2.

Objective function and location of the maximum. The above numerical example illustrates how the method can be applied to identify a single blockage in a pipe. This proposed method differs from most of the frequency domain methods for blockage detection (DeSalis and Oldham 2001; [START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF][START_REF] Louati | In-depth study of plane wave-blockage interaction and analysis of high frequency waves behaviour in water-filled pipe systems[END_REF] by the fact that it uses the frequency response of any set of frequencies rather than using only the resonant and the anti-resonant frequencies locations. Using the frequency response at all frequencies gives more accuracy in the blockage identification especially in the cases (i) presence of noise (ii) inaccurate identification of the eigen-frequencies. These two advantages will be discussed in what follows.

(

i) Presence of noise

To illustrate the advantage of using the response of frequencies other than the resonant and anti-resonant frequencies, the results for the blockage identification are compared in the presence of noise. A zero-mean Gaussian white noise (the noise power is the same at all frequencies) of different values of Signal-to-Noise Ratio (SNR) is added to the deterministic computed signals and . The SNR is herein defined in decibel (dB) as the ratio between the power of the signals (i.e.

) to the power of the noise (i.e. the variance of the random variable noise ). That is SNR

, where is the average of the magnitude of the signal at different frequencies. A set of numerical simulations are conducted in which two cases are considered, the first one is the case where "all frequencies" are used for blockage identification (i.e. frequency response for the frequencies in the range with a step size ). Whereas the second one is the case where only the "resonant and anti-resonant frequencies" are used (i.e. 8 resonant and 8 anti-resonant frequencies contained in the interval

). The simulation is repeated 30 times for different values of SNR from -10dB to 40dB and Algorithm ( 1) is again used for blockage identification. The ensemble average and the standard deviation of the estimation error in the two cases are summarized in Fig. (3). Here, the estimation error is defined by [10] in which is the number of feasible estimates of the change in area . The estimates and are found by enumeration of the entire feasible region with a discretization m. 3) shows that the estimation error, is about 5 times lower when the whole frequency bandwidth is used than for the case where only the first 8 resonant and anti-resonant frequencies are used. This shows that using more frequencies increase the accuracy of the blockage detection results. This is because, for the case of zero-mean Gaussian white noise, the matched-field processing (MFP) approach has the ability of identifying the blockage parameters which maximize the SNR [START_REF] Wang | Pipeline leak detection using the matched-field processing method[END_REF]. In other words, the blockage parameters are estimated such that the corresponding signal generated from the model (i.e.

) fits the estimated mean value of the noisy measured response (i.e.

). However, the mean is only well estimated if a large enough sample size is used for its estimation. Therefore, in cases where only the resonant and anti-resonant frequencies are used for blockage identification, the mean value of the "measured" signal is not well estimated because of the lack of "measurement" data. However, for the case where more frequencies ( ) are used for the estimation of the blockage parameters, the mean value of the measured response can be better estimated and therefore a more accurate results of the blockage parameters are obtained. Note that, in this case the noise is added to the amplitude of the frequency response. Therefore, blockage detection methods which do not use the amplitude of the frequency response, (e.g. [START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF](DeSalis et al., 2000)), are not directly affected by the presence of such type of noise. However, the presence of noise can make the identification of the eigen-frequencies inaccurate. In this case, MFP method has advantage comparing to the methods based on the eigen-frequencies.

(ii) Inaccurate identification of the eigen-frequencies

There are several reasons which can make the identification of the resonant and the anti-resonant frequencies inaccurate, for instance, the measurement of the frequency response with a large frequency step. In this case, methods based on the eigen-frequencies identification become less accurate. However, since the MFP method does not rely on the identification of the resonant and anti-resonant frequencies. It can give more accurate results. To illustrate this, consider the same numerical example given above, the pressure heads and are simulated using the transfer matrix method (i.e. using Eqs ( 1) and ( 2)) in the frequency range (

) with a frequency step size . The resonant frequencies and the anti-resonant frequencies are estimated by a linear interpolation and are found to be and . Under these conditions, the blockage is identified using the MFP method and results are compared with two other methods which only use the resonant and the anti-resonant frequencies. The first method used for comparison, estimates the blockage characteristics using a minimum error in the dispersion relation (MEDR) approach [START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF]. The main advantage of this method is that it does not require the knowledge of the input signal (i.e. ) as it only uses the eigen-frequencies. However, this method relies on the accurate E-proceedings of the 37 th IAHR World Congress August 13 -18, 2017 -Kuala Lumpur, Malaysia identification of the eigen-frequencies and involves an optimization problem with at least three variables. The second method used for comparison, uses an explicit area function based on the perturbation method (AFPM) [START_REF] Desalis | The development of a rapid single spectrum method for determining the blockage characteristics of a finite length duct[END_REF]. The main advantages of this method comes from the fact that it does not assume the number of blockages and is independent from the type of excitation (i.e. ). However, the method assumes the change in area to be small. Moreover, using this method requires a large number of eigenfrequencies to guarantee an accurate area reconstruction.

The comparison between the results of the blockage identification methods in the case of inaccurate identification of the eigen-frequencies is shown in 4) shows that the MFP method has more accurate result than the other two methods. The results found using the AFPM [START_REF] Desalis | The development of a rapid single spectrum method for determining the blockage characteristics of a finite length duct[END_REF] shows large discrepancy between the reconstructed crosssectional area and the reference pipe's area. This is because only 3 resonant and 3 anti-resonant frequencies are used for the area reconstruction. Furthermore, the results of MEDR method [START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF] shows a non-accurate prediction of the blockage location and length. This is because the method is sensitive to the inaccuracy of the resonant and anti-resonant frequency identification (Duan, 2015). In contrast, the estimate of the blockage parameters was more accurate using the MFP method, this is because MFP has the ability of using frequencies other than the eigen-frequencies.

Despite the advantages of using MFP for blockage parameters estimation, this method has some limitations. The first being that the method requires a relatively accurate measurement of the frequency response in at least two measurement locations, it also requires a good knowledge of the input and of the boundary conditions. The second limitation is that the general form of the blockage size used in this method (Eq. ( 8)), to solve for the blockage location and length independently from the blockage size, does not necessarily satisfy the feasibility condition. This is because the least square solution does not necessarily satisfy the condition , and this is the reason why more than one estimate for the area is obtained by this method. In progress research is conducted to solve these limitations.

CONCLUSION

A new method for blockage detection is proposed. The method has the advantage of using information not only from the resonant frequencies but also from other frequencies, this is a key advantage especially when the signal is contaminated with a zero-mean Gaussian noise, and when the identification of eigen-frequencies is not accurate. The proposed MFP-based method solves for the location and the length of the blockage independently from the blockage's area. This gives the advantage of visualizing the overall objective function for all possible blockage parameters which cannot be achieved by other methods. However, an accurate knowledge of the input and the boundary conditions is required for the robustness of the proposed modelbased MFP method.

Figure 1 .

 1 Figure 1. Pipe system layout and characteristic lengths.

Figure 3 .

 3 Figure 3. Error bar of the estimation error with increasing SNR.

Fig. (

  Fig. (3) shows that the estimation error, is about 5 times lower when the whole frequency bandwidth is used than for the case where only the first 8 resonant and anti-resonant frequencies are used. This shows that using more frequencies increase the accuracy of the blockage detection results. This is because, for the case of zero-mean Gaussian white noise, the matched-field processing (MFP) approach has the ability of identifying the blockage parameters which maximize the SNR[START_REF] Wang | Pipeline leak detection using the matched-field processing method[END_REF]. In other words, the blockage parameters are estimated such that the corresponding signal generated from the model (i.e.) fits the estimated mean value of the noisy measured response (i.e.). However, the mean is only well estimated if a large enough sample size is used for its estimation. Therefore, in cases where only the resonant and anti-resonant frequencies are used for blockage identification, the mean value of the "measured" signal is not well estimated because of the lack of "measurement" data. However, for the case where more frequencies ( ) are used for the estimation of the blockage parameters, the mean value of the measured response can be better estimated and therefore a more accurate results of the blockage parameters are obtained. Note that, in this case the noise is added to the amplitude of the frequency response. Therefore, blockage detection methods which do not use the amplitude of the frequency response, (e.g.[START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF] (DeSalis et al., 2000)), are not directly affected by the presence of such type of noise. However, the presence of noise can make the identification of the eigen-frequencies inaccurate. In this case, MFP method has advantage comparing to the methods based on the eigen-frequencies.

  Fig. (4).

Figure 4 .

 4 Figure 4. Comparison of the Matched-field processing with other methods in the literature.

Fig. (

  Fig. (4) shows that the MFP method has more accurate result than the other two methods. The results found using the AFPM[START_REF] Desalis | The development of a rapid single spectrum method for determining the blockage characteristics of a finite length duct[END_REF] shows large discrepancy between the reconstructed crosssectional area and the reference pipe's area. This is because only 3 resonant and 3 anti-resonant frequencies are used for the area reconstruction. Furthermore, the results of MEDR method[START_REF] Duan | Extended blockage detection in pipelines by using the system frequency response analysis[END_REF] shows a non-accurate prediction of the blockage location and length. This is because the method is sensitive to the inaccuracy of the resonant and anti-resonant frequency identification(Duan, 2015). In contrast, the estimate of the blockage parameters was more accurate using the MFP method, this is because MFP has the ability of using frequencies other than the eigen-frequencies.Despite the advantages of using MFP for blockage parameters estimation, this method has some limitations. The first being that the method requires a relatively accurate measurement of the frequency response in at least two measurement locations, it also requires a good knowledge of the input and of the boundary conditions. The second limitation is that the general form of the blockage size used in this method (Eq. (8)), to solve for the blockage location and length independently from the blockage size, does not necessarily satisfy the feasibility condition. This is because the least square solution does not necessarily satisfy the condition , and this is the reason why more than one estimate for the area is obtained by this method. In progress research is conducted to solve these limitations.

ACKNOWLEDGEMENTS

This study is supported by the Hong Kong Research Grant Council (project T21-602/15R) and by the Postgraduate Studentship.


Stef Graillat
email: stef.graillat@upmc.fr

An accurate algorithm for evaluating rational functions

Keywords: oating-point, error-free transformation, rational function, Horner scheme, accuracy, rounding errors AMS Subject Classi cations: -, G, -

An accurate algorithm for evaluating rational functions

Introduction

Evaluating a polynomial or a rational function is ubiquitous in computational sciences and their applications. For example, in signal processing, transfer functions are very o en rational functions. Moreover, real functions are o en approximated by polynomials or rational functions.

In this paper, we present fast and accurate algorithms to compute the evaluation of a rational function. Our aim is to increase the accuracy at a xed precision. We show that the results have the same error estimates as if computed in twice the working precision and then rounded to working precision. is paper was motivated by papers [ , , ], where similar approaches are used to compute summation, dot product, and polynomial evaluation.

is outline of this article is as follows. In Section , we quickly recall some information on oating-point arithmetic and we give some de nitions and notations used in the sequel. In Section , we recall the compensated Horner scheme [ , ]. is algorithm makes it possible to evaluate a polynomial whose accuracy of the computed result is similar to the one given by the classical algorithm computed in twice the working precision and then rounded to the current working precision. Section is devoted to the study of the accuracy of the classic algorithm to evaluate a rational function with Horner scheme. We also de ne and compute a closed formula for the condition number of rational function evaluation. A compensated algorithm for evaluating rational functions is presented in Section . is algorithm evaluates a fractional function and gives an accuracy of the computed result that is similar to the one given by the classical algorithm computed in twice the working precision and then rounded to the current working precision. Finally, numerical experiments showing the accuracy and the performance of our new compensated algorithm to evaluate fractional functions are presented in Section .

Floating-point arithmetic

roughout the paper, we assume to work with a oating-point arithmetic adhering to IEEE oating-point standard [ ]. We assume that no over ow nor under ow occur. e set of oating-point numbers is denoted by F, the relative rounding error by u. For IEEE double precision, we have u = -and for single precision u = -.

We denote by (⋅) the result of a oating-point computation, where all operations inside parentheses are done in oating-point working precision. Floating-point operations in IEEE satisfy [ ]

(a ○ b) = (a ○ b)( + ε ) = (a ○ b) ( + ε ) for ○ = {+, -, ⋅, } and ε ν ≤ u. is implies that a ○ b -(a ○ b) ≤ u a ○ b and a ○ b -(a ○ b) ≤ u (a ○ b) for ○ = {+, -, ⋅, }. ( . )
We use standard notation for error estimations. e quantities γ n are de ned as usual [ ] by

γ n ∶= nu -nu for n ∈ N,
where we implicitly assume that nu ≤ .

Following [ ], we also use the following classic properties in error analysis (we always assume that nu < ): γ k < γ k+ and ( + u)γ k ≤ γ k+ .

One can notice that a ○ b ∈ R and a ⊚ b ∶= (a ○ b) ∈ F but in general we do not have a ○ b ∈ F. It is known that for the basic operations +, -, ×, the rounding error of a oating-point operation is still a oating-point number (see for example [ ]):

x = a ⊕ b ⇒ a + b = x + y with y ∈ F, x = a ⊖ b ⇒ a -b = x + y with y ∈ F, x = a ⊗ b ⇒ a × b = x + y with y ∈ F. ( . )
ese are error-free transformations of the pair (a, b) into the pair (x, y).

Fortunately, the quantities x and y in ( . ) can be computed exactly in oating-point arithmetic. For the algorithms, we use Matlab-like notations. For addition, we can use the following algorithm by Knuth [ , m B. p. ].

Algorithm . (Knuth [ ]). Error-free transformation of the sum of two oating-point numbers

function [x, y] = TwoSum(a, b) x = a ⊕ b z = x ⊖ a y = (a ⊖ (x ⊖ z)) ⊕ (b ⊖ z)
Another algorithm to compute an error-free transformation is the following algorithm from Dekker [ ]. e drawback of this algorithm is that we have x + y = a + b provided that a ≥ b .

Algorithm . (Dekker [ ]). Error-free transformation of the sum of two oating-point numbers.

function

[x, y] = FastTwoSum(a, b) x = a ⊕ b y = (a ⊖ x) ⊕ b
For the error-free transformation of a product, we rst need to split the input argument into two parts. Let p be given by u = -p and de ne s = ⌈p ⌉. For example, if the working precision is IEEE double precision, then p = and s = . e following algorithm by Dekker [ ] splits a oating-point number a ∈ F into two parts x and y such that a = x + y and x and y nonoverlapping with y ≤ x .

Algorithm . (Dekker [ ]). Error-free split of a oating-point number into two parts

function [x, y] = Split(a) factor = s + c = factor ⊗ a x = c ⊖ (c ⊖ a) y = a ⊖ x
With this function, an algorithm from Veltkamp (see [ ]) makes it possible to compute an error-free transformation for the product of two oating-point numbers. is algorithm returns two oating point numbers x and y such that a × b = x + y with x = a ⊗ b.

Algorithm . (Veltkamp [ ])

. Error-free transformation of the product of two oating-point numbers

function [x, y] = TwoProduct(a, b) x = a ⊗ b [a , a ] = Split(a) [b , b ] = Split(b) y = a ⊗ b ⊖ (((x ⊖ a ⊗ b ) ⊖ a ⊗ b ) ⊖ a ⊗ b )
e TwoProduct algorithm can be re-written in a very simple way if a Fused-Multiply-and-Add (FMA) operator is available on the targeted architecture [ , ]. is means that for a, b, c ∈ F, the result of FMA(a, b, c) is the nearest oating-point number of a ⋅ b + c ∈ R. e FMA satis es

FMA(a, b, c) = (a ⋅ b + c)( + ε ) = (a ⋅ b + c) ( + ε ) with ε ν ≤ u.
Algorithm . (Ogita, Rump and Oishi [ ]). Error-free transformation of the product of two oating-point numbers using an FMA.

function [x, y] = TwoProductFMA(a, b) x = a ⊗ b y = FMA(a, b, -x)

Compensated Horner scheme

We recall herea er the compensated algorithm for Horner scheme. One can nd a more detailed description of the algorithm in [ , ]. We rst recall the classic algorithm for Horner scheme and give an error bound. We then present the compensated Horner scheme together with an error bound.

e classical method for evaluating a polynomial

p(x) = n i= a i x i
is the Horner scheme which consists in the following algorithm.

Algorithm . . Polynomial evaluation with Horner's scheme

function res = Horner(p, x) s n = a n for i = n -∶ -∶ s i = s i+ ⊗ x ⊕ a i end res = s A forward error bound for the result of Algorithm . is (see [ , p. ]): p(x) -res ≤ γ n n i= a i x i = γ n p( x ) ( . )
where p(x) = ∑ n i= a i x i . It is very interesting to express and interpret this result in terms of the condition number of the polynomial evaluation de ned by cond(p, x) ∶= lim

ε→ sup p(x) -p(x) ε p(x) ∶ a i -âi ≤ ε a i , i = , . . . , n . It is well-known that cond(p, x) = ∑ n i= a i x i p(x) = p( x ) p(x) . ( . ) us we have p(x) -res p(x) ≤ γ n cond(p, x). ( . )
We can modify the Horner scheme to compute the rounding error at each elementary operation that are a sum and a product. is is done in Algorithm . .

Algorithm . (Graillat, Langlois, Louvet [ , ]). Polynomial evaluation with a compensated Horner's scheme

function res = CompHorner(p, x) s n = a n r n = for i = n -∶ -∶ [p i , π i ] = TwoProduct(s i+ , x) [s i , σ i ] = TwoSum(p i , a i ) r i = r i+ ⊗ x ⊕ (π i ⊕ σ i ) end res = s ⊕ r
If we denote by p π and p σ the two following polynomials

p π = n- i= π i x i , p σ = n- i= σ i x i ,
then one can show, thanks to error-free transformations that p(x) = s + p π (x) + p σ (x).

If one looks at the previous algorithm closely, it is then clear that s = Horner(p, x). As a consequence, we can derive a new error-free transformation for polynomial evaluation p(x) = Horner(p, x) + p π (x) + p σ (x).

e compensated Horner scheme rst computes p π (x) + p σ (x) which corresponds to the rounding errors and then adds the obtained value to the result of the classic Horner scheme Horner(p, x).

We will show that the results computed by Algorithm . admits signi cantly better error-bounds than those computed with the classical Horner scheme. We argue that Algorithm . provides results as if they were computed using twice the working precision. is is summed up in the following theorem.

eorem . (Graillat, Langlois, Louvet[ , ]). Consider a polynomial p of degree n with oating-point coecients, and a oating-point value x. e forward error in the compensated Horner algorithm is such that

CompHorner(p, x) -p(x) ≤ u p(x) + γ n p(x).
( . )

It is interesting to interpret the previous theorem in terms of the condition number of the polynomial evaluation of p at x. Combining the error bound ( . ) with the condition number ( . ) for polynomial evaluation gives

CompHorner(p, x) -p(x) p(x) ≤ u + γ n cond(p, x). ( . )
In other words, the bound for the relative error of the computed result is essentially γ n times the condition number of the polynomial evaluation, plus the unavoidable term u for rounding the result to the working precision. In particular, if cond(p, x) < γ - n , then the relative accuracy of the result is bounded by a constant of the order of u. is means that the compensated Horner algorithm computes an evaluation accurate to the last few bits as long as the condition number is smaller than γ - n ≈ ( nu) -. Besides that, ( . ) tells us that the computed result is as accurate as if computed by the classic Horner algorithm with twice the working precision, and then rounded to the working precision.

Classic evaluation of rational functions

In this Section, we present a classic algorithm to evaluate a rational function. It is based on the evaluation of the numerator and denominator (which are polynomials) with the Horner scheme. We then give a de nition of the condition number of the rational functions evaluation and give an explicit formula to compute it. We then study the numerical stability of the algorithm in terms of the condition number.

Let p, q two polynomials of degree n (it is not complicated to deal with polynomials with di erent degrees but for simplicity, we assume they both have the same degree). ey are denoted by

p(x) = n i= a i x i and q(x) = n i= b i x i .
e rational fraction f (x) is f (x) = p(x) q(x). A classic way to compute f (x) is via the evaluation of p(x) and q(x) with Horner scheme as explained in Algorithm .

Algorithm . . Rational function evaluation with Horner scheme function res = RatEval(p, q, x) res = Horner(p, x) ⊘ Horner(q, x)

e condition number of the evaluation of a rational function measures the sensitivity of the evaluation with respect to perturbations on the coe cients of the rational function. It can be de ned as follows.

De nition . . Let f (x) = p(x) q(x) be a rational function. e condition number of the evaluation of f in x is de ned by

cond( f , x) ∶= lim ε→ sup ( p q)(x) -(p q)(x) ε(p q)(x) ∶ âi -a i ≤ ε a i , bj -b j ≤ ε b j for i = , . . . , n, j = , . . . , m .
where âi and bi are respectively the coe cients of p and q.

It is possible to obtain an explicit expression for this condition number.

eorem . . Let f (x) = p(x) q(x) be a rational function. e condition number of the evaluation of f at x satis es cond( f , x) = cond(p, x) + cond(q, x).

Proof. It is easy to show that

( p q)(x) -(p q)(x) ε(p q)(x) = p(x)(q(x) -q(x)) -q(x)(p(x) -p(x)) εp(x)q(x) , = ε p(x) p(x) ⋅ q(x) -q(x) q(x) - p(x) -p(x) p(x) .
As a consequence, we have

( p q)(x) -(p q)(x) ε(p q)(x) ≤ ε p(x) p(x) ⋅ q(x) -q(x) q(x) + p(x) -p(x) p(x)
.

By de nition of p and p, we have p(x) -p(x) ≤ ε p( x ) and q(x) -q(x) ≤ εq( x ).

It follows that

( p q)(x) -(p q)(x) ε(p q)(x) ≤ p(x)q(x) p(x)q(x) ⋅ cond(q, x) + cond(p, x).
By taking the supremum and the limit for ε → , we obtain cond( f , x) ≤ cond(p, x) + cond(q, x).

To show that we, indeed, have equality, let us de ne p(x) = p(x) + ε ∑ n i= sign(a i x i )a i x i and q(x) = q(x)ε ∑ n i= sign(b i x i )b i x i . In that case, we have

( p q)(x) -(p q)(x) ε(p q)(x) = p(x)q(x) p(x)q(x) ⋅ q( x ) q(x) + p( x ) p(x) .
Taking the supremum and the limit for ε → proves that it is possible to construct a sequence of polynomials such the limit in the de nition of the condition number converges to cond(p, x) + cond(q, x).

It is now possible to evaluate the numerical stability of the classic algorithm. It is shown in eorem . that this algorithm is backward-stable.

eorem . . Let f (x) = p(x) q(x) be a rational function with oating-point coe cients, and x be a oatingpoint value. en if no under ow occurs, and res = RatEval(p, q, x),

res -f (x) f (x) ≤ u + [γ n + O(u )] cond( f , x).
Proof. We have

( p(x) q(x) ) - p(x) q(x) = ( + ε) (p(x)) (q(x)) - p(x) q(x)
with ε ≤ u. As a consequence,

( p(x) q(x) ) - p(x) q(x) ≤ u (p(x)) (q(x)) + (p(x)) (q(x)) - p(x) q(x) . ( .
)
anks to Equation ( . ), we know that (p(x)) ≤ p(x) + γ n p( x ) and q(x)γ n q( x ) ≤ (q(x)).

We can then deduce that

(p(x)) (q(x)) ≤ p(x) + γ n p( x ) q(x) -γ n q( x ) = f (x) + γ n cond(p, x) -γ n cond(q, x) , ≤ f (x) ( + γ n cond( f , x) + O(u )). ( . )
Moreover, using Equation ( . ), we have

(p(x)) (q(x)) - p(x) q(x) = q(x) (p(x)) -p(x) (q(x)) q(x) (q(x)) , = [q(x) -(q(x))] (p(x)) -[p(x) -(q(x))] (q(x)) q(x) (q(x)) , ≤ (p(x)) (q(x)) ⋅ q(x) -(q(x)) q(x) + p(x) q(x) ⋅ p(x) -(p(x)) p(x) , ≤ (p(x)) (q(x)) γ n cond(q, x) + p(x) q(x) γ n cond(p, x) .
Now, using Equations ( .) and ( . ), we can deduce that

( f (x)) -f (x) ≤ u f (x) ( + γ n cond( f , x) + O(u ))+ f (x) ( + γ n cond( f , x) + O(u ))γ n cond(q, x) + f (x) γ n cond(p, x),
that can be simpli ed into

( f (x)) -f (x) ≤ f (x) ⋅ [u + (γ n + O(u )) cond( f , x)],
which concludes the proof.

A compensated algorithm for evaluating rational functions

In this Section, we present a compensated version of the classic algorithm for evaluating rational functions. e idea is to replace the Horner scheme used to evaluate the numerator and denominator by the compensated Horner scheme followed by a oating-point division. is is Algorithm . . Algorithm . . Rational function evaluation with compensated Horner scheme function res = CompRatEval(p, q, x) res = CompHorner(p, x) ⊘ CompHorner(q, x)

e following eorem shows that the error bound on the accuracy of the computed result given by the compensated algorithm is improved compared to the one for the classic algorithm. eorem . . Let f (x) = p(x) q(x) be a rational function with oating-point coe cients, and x be a oatingpoint value. en if no under ow occurs, and res = CompFracEval(p, q, x),

res -f (x) f (x) ≤ u + O(u ) + [ γ n+ + O(u )] cond( f , x). ( . )
Proof. For better readability, we will denote CH(p, x) the computed result of CompHorner(p, x). By de nition we have

CompRatEval(p, q, x) -f (x) = ( CH(p, x) CH(q, x) ) -f (x) = ( + ε) CH(p, x) CH(q, x) -f (x)
with ε ≤ u. As a consequence,

CompRatEval(p, q, x) -f (x) ≤ u CH(p, x) CH(q, x) + CH(p, x) CH(q, x) -f (x) . ( . )
Moreover, using Equation ( . ), we have CH(p, x) ≤ p(x) ( + u) + γ n p( x ) and (u) q(x)γ n q( x ) ≤ CH(q, x) .

As a consequence, we can bound

CH(p, x) CH(q, x) ≤ p(x) q(x) ⋅ + u + γ n cond(p, x) -u -γ n cond(q, x) , ≤ f (x) ⋅ [ + u + γ n cond(p, x)] ⋅ [ + u + γ n cond(p, x) + O(u )], ≤ f (x) ⋅ [( + u) + γ n cond( f , x) + O(u )]. ( . ) Moreover, CH(p, x) CH(q, x) -f (x) = CH(p, x)[q(x) -CH(q, x)] -CH(q, x)[p(x) -CH(p, x)] q(x) CH(q, x) , ≤ CH(p, x) CH(q, x) ⋅ q(x) -CH(q, x) q(x) + p(x) q(x) ⋅ p(x) -CH(p, x) p(x) , ≤ CH(p, x) CH(q, x) (u + γ n cond(q, x)) + f (x) (u + γ n cond(p, x)).
Using Equation ( . ), we deduce that

CH(p, x) CH(q, x) -f (x) ≤ f (x) ⋅ [( + u) + γ n cond( f , x) + O(u )] ⋅ (u + γ n cond(q, x))+ f (x) (u + γ n cond(p, x)).
which can be simpli ed into

CH(p, x) CH(q, x) -f (x) ≤ f (x) ⋅ [ u( + u) + γ n+ cond( f , x) + O(u )]. ( . )
Finally combining Equations ( . ) and ( . ) with Equation ( . ), we obtain

CompRatEval(p, q, x) -f (x) ≤ u f (x) ⋅ [( + u) + γ n cond( f , x) + O(u )]+ f (x) ⋅ [ u( + u) + γ n+ cond( f , x) + O(u )].
that can be rewritten as

CompRatEval(p, q, x) -f (x) ≤ f (x) [ u( + u) + ( γ n+ + O(u )) cond( f , x)].
and so concludes the proof.

In other words, the bound for the relative error of the computed result is essentially γ n+ times the condition number of the rational function evaluation, plus the term u for rounding the result to the working precision.

In particular, if cond( f , x) < ( γ n+ ) -, then the relative accuracy of the result is bounded by a constant of the order of u. is means that the compensated algorithm computes an evaluation accurate to the last few bits as long as the condition number is smaller than ( γ n+ ) -≈ (( n + )u) -. Besides that, ( . ) tells us that the computed result is as accurate as if computed by the classic rational function evaluation algorithm with twice the working precision, and then rounded to the working precision.

Numerical experiments

e numerical experiments have been done on a laptop with an Intel Core i processor at . GHz with Gb of RAM. We used MATLAB R b.

We test the rational function f n (x) = p n (x) q n (x) where p n is a random polynomial of degree n and q n is the expand form of the polynomial (x -) n . e argument x is chosen near to the unique real root of q n , and with many signi cant bits so that a lot of rounding errors occur during the evaluation of q n (x). We increment the degree n from until a su ciently large range has been covered by the condition number cond( f n , x).

Here we have

cond( f n , x) = cond(p n , x) + qn ( x ) q n (x) = cond(p n , x) + + x -x n ,
for x close to and cond( f n , x) grows exponentially with respect to n. In the experiments reported on Figure , cond( f n , x) varies from to (for x = ( . ), that corresponds to the degree range n = , . . . , ). ese huge condition numbers have some meaning since here the coe cients of p n and q n and the value x are chosen to be exact oating-point numbers.

We experiment both RatEval and CompRatEval. For each rational function f n , the exact value f n (x) is approximate with a high accuracy thanks to the Symbolic Math Toolbox of MATLAB. We observe that the compensated algorithm exhibits the expected behavior. e full precision solution is computed as long as the condition number is smaller than u -≈ . en, for condition numbers between u -and u -≈ , the relative error degrades to no accuracy at all.

We now demonstrate the practical e ciency in terms of running time comparing our algorithm and up-to-date challenger.

Since Bailey's double-double are usually considered as the most e cient portable library to double the IEEEdouble precision, we consider it as a reference in the comparisons. Double-double numbers are represented as an unevaluated sum of a leading double and a trailing double. More precisely, a double-double number a is the pair (a h , a l ) of oating-point numbers with a = a h + a l and a l ≤ u a h .

In the sequel, we present two algorithms to compute product of two double-double or a double times a double-double. ose algorithms are taken from [ ].

Algorithm . . Product of the double-double number (a h , a l ) by the double number b

function [c h , c l ] = prod dd d(a h , a l , b) [s h , s l ] = TwoProduct(a h , b) [t h , t l ] = FastTwoSum(s h , (a l ⊗ b)) [c h , c l ] = FastTwoSum(t h , (t l ⊕ s l ))
Algorithm . . Addition of the double number b to the double-double number (a h , a l )

function [c h , c l ] = add dd d(a h , a l , b) [t h , t l ] = TwoSum(a h , b) [c h , c l ] = FastTwoSum(t h , (t l ⊕ a l ))
e double-double library can be used to implement an Horner scheme in quadruple precision like DDHorner.

Algorithm . . Horner scheme with internal double-double computations

function res = DDHorner(p, x) s h = a n s l = for i = n -∶ -∶ [p h , p l ] = prod dd d(s h , s l , x) [s h , s l ] = add dd d(p h , p l , a i ) end res = s h
With DDHorner, we can evaluate a rational function in quadruple precision.

Algorithm . . Rational function evaluation with double-double Horner scheme function res = DDRatEval(p, q, x) res = DDHorner(p, x) ⊘ DDHorner(q, x)

We have implemented the three algorithms RatEval, CompRatEval, and DDRatEval in a C code to measure their overhead compared to the RatEval algorithm. We have programmed these tests straightforwardly with no other optimization than the ones performed by the compiler.

For each algorithm, we measured the ratio of its computing time over the computing time of the classic rational function evaluation algorithm. It turned out that our compensated algorithm CompRatEval is about times slower than the classic Horner scheme. e same slowdown factor is about for algorithm DDRatEval.

From a practical point of view, we can state that our algorithm is about times faster than the algorithm with double-doubles.

We compared RatEval, CompRatEval and DDRatEval in term of measured computing time. We tested with random rational functions where the degree of numerator and denominator vary from to .

Conclusion

We presented a fast algorithm for the evaluation of rational functions in oating-point arithmetic. We have proved that the accuracy of the result computed by our compensated algorithm is similar to the one given by the classic algorithm performed in doubled working precision. e only assumption we made is that the oating-point arithmetic available on the computer is conformed to the IEEE-Standard. Its low requirement make it highly portable, and our compensated algorithm could be easily integrated into numerical libraries. Our algorithm uses only basic oating-point operations, and only the same working precision as the data. Finally, our compensated algorithm runs much more faster than existing implementation producing the same output accuracy. is approach can easily be generalized to rational functions where numerators and denominators are written in other basis and for bivariate rational functions (see [ , , , ] for example).

  Figure presents the relative accuracy yf n (x) f n (x) of the evaluation y computed by the two algorithms.

Figure :

 : Figure : Comparison of the accuracy of RatEval and CompRatEval

Table :

 : Measured computing times with RatEval normalised to .

	n RatEval CompRatEval DDRatEval
	.	.	.
	.	.	.
	.	.	.
	.	.	.
	.	.	.

Acknowledgement

is work was partly supported by the project FastRelax ANR--CE --.