J Ragot

D Maquin

G Bloch

DATA RECONCILIATION FOR PROCESS FLOW

Keywords: data reconciliation, redundancy, observability, error detection, variance estimation, model based fault detection

Before improving the control of a plant, we must make sure of information coherency issued from instrument lines or sensors. In fact, this information can be corrupted by errors and can also deviate from the optimum functioning range. Various methods for the detection and location of gross errors in process data have been proposed in recent years including the parity space approach, the standardized least square residuals approach and the standardized imbalance residuals approach.

In this survey, the authors try to summarize the various aspects of data reconciliation, to point out the main difficulties and to present the state of the art in this field. The authors present the steps of the data reconciliation problem in the following order : techniques of data reconciliation, classification of the data by the observability concept, gross error detection and localisation, variance of measurement error estimation and sensor positioning. The whole presentation is applied to linear systems nevertheless some extensions are given for non linear systems.

INTRODUCTION

If the measurement errors are normally distributed with zero mean and known diagonal variance matrix V, the probability density function of the measurements Z can be written as : The validation of signals is a technique which integrates information from redundant and from functionally diverse sensors to provide highly reliable information to operating crews and to automatic controllers. Signal validation is generally performed by like-sensor comparisons (direct redundancy). When the increase of sensors is impossible, we prefer to use analytical redundancy ; analytical redundancy refers to the physical relationships, such as conservation of mass or conservation of energy, that exist among the many variables being measured in a system. P(Z) = 1 (2π) m/2 |V| 1/2 exp( -

1 2 || Z -H(X * ) || 2 V -1 ) (1.3)
The maximum likelihood estimator (X ^, θ ^) of the true values (X * , θ * ) maximises the probability density function P, subject to the model constraints.

The difficulty in solving this problem is closely related to the structure and the dimension of model equations as well as to the number of observations. In the case of linear constraints and measurement equations, analytical solutions can be found. In the other cases, one must use techniques of hierarchical calculus, of linearization, of variable substitution or even approximate solutions.

-DATA RECONCILIATION

The data reconciliation problem can be formulated in very simple terms. A set of measurements does not check process functioning equations. How may one correct (or reconcile) the measurements in order to force them to verify this set of equations which are supposed to be exact ? 1.2 -Linear models case In practise, the formulation of the problem is not so easy. Indeed one is faced with the choice of hypothesis about the statistical distribution of the measurement errors, the system dimension (number of variables and constraints), the nature of the constraints (static or not, linear or not), the process state (steady or not), and incomplete or imperfect knowledge of the model structure or of the parameters.

The linear models case can seem to be very restrictive. Nevertheless, industrial applications frequently call for this class of systems, which represents systems described physically by balance equations in total material or energy flows. Moreover, an approximate solution to problems described by multilinear equations can be found, using linear formulation, by decoupling multilinear equations in linear equations. Generally, a process is made up of interconnected subsets and can be represented by an oriented graph, where the v streams represent the flows and the n nodes processing units or junction points or several streams. Each node is associated to a balance equation, drawn from the conservation laws of material or energy. Thus, the complete set of the balance equations can be written as :

Reconciliation techniques have been applied in various fields, including chemical processes (see [START_REF] Clair | Measurement data set improvment through material balance calculation. Application to an industrial chemical plant[END_REF], [START_REF] Holly | Reconciliation of mass flow rate measurement in a chemical extraction plant[END_REF]) and mineralurgy [START_REF] Cutting | Estimation of interlocking mass balances on complex mineral beneficiation plants[END_REF]).

-Problem formulation

Generally, the systems in question can be described in two parts. The first is the set of static constraints and the second is the measurement equation :

M X * = b (1.4)
where M is the (n.v) incidence matrix, associated to the system.

F(X * , θ * ) = 0 (1.1)
If the variables are measured with linear devices, the measurement equation is simply :

Z = H(X * ) + ε (1.2) Z = H X * + ε (1.5) I II III 1 3 7 6 4 IV 9 8 x x x x x 2 5
A global resolution of problem (1.6) by the Lagrange multipliers technique leads to the unbiased estimator :

X ^ = P X = ( I -G -1 M T (M G -1 M T ) -1 M ) X (1.7) Figure 1 : a simple network G -1 = H T V -1 H + M T M (1.8)
In this commonplace example, we easily find that variables 1, 2, 4 and 5 are bound by a redundancy equation (an equation with only measured variables). These variables are said to be overdetermined and will be corrected by the previous techniques. On the contrary, variable 7 could not be estimated. Variables 3 and 6 can be easily deduced from equation of nodes I and II for the first one and III for the second one. The lack of information prohibits any estimation of variables 8 and 9. In fact, there is only one equation to determine two unknown quantities.

In practise, for low dimension and well-conditioned systems, the computation of X ^ by equation (1.7) does not pose any numerical difficulty. For large scale systems (for instance, more than hundred nodes), others techniques must be employed, including constraints elimination, recursive estimation, decentralized calculus, sparse matrices techniques, orthogonal transformations... For instance, the search for the minimum of a function subject to constraints can be replaced by the search for the minimum of a function not subject to constraints, by reducing the number of variables. Another technique consists of progressively incorporating the constraints.

The measurements are projected in the first constraint subspace ; the estimation obtained is then projected in the first two constraints subspace ...

In general the variables can be classified into two essential categories (i) the observable variables : this group contains the measured and the deducible variables, (ii) the unobservable variables : the lack of information prevents their estimation.

Observability

We consider a linear system under steady-state condition defined by equations (1.4) and (1.5). This system is said to be globally observable if the knowledge of the measurements Z and of the model constraints allow one to find a single estimation X ^ when ε = 0. This system of equations has a single solution only if :

-Bilinear models case

In this case, the system is defined by an incidence matrix M and two kinds of values X * and Y * i (i = 1, ..., nc) to be estimated (for instance, in mineralurgical processes, material flow rates and concentrations or contents, respectively, linked by total and partial material balance equations) : rank

      H M = dim(X * ) = v (2.1) M X * = 0 (1.9a)
If all the variables are not observable, it is necessary to determine the observable subsystem to which data reconciliation can be applied.

M X * * Y * i = 0 i = 1, ..., nc (1.9b)

where the symbol * denotes the Hadamard product : the jth component of the vector X*Y is the product of the jth components of X and Y.

Redundancy

The linear system described by (1.4) and (1.5) is said to be redundant when the data exceeds the minimum necessary for a single determination of the system state X * . This definition can be translated by an inequality on the dimensions of the different matrices. The system is redundant if n + m > v and its degree of redundancy r is equal to n + m -v.

The validation problem consists of seeking the minimum with respect to X ^ and Y ^i (i = 1, ..., nc) of the probability density function in X and Y.

Several methods can be used to solve this optimisation problem : an approximate suboptimal method which neglects the coupling between balance equations, a method based on the linearization of balance equations, a method using estimation by relaxation. In the case of multistage operations, the constraint equations take the form of a tridiagonal matrix and this particular sparse form is used to reduce the number of equations to be solved simultaneously in an iterative procedure. We note that the bilinear formulation can be extented to multilinear cases. But even in the bilinear case, the linearized systems become rapidly of large scale. So, for multilinear systems, one is well advised to solve the optimisation problem by a hierarchical decomposition of the calculations, for instance by variable type. The resulting direct iterative scheme, which is rather difficult to present with conciseness, can be found in Ragot (1990).

Practical method of decomposition

Several methods can be used to perform this decomposition (Darouach, 1986) singular values decomposition, reduction of the incidence matrix to an echelon form etc... This last method is the simplest and the most efficient because it uses the particular structure of the incidence matrices which contain only ± 1 and 0 elements (Maquin, 1987). After some elementary row and column operations, the incidence matrix can be written as :

A I 0 0 A 0 I B M 0 0 0 r n -r m r v-m-r r r 1 2 r 1 2

Deduction equations

Redundancy equations

Observable variables

Non observable variables

-OBSERVABILITY

-Introduction

The different techniques which have been previously presented can be applied to systems described by algebraic linear or non linear equations when all the variables are measured. Technological and economical constraints generally prevent the measurement of all the variables. The fragmentary character of information prohibits the immediate usage of the previous methods. It is then necessary to do a preliminary quantitative analysis of all the available information in order to determine, before any calculus, those which can The examination of this matrix shows : redundancy equations to null columns of M r and cannot be corrected), unmeasured and deducible variables (these variables are deduced from those of the previous classes), unmeasured and non estimable variables : these variables do not appear in the deduction equations and cannot be corrected ; further measurements are necessary to make these variables observable.

(Q M 1 ⊗Y m1 ) X m _ 1 + (Q M 3 ⊗X m1 ) Y m _ 1 = Q d 2 (2.8b)
With R, a regular matrix defined by : R Q M 3 = 0, the system (2.8) may be transformed to :

Q M 1 X m _ 1 = Q d 1 (2.9a) (R Q M 1 ⊗Y m1 ) X m _ 1 = R Q d 2 (2.9b) 2.

-The bilinear case

In the case of only one component, the bilinear system can be described by equations (1.9) with i=1. Equation (1.9b) can also be written as :

Four steps are then necessary to study the observability : a) System (2.9) enables one to extract the observable part of X m

_ 1 . M X * Y = (M ⊗ X) Y = (M ⊗ Y) X b)
Using the already known observable part of X m _ 1 , equation (2.8) allows the determination of the observable part of Y m _ 1 . where the ⊗ operator is applied between a matrix and a vector in order to form a matrix in which the line k is obtained by multiplying each of the terms of the line k of the matrix by each term of the vector.

c) The observable part of X m _ 2 is obtained from the knowledge of observable parts of X m _ 1 and Y m _ 1 , and equation (2.4a). The two measurement vectors X and Y can be partitioned into measured and unmeasured parts : 1986) have developed a quite similar approach to solve this problem. They proposed the construction of two projection matrices in order to decompose the problem into subproblems to be solved in sequence. The first matrix eliminates all unmeasured component flow rates and concentrations from the equations ; the second then removes the unmeasured total flow rates. To solve this problem, [START_REF] Romagnoli | On the rectification of measurement errors for complex chemical plants, steady-state analysis[END_REF] introduced a structural representation of the system. Using the properties of the structural matrices, which only take into account the occurrences of a variable in the equations, they proposed a classification of the variables. Their analysis uses the generic rank of the structural matrices and an algorithm to assign one equation to any unmeasured variables. It is important to note that these previous works require the composition of each stream to be either completely measured or not measured at all.

d) Similarly, equation (2.4b) is used to extract the observable part of Y m _ 2 . X =           X m X m _ Y =           Y m Y m _ (2.2) Crowe et al. (
These decompositions allow the classification of the variables into four distinct groups. According to the above partition (2.2), the incidence matrix M is partitioned as well :

M = [ M 1 M 2 M 3 M 4 ] (2.3)
The system (1.9) can then be written as : They make use of graph theoretical concepts and proceed in a layered approach through the various graphs derived from the process graph.

M 1 X m _ 1 + M 2 X m _ 2 = -(M 3 X m1 + M 4 X m2 ) = d 1 (2.4a) M 1 X m _ 1 *Y m1 + M 2 X m _ 2 *Y m2 + M 3 X m1 *Y m _ 1 = -M 4 X m2 *Y m2 = d 2 (2.
or in matrix notation :

O (X m ,Y m )             X m _ 2 X m _ 1 X m _ 2 * Y m _ 2 Y m _ 1 =         d 1 d 2 (2.5)

-GROSS ERROR DETECTION with

This part of the paper presents methods of detection and location of measurement failures which are also called gross or large errors. Some of the techniques presented use linear models but some others can be applied to non-linear models. We assume that all variables are measured (redundant system). The general procedure of error detection is divided into two parts :

O (X m ,Y m ) =         M 2 M 1 0 0 0 M 1 ⊗Y m1 M 2 M 3 ⊗X m1 (2.6)
The system will be observable if :

rank O (X m ,Y m ) = dim             X m _ Y m _ (2.7)
the generation of so-called residuals, which are functions of measurements that are accentuated by the errors, the detection, the isolation and the estimation of the error.

In the following, we present the balance residuals analysis (known as method of pseudonodes), the parity space approach, the analysis of residuals (known as method of measurement test), the analysis of the sum of squares deviations. Most of the methods which have been developed for gross error detection involve the use of statistical tests based on the assumption that the measurement errors are gaussian. After generating the set of residuals, tests for outliers are applied ; any measurement for which the residual fails the test is considered to contain a gross error.

If the system is not globally observable, a decomposition is necessary. We will show that the study of observability of bilinear systems can be reduced to a sequential study of observability of linear systems (Maquin, 1987).

For the study of the observability of X m _ 1 , let us introduce two matrices Q and R which allow the elimination of the unmeasured

terms X m _ 2 and Y m _ 2 in equation (2.4). Let us consider Q, an
With the previous hypothesis of a gaussian distribution of the measurement errors, one shows that the vector R follows a normal distribution with zero mean and covariance V R : V R = M V M T . In order to compare the components of the R vector, let us define a normalized imbalanced vector R N whose component R N (i) is defined by :

P j = W j T P || W j || (3.4)
When the kth sensor is faulty, then the parity vector P follows the direction of the kth column of W which is the greatest projection of P on the axis. Then after the detection of one or several failed sensors, we locate those with the greatest projection. Next, we delete the suspect sensor and calculate the detection test after the deletion of each sensor. We stop the procedure when the magnitude of the parity vector P corresponding to the remaining sensors no longer satisfies the detection test.

R N (i) = R(i) V R (i,i) for i = 1, ..., n (3.1)
Each component R N (i) has a normal distribution with a zero mean and unity variance. Then, a simple statistical test criterion of data inconsistency can be used. From a cumulative normal distribution table the probability of R N (i) being, for example, in the interval of 1.96 to 1.96 is read to be 0.95. Therefore, when |R N | > 1.96, we might say that the inconsistency is significant with a probability of 0.95.

-Standardized least square residuals analysis

The so-called measurement test is based on statistical tests for outliers applied to the residuals between measurement and estimations. The different steps are :

If R N (i) exceeds the critical value, this denotes that node i is a bad node. If we assume the presence of only one gross error which affects only one stream of the process, it has been established that the node which contains this stream has the bigger residual. Then, a simple examination of the components of R shows the suspect stream.

Step 1. Apply the reconciliation procedure to the data of the process and compute the estimate vector X ^, the least square residuals vector E and their variance matrix V E .

Step 2. Compute the standardized least square residuals by : E N = diag(V E ) -1 E (3.5) When several gross errors are present, their location is more complicated. It has been proposed (Mah, 1982) to apply the preceding test to each node and also to the aggregates of two or more nodes, which are called pseudonodes. Note that the main assumption underlying the method is that the gross errors in two or more measurements do not cancel each other. The location of the suspect stream necessitates the fusion of not necessarily linear equations. Notice that these fusions, which correspond to the elimination of one variable between two equations, are not easy from an analytical point of view and, what is more, these fusions can not always be achieved, nor are they unique. See Kratz (1990) for an example.

Under the hypothesis that the measured values contain no gross error, E N has a unit variance matrix.

Step 3. Detect and locate the failed measurement : as for the residuals of the imbalance, each E N (i) is compared with a critical test value (defined by the overall probability of type I error). If at least one component of E N is out of the confidence interval, this denotes the presence of a bad stream. Locating the failed sensor one proves that it corresponds to the greatest standardized least square residual.

For the linear case, instead of (3.5), [START_REF] Tamhane | A note on the use of residuals for detecting an outlier in linear regression[END_REF] has shown that for a non diagonal covariance matrix V, a vector of test statistics with the maximal power for detecting a single gross error is obtained by premultiplying R by V -1 . The power of the test (the probability of correctly detecting and identifying gross errors when they are present in the process data) has been established and discussed by [START_REF] Iordache | Performance studies of the measurement test for detection of gross errors in process data[END_REF] under different conditions (various networks, position of the error, variance values ...) using a Monte Carlo simulation. Note that Jongenelen (1988) pointed out the case where the variance V depends on an unknown scale factor σ 2 .

-Parity space approach

In the absence of gross errors, the measurements depend on the true values following the linear relation (1.5). The so-called parity space is orthogonal to the space taken from the columns of H. The parity vector is related to the measurement vector Z through a projection matrix W [START_REF] Potter | Thresholds redundancy management with arrays of skewed instruments[END_REF] of dimension n.v (n=v-m) :

P = W Z (3.2)
From (1.5), written with direct measurement (H=I), the expression (3.2) yields : Some of the preceding results can be applied for non linear systems. P = W ε (3.3) 3.4 -Residual criterion analysis Consequently, the parity vector is independent of the measurement vector and contains only the errors due to the faults. As an extension it is possible to consider the constrained case, where the true values are constrained by the linear model.

In 1988, [START_REF] Crowe | Recursive identification of gross error in linear data reconciliation[END_REF] identified the gross errors (in the linear case) by examinating the objective function. A recursive algorithm has been proposed which results from the deletion of suspect measurements. It is shown that if the deletion of a single measurement decreases the objective function (in respect to a statistical test), this measurement corresponds to a gross error. Equation (3.3) shows that for normal functioning, the magnitude of parity vector is small (presence of measurement noise). If a failure occurs in only one of the sensors, then the parity vector may grow in a fixed direction associated with the failed sensor.

As was pointed out by [START_REF] Reilly | Application of statistical theory of adjustment to material balances[END_REF], the quantity Φ r = E T V -1 E has a chi-square distribution with the number of degrees of freedom equal to the rank ok M. Thus the imbalances can be globally tested against tabulated values of chi2.

Moreover the components of the parity vector have the same probability distribution as the measurement errors which are independent gaussian with a zero mean value. As the variable c 2 = P T P is the sum of the square of (v-m) normally distributed variables, it has a chi-square probability distribution with (v-m) degrees of freedom and may be compared to the threshold c 2 1-α where c 2 1-α is the value A difficulty with this global test is that while it indicates well the presence of gross errors it is not able to identify the source of these errors. However, the use of a sequential procedure allows one to [START_REF] Rosenberg | Evaluation of schemes for detecting and identifying gross errors in process data[END_REF] which also proposes two composite tests (extended measurement test, dynamic measurement test).

fying which of the measurable parameters have to be measured. Another objective is to select an optimal measurement structure to enhance estimation accuracy. A comparison of the different methods and with several others can be found in [START_REF] Serth | Gross error detection and data reconciliation in steam-metering system[END_REF]. Another approach has been formulated by Yamamura (1988) also using a serial elimination ; a new criterion test has been tested which is based on the Akaike's information function and in order to minimize the amount of computation a branch-and-bound algorithm is used.

More recently, Kretsovalis (1987) developed a new approach for steady-state systems with measurement noise only. This work point out the problem on the accuracy of the estimates obtained by data reconciliation. Quantitative relations have been developed for the effects of adding and removing a single measurement on the estimation accuracy. These relationships have been utilized in developing evolutionary strategies for selecting an optimal measurement structure. For the sake of brevity, all these different approaches cannot be covered. We restricted ourselves to the presentation of an original method which can be used for the placement of sensors in systems described by linear and, by extension, bilinear equations.

-ESTIMATION OF MEASUREMENT VARIANCE

Most of the reconciliation techniques are based on the assumption that the measurement errors are random variables obeying a known statistical distribution. Almost without exception the techniques start with a given known covariance matrix of measurement errors. [START_REF] Almasy | Estimation of measurement error variances from process data[END_REF] has proposed a method of estimating this matrix which makes use of available data and takes into account the reconciliation point of view. Recently, Ragot (1990) proposed also a method which allows simultaneously the estimation of the variance matrix and the reconciliation of the data on several time intervals.

-The linear case

Two distinct questions can be formulated : (i) knowing that different system measurement points exist and that some streams are unmeasurable, is it possible to attain total observability by adding further sensors ? and (ii) fixing a cost function that attributes to each streams a given weight, which group of streams then must be measured to obtain total observability by minimizing the measurement cost ? This given procedure is referred as a direct method due to the fact that the covariance matrix and the reconciliation of data are estimated at the same time. For a series of p steady state point X * j :

The first problem is easily solved by fixing boolean constraints which specify absolutely whether or not measurements can be made. The second one requires the introduction of continuous constraints.

X ij = X * j + ε ij j = 1, ..., p i = 1, ..., m j (4.1)

M X * j = 0 j = 1, ..., p (4.2) m = Σ p j=1 m j
The vector of true values X * is decomposed into three distinct categories : X * m : vector of measured streams, X * f : vector of unmeasured but measurable streams (free stream), X * m _ : vector of unmeasured and unmeasurable streams. The maximisation of the likelihood function of the measurements gives the optimality equations :

We try to solve the following problem : is the knowledge of the measurements X * m sufficient to obtain total observability of the process ? Elsewhere, can this observability be obtained in placing supplementary sensors on streams defined by X * f ?

V = 1 m diag ( Σ p j=1 Σ m j i=1 (X ^j -X ij ) (X ^j -X ij ) T ) (4.3) X ^j = (I v -VM T (MVM T ) -1 M) 1 m j Σ m j i=1 X ij j = 1, ..., p (4.4)
According to the X * decomposition, the incidence matrix of the network may also be decomposed as :

A hierarchical solution for estimating V and X ^j has been proposed (Ragot 1990) : it uses a structure with two levels of hierarchy :

M = ( M m | M f | M m _ ) (5.1)
The basic idea (Maquin, 1986) is to apply the transformation described in the second part of this paper independently to M f and to M m _ . After extracting regular parts from M m _ and applying a suitable transformation to the result, in order to obtain an identity matrix, the incidence matrix can be written as :

-The first level calculates X ^j by taking V to be known (evalued in the second level of calculation).

-The second level calculates V as a function of the estimation X ^j transmitted by the first level.

The two levels of calculation are repeated until the estimations converge.

M' =           A 1m A 1f 0 A 1m _ I A 2m A 2f I 0 0 A 3m 0 0 0 0
(5.2) 5 -SENSOR POSITIONING Data validation, as previously described, is considered as an analysis stage and it gives the user coherent statistical information. Additionally, the above may be completed by a more ambitious study introducing modifications of the instrumental scheme. Firstly, it is necessary to analyse the lack of information and to make the user aware of the instrumental inadequacy of certain parts of the process. Secondly, in order to make the process information fuller, an understanding of its deficiency enables the satisfactory localization of sensors. However, in practice, the introduction of such supplementary sources of information must take into account economic and technical constraints. This matrix can be read in the following manner. Remembering that the subscripts m, f and m _ indicate respectively the measured streams, the free one and the unmeasurable one, line 3 of M' corresponds to the redundancy equations. Line 2 does not contain unmeasured variables. Therefore it is sufficient to fix the "f" variables intervening in the second column, to make the variable in the third column deducible. Line two is then only composed by deduction equation. Then, the first line only contains unmeasurable varia-_

-The bilinear case

[KRET 87] A. KRETSOVALIS, R.S.H. MAH, "Observability and redundancy classification in multicomponent process networks", AIChE J., vol. 33, n° 1, p. 70-82, 1987. The optimal location of measurements in linear systems is based on the counting of equations which contain more than two unknown quantities ; according to the nature of these unknown quantities (measurable or not measurable), supplementary sensors can be placed to render the system observable.

[KRET 88] A. KRETSOVALIS, R.S.H. MAH, "Observability and redundancy classification in generalized process networks, 1 : theorems, 2 : algorithms", Comp. Chem. Eng., vol. 12, n° 7, p. 671-703, 1988. The bilinear case can be treated according to this principle. However, this is more complicated in practise because of the interaction of X and Y. This means that the observability of X or Y variables are simultaneously modified by the addition of a sensor for X or Y variables. This interaction between the different variables can be summarized by the following three rules :

[MACD 88] R.J. MACDONALD, C.S. HOWAT, "Data Reconciliation and Parameter Estimation in Plant Performance Analysis", AIChE J., vol. 34, n°1 p. 1-9, 1988.

[MAH 82] R.S.H. MAH, A.C TAMHANE, "Detection of gross errors in process data", AIChE J., vol. 28, n° 5, p. 828-830, 1982. 1-An equation which contains only one unknown X quantity is observable. The X unknown quantity is directly deduced from a linear equation.

[MAQU 86] D. MAQUIN, J. FAYOLLE, M. DAROUACH, J. RAGOT, "Le positionnement de capteurs dans les systèmes industriels de grande dimension", Proc. of IMACS cong., Villeneuve d'Ascq, 1986. 2-An equation which contains only one unknown Y quantity (with all the X variables known) is observable. The unknown Y quantity can be deduced from a bilinear equation.

[MAQU 87] D. MAQUIN, "Observabilité, diagnostic et validation de données des procédés industriels", Thèse de l'Université de Nancy I, novembre 1987. 3-An equation with all the X variables known except two of them and all the Y variables known is observable. The two unknown quantities can be deduced from the set of linear and bilinear equation.

[NARA 87] S. NARASIMHAN, R.S.H. MAH, "Generalized likelihood ratio method for gross error identification", AIChE J., vol. 33, n° 9, p. 1514-1521, 1987. The principle of the location of measurements is directly derived from this analysis. The addition of a new sensor is carried out in order to obtain one of the three previous situations. As the addition of a sensor modifies the observability of the system, this procedure had to be repeated until total observability is obtained.

[NOGI 72] S. NOGITA, "Statistical test and adjustment of process data", Ind. Eng. Chem. Proc. Des. Dev., vol. 11, n° 2, p. 197-200, 1972. [PATT 89] R.J. PATTON, P.M. FRANK, R.N. CLARK, "Fault diagnosis in dynamic systems : theory and application". Prentice Hall, 1989.

  4b) More recently, Kretsovalis et al. (1988) presented a classification algorithm which does not impose any restriction on the number and the localization of measurements. They included in their treatment chemical reactions, flow splitters and pure energy flows.


Qian Feng

Sing Kiong Nguang
email: nguang@auckland.ac.nz.

Alexandre Seuret
email: aseuret@laas.fr.

Stability analysis of linear coupled differential-difference systems with general distributed delays

Keywords: Distributed Delay, Integral Inequality, Coupled differential-difference systems, Approximations. I. INTRODUCTION

We present a new approach for the stability analysis of linear coupled differential-difference systems (CDDS) with a general distributed delay. The distributed delay term in this note can contain any L 2 function which is approximated via a class of elementary functions including polynomial, trigonometric and exponential functions etc. Through the application of a new proposed integral inequality, sufficient condition for the stability of the system is derived in terms of linear matrix inequalities based on the construction of a Liapunov Krasovskii functional. The methods proposed in this note can handle problems which cannot be deal with by existing approaches. Two numerical examples are presented to show the effectiveness of our proposed stability condition.

specific class of elementary functions. This motives us to develop a more efficient and robust approach to handle the presence of general distributed delay terms in a system.

In this paper, we propose a new approach to handle the problem of the stability analysis of a linear CDDS with a general distributed delay term. Specifically, the distributed delay kernels can be any L 2 function over an interval. Moreover, the models of many existing delay systems, such as the ones in [START_REF] Gu | Lyapunov Krasovskii functional for uniform stability of coupled differential-functional equations[END_REF], [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF], [START_REF] Feng | Stabilization of uncertain linear distributed delay systems with dissipativity constraints[END_REF], [START_REF] Münz | Robust stabilization and H ∞ control of uncertain distributed delay systems[END_REF], are generalized by the model we considered. The distributed delay kernels are handled by a new approximation scheme in this paper based on a class of elementary functions which includes the Legendre polynomials approximation scheme proposed in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF]. Moreover, our approximation scenario also provides a unified matrix framework which can easily handle the situation of a CDDS having multiple distributed delay kernels. It shows that, if orthogonal functions {ϑ i (•)} d i=1 are chosen as the approximator, then the resulting approximation-error-vector converges to zeros in a L 2 sense as d → ∞. Meanwhile, a new general integral inequality is derived which has an approximation-error-related term at its lower bound. By constructing a Liapunov-Krasovskii functional via the application of this inequality, sufficient conditions for the asymptotic (exponential) stability of the linear CDDS can be derived in terms of LMIs where the corresponding approximation error is expressed by a positive definite matrix. Our stability condition is further proved to have a hierarchical feasibility enhancement under the framework of orthogonal functions, which generalizes the results in [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF]. Finally, two numerical examples are given to demonstrate the effectiveness and capacity of the proposed methodologies.

The novelties in this paper are summarized as follows:

• We proposed a much wider class of functions for the approximation of delay kernel functions compared to polynomials approximation. The idea behind the proposed approximation is based on the application of the Hilbert projection theorem (least square approximation in Hilbert space) [START_REF] Muscat | Functional Analysis: An Introduction to Metric Spaces, Hilbert Spaces, and Banach Algebras[END_REF]. Moreover, this also allows one to construct a Liapunov-Krasovskii functional with non-polynomials kernels when distributed delay terms are handled via approximations. • A very general inequality is derived in this paper where the approximator functions belong to a weighted L 2 space. This means that the proposed integral inequality has tremendously general structure and it has very good potential to be used for the situations where a general approximator is employed to deal with distributed delay kernels. Thus many future works might be done based on the proposed inequality. Notations and rules: Empty matrices [], which follows the rules in Matlab environment, is applied in this paper. N stands for the set of natural numbers and N 0 := N ∪ {0}. We also define R ≥a := {x ∈ R :

x ≥ a} and S n = {X ∈ R n×n : X = X ⊤ }. M (L (X ) /B(R)) := f (•) ∈ R X : ∀Y ∈ B(R), f -1 (Y) ∈ L(X )
denotes the space of all L (X ) /B(R) measurable functions from X onto R, where L (X ) contains all the subsets of X which are Lebesgue measurable with X ∈ L (R), and B(R) is the Borel σ-algebra on R. For any X ∈ R n×n , we define Sy(X) := X + X ⊤ . A column vector containing a sequence of mathematical objects (scalars, matrices, functions etc) is defined as

Col n i=1 x i := Row n i=1 x ⊤ i ⊤ = x ⊤ 1 • • • x ⊤ i • • • x ⊤ n ⊤ . In addition, it is defined that Col n i=1 = [] with n < 1, where [] is an empty matrix with appropriate column dimensions. The symbol * is used to indicate [ * ]Y X = X ⊤ Y X or X ⊤ Y [ * ] = X ⊤ Y X or A B * C = A B B ⊤ C .
O n×m stands for a n × m zero matrix and the notation O n×n is abbreviated as O n , whereas 0 n denotes a n × 1 column vector. The diagonal sum of two matrices is defined as X ⊕ Y = X O * Y . ⊗ denotes the Kronecker product. Finally, we assume the operator precedence of matrices to be matrix (scalars) multiplications > ⊗ > ⊕ > +.

II. PRELIMINARIES AND PROBLEM FORMULATION

In this paper, we consider the stability analysis of a CDDS

ẋ(t) = A 1 x(t) + A 2 y(t -r) + 0 -r A 3 (τ )y(t + τ )dτ y(t) = A 4 x(t) + A 5 y(t -r) x(t 0 ) = ξ ∈ R n , y(t 0 + θ) = ϕ(θ), ∀θ ∈ [-r, 0] (1)
with a general distributed delay, where t 0 ∈ R and r > 0 is given and

ϕ(•) ∈ C ([-r, 0) R ν ). C([-r, 0) R n )
here stands for the space of the functions which are bounded, right-continuous and piecewise continuous functions, endowed with the norm ∥ϕ(•)∥ ∞ = sup τ ∈X ∥ϕ(τ )∥ 2 . Furthermore, x(t) ∈ R n ; y(t) ∈ R ν are the solution of (1) with n; ν ∈ N which determine the size of A 1 ; A 2 ; A 3 (τ ) and A 4 ; A 5 , where

A 3 (•) ∈ L 2 ([-r, 0] R n×ν ). Meanwhile, we define f (•) := Col d i=1 f i (•) ∈ C R R d with d ∈ N 0 and we assume f (•) satisfies the property ∃M ∈ R d×d , df (τ ) dτ = M f (τ ). (2)
It is obvious that given a

f (•) ∈ C R R d satisfying (2), then for any A 3 (•) ∈ L 2 ([-r, 0] R n×ν ) in (1) there exist A 3 ∈ R n×ρν with ρ = d + δ, and φ(•) ∈ L 2 [-r, 0] R δ with δ ∈ N 0 such that ∀τ ∈ [-r, 0], A 3 (τ ) = A 3 φ(τ ) f (τ ) ⊗ I ν (3) 0 -r φ(τ ) f (τ ) ⊤ φ ⊤ (τ ) f ⊤ (τ ) dτ ≻ 0 (4)
Remark 1. The condition in [START_REF] Rsvan | Functional differential equations of lossless propagation and almost linear behavior[END_REF] indicates that the functions in f (•) are the solutions of linear homogeneous equations with constant coefficients, such as polynomials, trigonometric and exponential functions etc. Moreover, the matrix inequality in (4) indicates that the functions in φ ⊤ (τ ) f ⊤ (τ ) ⊤ are linearly independent in a Lebesgue sense. (See Theorem 7.2.10 in [START_REF] Horn | Matrix Analysis[END_REF]) Finally, we stress here that the decomposition in (3) with ( 4) is always achievable since φ(

•) ∈ L 2 [-r, 0] R δ can cover any function in A 3 (•) ∈ L 2 ([-r, 0] R n×ν ) if such function is not covered by f (•).
Remark 2. The distributed delay A 3 (•) in ( 1) can be equivalently denoted by the decomposition in (3) which is employed in this paper so that well-posed stability condition can be derived. This will be illustrated later in light of the results in Lemma 3 and Theorem 1. Note that φ(•) can be a 0 × 1 empty vector if δ = 0. Finally, [START_REF] Briat | Linear Parameter Varying and Time-Delay Systems[END_REF] generalizes almost all the models of linear distributed delay (neutral) systems in the literature [START_REF] Fridman | H ∞ control of distributed and discrete delay systems via discretized Lyapunov functional[END_REF], [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF], [START_REF] Feng | Stabilization of uncertain linear distributed delay systems with dissipativity constraints[END_REF] without considering uncertainties or disturbances or inputs or outputs.

To tackle the function φ(•) ∈ L 2 [-r, 0] R δ in (3), which may not satisfy (2), we use f (τ ) to approximate φ(τ ) over [-r, 0]. This is denoted mathematically as

∀τ ∈ [-r, 0], φ(τ ) = Γ d f (τ ) + ε d (τ ) (5)
where

Γ d := 0 -r φ(τ )f ⊤ (τ )dτ F d , F -1 d := 0 -r f (τ )f ⊤ (τ )dτ (6)
and ε d (τ ) := φ(τ ) -Γ d f (τ ) defines the error of approximations which gives

E d := 0 -r ε d (τ )ε ⊤ d (τ )dτ ∈ S δ (7)
to measure the error residual of (5).

Remark 3. Note that using the approximation scheme in [START_REF]Reducing the complexity of the sum-of-squares test for stability of delayed linear systems[END_REF] does not put extra constraints to f (•) in ( 2)-(4) and f (•) can be any function as long as (2)-( 4) are satisfied. To obtain non-conservative results, we suggest that one can choose f (•) to be in line with the functions in A 3 (•) considering the structures of φ(τ ). Namely, f (•) should be "easy" to approximate φ(•) over [-r, 0] if δ ̸ = 0.

Remark 4. If f (•) contains only Legendre polynomials, then (5)-( 6) generalizes the polynomials approximation scenario in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF] via a matrix framework. The mathematical structure of (6) can be interpreted as a vector form of the Hilbert projection Theorem (Least square approximation) [START_REF] Muscat | Functional Analysis: An Introduction to Metric Spaces, Hilbert Spaces, and Banach Algebras[END_REF] via f (•) in (3). In addition, the approximator f (•) in ( 6) belongs to a class of elementary functions, not approximation by a specific type of functions such as polynomials. Thus the proposed approximation scheme in ( 5) is both conceptually and mathematically more general than the polynomials approximation in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF] as f (•) in [START_REF] Gu | Lyapunov Krasovskii functional for uniform stability of coupled differential-functional equations[END_REF] can contain functions which are not orthogonal polynomials or even mutually orthogonal.

The following lemma of the properties of the Kronecker product will be used throughout the rest of our paper. Lemma 1. For all X ∈ R n×m and Y ∈ R m×p and Z ∈ R q×r ,

(X ⊗ I q )(Y ⊗ Z) = (XY ) ⊗ (I q Z) = (XY ) ⊗ Z = (XY ) ⊗ (ZI r ) = (X ⊗ Z)(Y ⊗ I r ). (8)
Moreover, for all X ∈ R n×m , we have

A B C D ⊗ X = A ⊗ X B ⊗ X C ⊗ X D ⊗ X (9)
for any A, B, C, D with appropriate dimensions which make the block matrix at the left hand of the equality in [START_REF] Li | Discretized LKF method for stability of coupled differential-difference equations with multiple discrete and distributed delays[END_REF] to be compatible.

Considering ( 3) and ( 5), now (1) can be re-expressed as

ẋ(t) = Aϑ(t), y(t) = Ξϑ(t), Ξ = A 4 A 5 O ν×ρν x(t 0 ) = ξ ∈ R n , y(t 0 + θ) = ϕ(θ), ∀θ ∈ [-r, 0] (10)
where

A := A 1 A 2 A 3 Γ d I d ⊗ I ν A 3 E d O d×δ ⊗ I ν (11) ϑ(t) := Col x(t) y(t -r) , 0 -r F d (τ )y(t + τ )dτ 0 -r E d (τ )y(t + τ )dτ ∈ R κ (12)
with

κ := n + ν + ρν and E d (τ ) := E -1 d ε d (τ ) ⊗ I ν with ε d (τ
) and E d are defined in [START_REF] Pepe | On the asymptotic stability of coupled delay differential and continuous time difference equations[END_REF]. Note that E -1 d is well defined and we will explain it in Remark 6 in the next section. Note that also the matrix in [START_REF] Breda | Stability of linear delay differential equations: A numerical approach with MATLAB[END_REF] associated with the distributed delay terms are derived based on the identity

φ(τ ) f (τ ) ⊗ I ν y(t + τ ) = Γ d I d ⊗ I ν F d (τ )y(t + τ ) + E d O d×δ ⊗ I ν E d (τ )y(t + τ )
which itself is obtained via the property in [START_REF] Pepe | On the Liapunov-Krasovskii methodology for the ISS of systems described by coupled delay differential and difference equations[END_REF].

III. MAIN RESULTS ON STABILITY ANALYSIS

In this section, the main results on the stability analysis of (1) are presented in subsection III-B preceded by the presentation of some important Lemmas in Subsection III-A where we present a novel integral inequality and the Liapunov-Krasovskii stability criteria for [START_REF] Briat | Linear Parameter Varying and Time-Delay Systems[END_REF].

A. Mathematical prelimitaries

Lemma 2. Given r > 0, the trivial solution x(t) ≡ 0 n , y(t) ≡ 0 ν of (1) is globally uniformly asymptotically stable if there exist ϵ 1 ; ϵ 2 ; ϵ 3 > 0 and a differentiable functional v :

R n × C([-r, 0) R ν ) → R such that v(0 n , 0 ν (•)) = 0 and ϵ 1 ∥ξ∥ 2 2 ≤ v(ξ, ϕ(•)) ≤ ϵ 2 [max(∥ξ∥ 2 , ∥ϕ(•)∥ ∞ )] 2 (13) v(ξ, ϕ(•)) = d + dt v(x(t), y t (•)) t=t0,x(t0)=ξ,yt 0 (•)=ϕ(•) ≤ -ϵ 3 ∥ξ∥ 2 2 ( 14
)
for any ξ ∈ R n and ϕ(•) ∈ C([-r, 0) R ν ) in (1), where t 0 ∈ R and d + dx f (x) = limsup η↓0 f (x+η)-f (x) η
. Furthermore, y t (•) in ( 14) is defined by ∀t ≥ t 0 , ∀θ ∈ [-r, 0), y t (θ) = y(t + θ) where x(t) and y(t) satisfying [START_REF] Briat | Linear Parameter Varying and Time-Delay Systems[END_REF]. [START_REF] Gu | Lyapunov Krasovskii functional for uniform stability of coupled differential-functional equations[END_REF] be quadratic functions with the multiplier factors ϵ 1 ; ϵ 2 ; ϵ 3 > 0. Since (1) is a particular case of the general system considered in Theorem 3 of [START_REF] Gu | Lyapunov Krasovskii functional for uniform stability of coupled differential-functional equations[END_REF], then Lemma 2 is obtained. Now we derive the following new integral inequality which will be employed later in deriving our stability condition. First of all, we define the weighted Lebesgue function space

Proof: Let u(•), v(•), w(•) in Theorem 3 of
L 2 ϖ K R d := ϕ(•) ∈ M L (K) /B(R d ) : ∥ϕ(•)∥ 2,ϖ < ∞ (15) with d ∈ N 0 and ∥ϕ(•)∥ 2,ϖ := K ϖ(τ )ϕ ⊤ (τ )ϕ(τ )dτ where ϖ(•) ∈ M (L (K) /B(R ≥0
)) and the function ϖ(•) has only countably infinite or finite numbers of zero values. Furthermore, K ⊆ R and the Lebesgue measure of K is non-zero.

Lemma 3. Given K and ϖ(•) in [START_REF] Münz | Robust stabilization and H ∞ control of uncertain distributed delay systems[END_REF] and

U ∈ S n ⪰0 := {X ∈ S n : X ⪰ 0} with n ∈ N. Let f(•) := Col d i=1 f i (•) ∈ L 2 ϖ K R d and g(•) := Col δ i=1 g i (•) ∈ L 2 ϖ K R δ with d ∈ N and δ ∈ N 0 , in which the functions f(•) and g(•) satisfy K ϖ(τ ) g(τ ) f(τ ) g ⊤ (τ ) f ⊤ (τ ) dτ ≻ 0. ( 16
)
Then the inequality

K ϖ(τ )x ⊤ (τ )U x(τ )dτ ≥ [ * ] (F d ⊗ U ) K ϖ(τ )F(τ )x(τ )dτ + [ * ] E -1 d ⊗ U K ϖ(τ )E(τ )x(τ )dτ ( 17
)
holds for all x(•) ∈ L 2 ϖ (K R n ), where F(τ ) = f(τ ) ⊗ I n ∈ R dn×n and F -1 d = K ϖ(τ )f(τ )f ⊤ (τ )dτ ∈ S d ≻0 . In addition, E(τ ) = e(τ ) ⊗ I n ∈ R δn×n where e(τ ) = g(τ ) -Af(τ ) ∈ R δ and A = K ϖ(τ )g(τ )f ⊤ (τ )dτ F d ∈ R δ×d and E d := K ϖ(τ )e(τ )e ⊤ (τ )dτ ∈ S δ .
Proof: The proof of Lemma 3 is inspired by the proofs of Lemma 2 in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF] and the proof of Lemma 5 in [START_REF] Feng | Stabilization of uncertain linear distributed delay systems with dissipativity constraints[END_REF]. Firstly, we have

E d := K ϖ(τ )e(τ )e ⊤ (τ )dτ = I δ -A K ϖ(τ ) g(τ ) f(τ ) g ⊤ (τ ) f ⊤ (τ ) dτ I δ -A ⊤ ≻ 0, (18)
where the positive definiteness in ( 18) can be established based on [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF] and the property of congruence transformations with the fact that rank

I δ -A = δ. Consequently, E -1 d is well defined. Let y(τ ) := x(τ ) -F ⊤ (τ )(F d ⊗ I n ) K ϖ(θ)F(θ)x(θ)dθ -E ⊤ (τ ) E -1 d ⊗ I n K ϖ(θ)E(θ)x(θ)dθ,
where the definitions of F(•) and E(•) have been given in Lemma 3.

By A = K ϖ(τ )g(τ )f ⊤ (τ )dτ F d and e(τ ) = g(τ ) -Af(τ ) ∈ R δ , we have K ϖ(τ )e(τ )f ⊤ (τ )dτ = K ϖ(τ ) [g(τ ) -Af(τ )] f ⊤ (τ )dτ = K ϖ(τ )g(τ )f ⊤ (τ )dτ -A K ϖ(τ )f(τ )f ⊤ (τ )dτ = K ϖ(τ )g(τ )f ⊤ (τ )dτ - K ϖ(τ )g(τ )f ⊤ (τ )dτ F d F -1 d = O δ×d . ( 19
)
Substituting the expression of y(•) into K ϖ(τ )y ⊤ (τ )U y(τ )dτ and considering [START_REF] Park | Auxiliary function-based integral inequalities for quadratic functions and their applications to time-delay systems[END_REF] yields

K ϖ(τ )y ⊤ (τ )U y(τ )dτ = K ϖ(τ )x ⊤ (τ )U x(τ )dτ -2 K ϖ(τ )x ⊤ (τ )U F ⊤ (τ )dτ (F d ⊗ I n )ζ + ζ ⊤ K ϖ(τ )(F d ⊗ I n ) ⊤ F(τ )U F ⊤ (τ )(F d ⊗ I n )dτ ζ -2 K ϖ(τ )x ⊤ (τ )U E ⊤ (τ )dτ (E -1 d ⊗ I n )ω + ω ⊤ K ϖ(τ )(E -1 d ⊗ I n ) ⊤ E(τ )U E ⊤ (τ )(E -1 d ⊗ I n )dτ ω ( 20
)
where ζ := K ϖ(θ)F(θ)x(θ)dθ and ω := K ϖ(θ)E(θ)x(θ)dθ. Now apply [START_REF] Pepe | On the Liapunov-Krasovskii methodology for the ISS of systems described by coupled delay differential and difference equations[END_REF] to the term U F ⊤ (τ ) and U E ⊤ (τ ) and consider F(τ ) = f(τ ) ⊗ I n and E(τ ) = e(τ ) ⊗ I n . Then we have

U F ⊤ (τ ) = F ⊤ (τ )(I d ⊗ U ), U E ⊤ (τ ) = E ⊤ (τ )(I δ ⊗ U ) (21)
given [START_REF] Chen | Two general integral inequalities and their applications to stability analysis for systems with time-varying delay[END_REF] with [START_REF] Pepe | On the Liapunov-Krasovskii methodology for the ISS of systems described by coupled delay differential and difference equations[END_REF] to the integral terms in (20) related to [START_REF] Chen | Two general integral inequalities and their applications to stability analysis for systems with time-varying delay[END_REF] yields

(X ⊗ Y ) ⊤ = X ⊤ ⊗ Y ⊤ . Now applying
K ϖ(τ )x ⊤ (τ )U F ⊤ (τ )dτ (F d ⊗ I n )ζ = ζ ⊤ (F d ⊗ U )ζ K ϖ(τ )x ⊤ (τ )U E ⊤ (τ )dτ (E -1 d ⊗ I n )ω = ω ⊤ (E -1 d ⊗ U )ω, (22)
K (F d ⊗ I n ) ⊤ ϖ(τ )F(τ )U F ⊤ (τ )(F d ⊗ I n )dτ = (F d ⊗ I n ) K ϖ(τ )F(τ )F ⊤ (τ )dτ (F d ⊗ U ) (23)
K (E -1 d ⊗ I n ) ⊤ ϖ(τ )E(τ )U E ⊤ (τ )(E -1 d ⊗ I n )dτ = (E -1 d ⊗ I n ) K ϖ(τ )E(τ )E ⊤ (τ )dτ E -1 d ⊗ U . ( 24
)
Meanwhile, since F(τ ) = f(τ ) ⊗ I n and E(τ ) = e(τ ) ⊗ I n , we have

K ϖ(τ )F(τ )F ⊤ (τ )dτ = K ϖ(τ )f(τ )f ⊤ (τ )dτ ⊗ I n = F -1 d ⊗ I n ( 25
) K ϖ(τ )E(τ )E ⊤ (τ )dτ = K ϖ(τ )e(τ )e ⊤ (τ )dτ ⊗ I n = E d ⊗ I n . ( 26
)
Moreover, ( 23)-( 24) can be reformulated into

K (F d ⊗ I n ) ⊤ ϖ(τ )F(τ )U F ⊤ (τ )(F d ⊗ I n )dτ = (F d ⊗ I n ) K ϖ(τ )F(τ )F ⊤ (τ )dτ (F d ⊗ U ) = F d ⊗ U, K (E -1 d ⊗ I n ) ⊤ ϖ(τ )E(τ )U E ⊤ (τ )(E -1 d ⊗ I n )dτ = (E -1 d ⊗ I n ) K ϖ(τ )E(τ )E ⊤ (τ )dτ E -1 d ⊗ U = E -1 d ⊗ U (27)
by ( 25) and ( 26). Now by ( 27) and ( 22), ( 20) can be simplified into

K ϖ(τ )y ⊤ (τ )U y(τ )dτ = K ϖ(τ )x ⊤ (τ )U x(τ )dτ -[ * ] (F d ⊗ U ) K ϖ(τ )F(τ )x(τ )dτ -[ * ] E -1 d ⊗ U K ϖ(τ )E(τ )x(τ )dτ . ( 28
)
Given U ⪰ 0, one can obtain ( 17) via (28).

Remark 5. By Theorem 7.2.10 in [START_REF] Horn | Matrix Analysis[END_REF], we know [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF] indicates that the functions in f(•) and g(•) are linearly independent in a Lebesgue sense. Since

f(•) ∈ L 2 ϖ K R d and g(•) ∈ L 2 ϖ K R δ with plenty options for ϖ(•) ∈ L ∫ (K R ≥0
), then it is certain that the structure of the inequality ( 17) is tremendously general. If f(•) contains only Legendre polynomials, then the result of Lemma 1 in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF] is covered by [START_REF] Muscat | Functional Analysis: An Introduction to Metric Spaces, Hilbert Spaces, and Banach Algebras[END_REF]. With δ = 0, then one can conclude that the inequalities in [START_REF] Park | Auxiliary function-based integral inequalities for quadratic functions and their applications to time-delay systems[END_REF], [START_REF] Gyurkovics | Multiple integral inequalities and stability analysis of time delay systems[END_REF], [START_REF] Chen | Two general integral inequalities and their applications to stability analysis for systems with time-varying delay[END_REF], [START_REF] Liu | Generalized Jensen inequalities with application to stability analysis of systems with distributed delays over infinite time-horizons[END_REF], [START_REF] Feng | Stabilization of uncertain linear distributed delay systems with dissipativity constraints[END_REF] are generalized by [START_REF] Muscat | Functional Analysis: An Introduction to Metric Spaces, Hilbert Spaces, and Banach Algebras[END_REF] with appropriate ϖ(•), f(•) and x(•). Note that we assume the inverse of a 0 × 0 matrix is still a 0 × 0 matrix which implies that E -1 d becomes a 0 × 0 matrix with δ = 0. Remark 6. In [START_REF] Muscat | Functional Analysis: An Introduction to Metric Spaces, Hilbert Spaces, and Banach Algebras[END_REF], f(•) can be interpreted as to approximate g(•). By letting f(τ ) = f (τ ) and g(τ ) = φ(τ ) with ϖ(τ ) = 1 in Lemma 3, then we have E d = E d where E d is given in [START_REF] Pepe | On the asymptotic stability of coupled delay differential and continuous time difference equations[END_REF]. This also indicates that E d is invertible. Furthermore, if f(•) contains only orthogonal functions, then the behavior of E d can be quantitatively characterized by the following corollary. Finally, since f(•) and g(•) in [START_REF] Muscat | Functional Analysis: An Introduction to Metric Spaces, Hilbert Spaces, and Banach Algebras[END_REF] belong to very general function space and ϖ(•) can be adjusted based on specific requirements, thus the results in Lemma 3 is much more superior both conceptually and mathematically than the result of Lemma 1 in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF].

Corollary 1. Given all the parameters defined in Lemma 3 and assume that

f(•) = Col d i=1 f i (•) where {f i (•)} ∞ i=1
contains only orthogonal functions. Then we have that for all d ∈ N 0

0 ≺ E d+1 = E d - K ϖ(τ )f 2 d+1 (τ )dτ a d+1 a ⊤ d+1 ⪯ E d ( 29
)
where

a d+1 := K ϖ(τ )g(τ )f d+1 (τ )dτ K ϖ(τ )f 2 d+1 (τ )dτ -1 ∈ R δ and f d+1 (•) ∈ L 2 ϖ (K R)
. Proof: Note that only the dimension of f(•) is related to d, whereas δ as the dimension of g(•) is independent from d. Since f(•) contains only orthogonal functions, it is obvious that

F d+1 = F d ⊕ K ϖ(τ )f 2
d+1 (τ )dτ -1 (See the Definition 1 in [START_REF] Feng | Orthogonal functions based integral inequalities and their applications to time delay systems[END_REF]). By using this property considering the definition of e(τ ) in Lemma 3, it follows that

e d+1 (τ ) = g(τ ) - K ϖ(τ )g(τ ) f ⊤ (τ ) f d+1 (τ ) dτ F d ⊕ K ϖ(τ )f 2 d+1 (τ )dτ -1 f(τ ) f d+1 (τ ) = g(τ ) -A d a d+1 f(τ ) f d+1 (τ ) = e d (τ ) -f d+1 (τ )a d+1 (30)
for all d ∈ N 0 , where a d+1 has been defined in ( 29) and e d (τ ) = g(τ ) -A d f(τ ). Note that the index d is added to the symbols A and e(τ ) in Lemma 3 without causing ambiguity. Considering (30) and ( 18), we have

0 ≺ E d+1 = K ϖ(τ )e d+1 (τ )e ⊤ d+1 (τ )dτ = E d -Sy a d+1 K ϖ(τ )f d+1 (τ )e ⊤ d (τ )dτ + K ϖ(τ )f 2 d+1 (τ )dτ a d+1 a ⊤ d+1 . (31)
By [START_REF] Park | Auxiliary function-based integral inequalities for quadratic functions and their applications to time-delay systems[END_REF] and the fact that K ϖ(τ )f d+1 (τ )f(τ )dτ = 0 d due to the orthogonality among {f i (•)} ∞ i=1 , we have

O δ×(d+1) = K ϖ(τ )e d+1 (τ ) f ⊤ (τ ) f d+1 (τ ) dτ = K ϖ(τ ) (e d (τ ) -a d+1 f d+1 (τ )) f ⊤ (τ ) f d+1 (τ ) dτ = K ϖ(τ ) e d (τ )f ⊤ (τ ) f d+1 (τ )e d (τ ) dτ -a d+1 K ϖ(τ ) f d+1 (τ )f ⊤ (τ ) f 2 d+1 (τ ) dτ = O δ×d K ϖ(τ )f d+1 (τ )e d (τ )dτ -O δ×d K ϖ(τ )f 2 d+1 (τ )dτ a d+1 = O δ×(d+1) . (32) which leads to the conclusion that K ϖ(τ )f d+1 (τ )e d (τ )dτ = K ϖ(τ )f 2 d+1 (τ )dτ a d+1 .
Substituting this equality into (31) yields (29) given K ϖ(τ )f 2 d+1 (τ )dτ > 0 and a d+1 a ⊤ d+1 ⪰ 0. Remark 7. The conclusion of Lemma 1 in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF] is generalized by Corollary 1 as f(•) can have Legendre polynomials.

B. Stability Analysis

The main result on the stability analysis of (1) is presented in the following theorem.

Theorem 1. Given (1) with r > 0 and f

(•) ∈ C R R d ,A 3 ∈ R n×ρν , φ(•) ∈ L 2 [-r, 0] R δ in (3) with d ∈ N
and δ ∈ N 0 and the matrices Γ d and E d in (5)- [START_REF] Pepe | On the asymptotic stability of coupled delay differential and continuous time difference equations[END_REF], then the trivial solution of (1) is globally asymptotically stable if there exist P ∈ S n+dν and Q; R ∈ S ν such that

P := P + O n ⊕ F d ⊗ Q ≻ 0 (33) Q ≻ 0, R ≻ 0, Φ ≺ 0 ( 34
)
where

F -1 d = 0 -r f (τ )f ⊤ (τ )dτ and Φ := Sy Θ ⊤ 2 P Θ 1 -O n ⊕ Q ⊕ F d ⊗ R ⊕ E d ⊗ R + Ξ ⊤ (Q + rR)Ξ (35) Θ 1 := Col A, (H ⊗ I ν ) Π , Θ 2 = Υ O (n+dν)×δν (36)
H = f (0) -f (-r) -M O d×δ (37) Υ := I n O n×ν O n×dν O dν×n O dν×ν I dν (38) Π = Col Ξ, O (ν+dν+δν)×n I ν+dν+δν (39)
with A in [START_REF] Breda | Stability of linear delay differential equations: A numerical approach with MATLAB[END_REF] and Ξ in [START_REF] Vyhlídal | QPmR -Quasi-Polynomial Root-Finder: Algorithm Update and Examples[END_REF] and M in [START_REF] Rsvan | Functional differential equations of lossless propagation and almost linear behavior[END_REF].

Proof: To analyze the stability of (1), we consider the Krasovskii functional

v(ξ, ϕ(•)) = η ⊤ P η + 0 -r ϕ ⊤ (τ ) Q + (τ + r)R ϕ(τ )dτ (40)
to be constructed, where P ∈ S n+dν and Q, R ∈ S ν and

η := Col ξ, 0 -r F d (τ )ϕ(τ )dτ (41) with F d (τ ) = f (τ ) ⊗ I ν ∈ R dν×ν .
Since (1) can be equivalently denoted by [START_REF] Vyhlídal | QPmR -Quasi-Polynomial Root-Finder: Algorithm Update and Examples[END_REF] via f (•), A 3 , φ(•), differentiate v(x(t), y t (•)) with x(t), y t (•) in Lemma 2 at t = t 0 in light of [START_REF] Vyhlídal | QPmR -Quasi-Polynomial Root-Finder: Algorithm Update and Examples[END_REF], it produces

d + dt v(x(t), y t (•)) t=t0,x(t0)=ξ,yt 0 (•)=ϕ(•) = ϑ ⊤ (t 0 ) Sy Θ ⊤ 2 P Θ 1 ϑ(t 0 ) -ϕ ⊤ (-r)Qϕ(-r) + ϕ ⊤ (0) (Q + rR) ϕ(0) - 0 -r ϕ ⊤ (τ )Rϕ(τ )dτ (42)
where ϑ(t) and Θ 1 ;Θ 2 have been defined in ( 12) and (36), respectively, and

ϑ(t 0 ) := Col ξ ϕ(-r) , 0 -r F d (τ )ϕ(τ )dτ 0 -r E d (τ )ϕ(τ )dτ ∈ R κ ( 43
)
given the initial conditions in [START_REF] Vyhlídal | QPmR -Quasi-Polynomial Root-Finder: Algorithm Update and Examples[END_REF]. Note that the relation

0 -r F d (τ ) ẏ(t + τ )dτ t=t0 = F d (0)ϕ(0) -F d (-r)ϕ(-r) -[M ⊗ I ν ] 0 -r F d (τ )ϕ(τ )dτ = (H ⊗ I ν ) Πϑ(t 0 ), (44)
has been applied to derive H in (37), where (44) itself can be derived via ( 1)-( 2) with ( 8)-( 9). To obtain a upper bound for (42

), let R ≻ 0, ϖ(τ ) = 1, f(τ ) = f (•) and g(τ ) = φ(τ ), then one can derive 0 -r ϕ ⊤ (τ )Rϕ(τ )dτ ≥ [ * ] (F d ⊗ R) 0 -r F d (τ )ϕ(τ )dτ + [ * ] (E d ⊗ R) 0 -r E d (τ )ϕ(τ )dτ (45)
via [START_REF] Muscat | Functional Analysis: An Introduction to Metric Spaces, Hilbert Spaces, and Banach Algebras[END_REF] given ( 4) and the forms in ( 5)- [START_REF] Pepe | On the asymptotic stability of coupled delay differential and continuous time difference equations[END_REF], where F d (τ ) and E d (τ ) are given in [START_REF] Fridman | H ∞ control of distributed and discrete delay systems via discretized Lyapunov functional[END_REF]. Now applying (45) with R ≻ 0 to (42) yields

d + dt v(x(t), y t (•)) t=t0,x(t0)=ξ,yt 0 (•)=ϕ(•) ≤ ϑ ⊤ (t 0 )Φϑ(t 0 ) (46)
for all t 0 ∈ R where Φ is given in (35). It is obvious that there exists ϵ 3 > 0 such that

d + dt v(x(t), y t (•)) t=t0,x(t0)=ξ,yt 0 (•)=ϕ(•) ≤ -ϵ 3 ∥ξ∥ 2 (47) if Φ ≺ 0.
Thus it follows that if Φ ≺ 0 and R ≻ 0 are satisfied then (40) satisfies [START_REF] Feng | Stabilization of uncertain linear distributed delay systems with dissipativity constraints[END_REF]. Hence the feasible solutions of (34) infers the existence of (40) satisfying ( 14). Now we start to prove that (33) with Q ≻ 0 and R ≻ 0 in (34) infers that (40) satisfies [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF]. First of all, given the structure of (40), it follows that ∃λ; η > 0 :

v(ξ, ϕ(•)) ≤ ξ 0 -r F (τ )ϕ(τ )dτ ⊤ λ ξ 0 -r F (τ )ϕ(τ )dτ + 0 -r ϕ ⊤ (τ )λϕ(τ )dτ ≤ λ∥ξ∥ 2 2 + 0 -r ϕ ⊤ (τ )F ⊤ (τ )dτ λ 0 -r F (τ )ϕ(τ )dτ + λr∥ϕ(•)∥ 2 ∞ ≤ λ∥ξ∥ 2 2 + λr∥ϕ(•)∥ 2 ∞ + 0 -r F (τ )ϕ(τ )dτ ⊤ (ηF ⊗ I n ) 0 -r F (τ )ϕ(τ )dτ ≤ λ∥ξ∥ 2 2 + λr∥ϕ(•)∥ 2 ∞ + 0 -r ϕ ⊤ (τ )ηϕ(τ )dτ ≤ λ∥ξ∥ 2 2 + (λr + ηr) ∥ϕ(•)∥ 2 ∞ ≤ (λ + λr + ηr) ∥ξ∥ 2 2 + ∥ϕ(•)∥ 2 ∞ ≤ 2 (λ + λr + ηr) [max (∥ξ∥ 2 , ∥ϕ(•)∥ ∞ )] 2 (48)
which is derived via the property of quadratic forms: ∀X ∈ S n , ∃λ > 0 : ∀x ∈ R n \ {0}, x ⊤ (λI n -X) x > 0 with [START_REF] Muscat | Functional Analysis: An Introduction to Metric Spaces, Hilbert Spaces, and Banach Algebras[END_REF] with f(τ ) = f (τ ) and δ = 0. This shows that (40) satisfies the upper bound property in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF]. Now to construct a lower bound for v(ξ, ϕ(•)), apply ( 17) to (40) with ϖ(τ ) = 1, δ = 0 and f(τ ) = f (τ ) matching the term F d (τ ) in (40), which gives

0 -r ϕ ⊤ (τ )Qϕ(τ )dτ ≥ [ * ] F d ⊗ Q 0 -r F d (τ )ϕ(τ )dτ. ( 49
)
Considering the conclusion in (48) and by further applying (49) to (40) with Q ≻ 0 and R ≻ 0 in (34), it yields that ( 13) is satisfied if (33) and Q ≻ 0 and R ≻ 0 hold. This shows that feasible solutions of (33)-(34) infers the existence of (40) satisfying all the conditions in ( 13) and ( 14).

Remark 8. If one wants to increase the dimension of f (•) in (40), then extra zeros need to be added to A 3 in (3) to ensure compatibility. Note that there is no upper bound imposed on the dimension of f (•). Furthermore, if (1) is a retarded system, then the functional in (40) generalizes the Krasovskii functional in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF] where f (•) contains only Legendre polynomials. Finally, the use of f (•) in ( 2) to approximate φ(•) allows one to construct functional with more sophisticated kernels in (40) other than only the option of polynomials for (40).

Remark 9. To apply ( 17) at (45) with ϖ(τ ) = 1 and f(τ ) = f (τ ) and g(τ ) = φ(τ ), the linear independent condition (4) must be satisfied which guarantees an invertible E d . Moreover, the absence of φ(•) in (40) enables us to avoid dealing with φ(•). These are the major motivations for using the structure of the decomposition in (3) to handle distributed delay terms in this paper. 3) and (40) contains only Legendre polynomials, then Theorem 1 with the approximation scheme in (6) cover the stability results in Theorem 4 of [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF]. Note that one can apply congruence transformations to Φ to avoid potential numerical problems caused by the position of the error matrix E d in (34) if some of its eigenvalues are too small. Namely, the matrix inequality

Remark 10. If f (•) in (
Φ = I n+ν+ϱ ⊕ ηE -1 2 d ⊗ I ν ⊤ Φ I n+ν+ϱ ⊕ ηE -1 2 d ⊗ I ν ≺ 0 (50)
is equivalent to (34), where η ∈ R and the block-diagonals of Φ are no longer related to E d .

The following theorem shows that the LMI conditions in Theorem 1 exhibit a hierarchy with respect to d when f (•) contains only orthogonal functions.

Theorem 2. Let f (•) := Col d i=1 f i (•) ∈ C R R d in (2)
where {f i (•)} ∞ i=1 contains only functions which are mutually orthogonal over [-r, 0]. (See the Definition in [START_REF] Feng | Orthogonal functions based integral inequalities and their applications to time delay systems[END_REF] ). Given Γ d , E d (5)-( 7), we have

∀d ∈ N 0 , F d ⊆ F d+1 ( 51
)
where

F d := r > 0 | (33)-(34) hold & P ∈ S n+dν , Q; R ∈ S ν .
Proof: The proof here is inspired by the proof of Theorem 8 in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF]. Note that P ≻ 0 in (33) and Φ ≻ 0 in (34) are indexed by the value of d. Hence d might be automatically attached to the related variables throughout the entire proof. Let r ∈ F d with F d ̸ = ∅ which infers that there exist feasible solutions for (33) and (34) at d. Moreover, let P d ∈ S n+dν and Q; R ∈ S ν to be a feasible solution of (33)-(34). We will show that the feasible solutions of (33)-(34) exist at d + 1 if the feasible solutions of (33) and (34) at d exist.

Assume that

P d+1 := P d ⊕ O ν , η d+1 (t) := η d (t) 0 -r f d+1 (τ )y(t + τ )dτ . ( 52
)
with reference to the structures in (40). Since {f i (•)} ∞ i=1 contains functions which are orthogonal over [-r, 0], we have

F d+1 = F d ⊕ f d+1 with f -1 d+1 = 0 -r f 2 d+1 (τ )
dτ . Now by ( 52) and ( 9) with F d+1 = F d ⊕ f d+1 , it follows that P d+1 = P d ⊕ f d+1 Q with Q ≻ 0 and f d+1 > 0. Thus the feasible solutions of P d ≻ 0 infer the existence of the feasible solutions of P d+1 ≻ 0. Now we start to prove the hierarchical results for Φ d ≺ 0 in (34) with respect to d. Note that Φ d ≺ 0 can be written as

Φ d = Φd T d E d ⊗ I ν * -E d ⊗ R ≺ 0 (53)
where Φd can be easily obtained based on the structure of Φ d , and

T d = Υ ⊤ d P d I n O (ν+dν)×n A 3 I δ O d×δ ⊗ I ν ( 54
)
with Υ d = Υ given in (38) at d. Now apply Schur complement to (53) considering the fact that E d ≻ 0 and R ≻ 0. Then it yields

Φd + T d E d+1 + f -1 d+1 γ d+1 γ ⊤ d+1 ⊗ R -1 T ⊤ d = Φd + T d [ * ] f -1 d+1 R -1 O ν×νδ * E -1 d+1 ⊗ R -1 γ ⊤ d+1 ⊗ I ν E d+1 ⊗ I ν T ⊤ d ≺ 0 (55)
based on ( 29) and ( 8), where γ d+1 = 0 -r φ(τ )f d+1 (τ )dτ f d+1 . Apply Schur complement again to (55) yields that given R ≻ 0, (55) holds if and only if

Θ =   Φd T d γ d+1 ⊗ I ν T d E d+1 ⊗ I ν * -f d+1 R O ν×νδ * * -E d+1 ⊗ R   ≺ 0. ( 56
)
Now the fact is

Φ d+1 = Θ = Φd+1 T d+1 E d+1 ⊗ I ν * -E d+1 ⊗ R ≺ 0 (57)
given the structure of Φ d at d in (53) and the relations

A 3 O n×ν Γ d+1 I d+1 ⊗ I ν = A 3 O n×ν     Γ d+1 ⊗ I ν γ d+1 ⊗ I ν I dν O dν×ν O ν×dν I ν     = A 3 Γ d ⊗ I ν I dν A 3 γ d+1 ⊗ I ν O dν×ν = A 3 Γ d ⊗ I ν I dν A 3 I δν O dν×δν (γ d+1 ⊗ I ν ) , ( 58
)
T d+1 = (Υ d ⊕ I ν ) ⊤ (P d ⊕ O ν ) I n O (2ν+dν)×n A 3 O n×ν I δ O (d+1)×δ ⊗ I ν = T d O ν×δν (59)
in view of (54), and the coefficient A 5 O n×ν is obtained based on the relation

A 3 (τ ) = A 3 O n×ν φ ⊤ (τ ) f ⊤ (τ ) f d+1 (τ ) ⊗ I ν . ( 60
)
Hence we have proved that given Q ≻ 0 and R ≻ 0 in (34), the feasible solutions of Φ d ≺ 0 infers the existence of the feasible solutions of Φ d+1 ≺ 0.

Remark 11. Theorem 2 generalizes the result of Theorem 8 in [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF] as the latter one only considers f (•) to contain Legendre polynomials. [START_REF] Briat | Linear Parameter Varying and Time-Delay Systems[END_REF] IV. NUMERICAL SIMULATION In this section, two numerical examples are presented to show the effectiveness of our proposed method. Numerical calculations are conducted in Matlab environment with Yalmip [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF] and SDPT3 [START_REF] Toh | On the implementation and usage of SDPT3-a Matlab software package for semidefinite-quadratic-linear programming, version 4.0[END_REF].

Consider the parameters in the following table for the CDDS system in (1). To the best of our knowledge, no existing methods may analyze the stability of (1) with the parameters in Example 1 and 2. This is even true for the frequency domain approaches in [START_REF] Vyhlídal | QPmR -Quasi-Polynomial Root-Finder: Algorithm Update and Examples[END_REF], [START_REF] Breda | Stability of linear delay differential equations: A numerical approach with MATLAB[END_REF] due to the CDDS structures and non-trivial distributed delay terms in Table I.

In order to demonstrate the impact of having different f (•) for Theorem 1, we consider the following functions:

ℓ d (τ ) = Col d i=0 ℓ i (τ ) ∈ R d+1 and h d (τ ) = Col 1, Col d/2 i=1 sin 12iτ, Col d/2 i=1 cos 12iτ ∈ R d+1
= 0 ⊕ Od/2 ⊕ d/2 i=1 12i - ⊕ d/2 i=1 12i Od/2
in (2) for f (τ ) = h d (τ ), and d must be a positive even number for h d (τ ). Note that also the functions in h d (τ ) are not orthogonal over [-r, 0]. Now apply Theorem 1 with (50) and η = 1 and a testing vector r = (500 : 2500)/1000 to (1) with the parameters in Table I and Remark 12. Note that when one applies Theorem 1, the numerical results produced by f (τ ) = h d (τ ) and f (τ ) = ℓ d (τ ) are not mutually exclusive. Namely, users can choose the best results considering both side (and potentially even more if one wants to use more options for f (τ )) as the final conclusion produced by Theorem 1. This clearly shows one of the advantages of the proposed method in this paper.

A 3 = -5 0 ⊤ d+1 (62) A 3 = 0 0.02 0 0 -4 0 0 0 0 0 0.01 0 0 0 0 -10 O 2×(2d+2) , ( 63
For Example 1 and 2 in Table I, we did not run our program with higher d > 24 for f (τ ) = ℓ d (τ ) via the testing vector r = (500 : 2500)/1000. This is because the expected computing time becomes too long due to the calculation of Γ d , E d in ( 5)-( 7) via the function vpaintegral with variable precision in Matlab. On the other hand, the numerical integration function integral in Matlab is not an alternative option in this case to calculate Γ d and E d because of its limited numerical accuracy. The results in Tables II-V can be explained by the fact that the functions in φ(τ ) over [-r, 0] is not "easy" to be approximated by polynomials if r becomes relatively large. Thus the numerical results in this section give examples demonstrating the advantage of our proposed approximation scheme over the existing polynomials only approximation approach [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF].

V. CONCLUSION

In this paper, a new method for the stability analysis of a linear CDDS (1) with a general distributed delay has been proposed. The kernel of A 3 (•) can be any L 2 function over the delay interval, which has been handled by the decomposition (3) with a general form of approximation (5) using a class of elementary functions. The structure of the proposed approximation in ( 5)-( 7) is based on the idea of least square approximation in Hilbert Space. Thanks to the newly proposed inequality, the approximation error has been included in the proposed stability condition. It is also proved that the stability condition possesses a hierarchy when ( 6) is utilized with a f (•) containing functions which are mutually orthogonal. Two numerical examples have been presented which have shown the strength and effectiveness of our proposed methodology.

  , where ℓ d(τ ) contains Legendre polynomials ℓ i (τ )ℓ d (τ ) := F d = r -1 d i=0 2i +1 and the corresponding M in (2) for f (τ ) = ℓ d (τ ) can be easily determined. Note that h d (τ ) corresponds to M

  ) respectively, where (62) corresponds to both f (τ ) = h d (τ ) and f (τ ) = ℓ d (τ ) for Example 1 with φ(τ ) = sin(cos(12τ )), and (63) corresponds to both f (τ ) = h d (τ ) and f (τ ) = ℓ d (τ ) with φ(τ ) = e sin(12τ ) e cos(12τ ) sin(sin(12τ )) sin(cos(12τ )) ⊤ for Example 2, respectively. The computing results of detectable delay margins are summarized in

TABLE I PARAMETERS

 I OF NUMERICAL EXAMPLES FOR

Table II -

 II V where NoV denotes the number of decision variables required by each optimization program for a pointwise value of r.

	f (τ )	First interval Second interval NoV
	h d (τ ), d = 6	[0.599, 0.71]	[1.123, 1.233]	38
	h d (τ ), d = 10	[0.599, 0.71]	[1.122, 1.234]	80
	ℓ d (τ ), d = 10	[0.599, 0.71]	[1.126, 1.237]	80
	ℓ d (τ ), d = 23	[0.599, 0.71]	[1.122, 1.234]	327
				TABLE II
	DETECTABLE STABLE DELAY MARGINS OF EXAMPLE 1
	f (τ )	Third Interval Fourth Interval NoV
	h d (τ ), d = 6	[1.647, 1.757]	[2.171, 2.28]	38
	h d (τ ), d = 10	[1.646, 1.758]	[2.168, 2.281]	80
	ℓ d (τ ), d = 10 not detectable	not detectable	80
	ℓ d (τ ), d = 23	[1.661, 1.739]	not detectable	327
				TABLE III
	DETECTABLE STABLE DELAY MARGINS OF EXAMPLE 1
	f (τ )	First interval Second interval NoV
	h d (τ ), d = 6	[0.535, 0.774]	[1.059, 1.297]	142
	h d (τ ), d = 10 [0.535, 0.774]	[1.059, 1.297]	306
	ℓ d (τ ), d = 10 [0.535, 0.773]	[1.062, 1.287]	306
	ℓ d (τ ), d = 24 [0.535, 0.774]	[1.058, 1.298]	1384
				TABLE IV
	DETECTABLE STABLE DELAY MARGINS OF EXAMPLE 2
	f (τ )	Third Interval Fourth Interval	Fifth Interval	NoV
	h d (τ ), d = 6	[1.583, 1.82]	[2.107, 2.168]	[2.283, 2.344]	142
	h d (τ ), d = 10	[1.582, 1.821]	[2.106, 2.168]	[2.283, 2.344]	306
	ℓ d (τ ), d = 10 not detectable	not detectable	not detectable	306
	ℓ d (τ ), d = 24	[1.59, 1.69] ∪ [1.717, 1.809]	not detectable	not detectable 1384

TABLE V DETECTABLE

 V STABLE DELAY MARGINS OF EXAMPLE 2

ACKNOWLEDGEMENT

The authors would like to thank Prof. Keqin Gu for his comments on the choice of the function space for the states of CDDS. In addition, we thank Prof. Johan Löfberg and Prof. Dimitri Breda for the remarks concerning their Matlab software packages. Finally, special thanks to Prof. Ji Li Feng who gave us permission to use the high-speed computer in his lab.