Abdouramane Moussa Ali

Cédric Join
email: cjoin@cran.uhp-nancy.fr

Frédéric Hamelin
email: fhamelin@cran.uhp-nancy.fr

Fault diagnosis of uncertain

Keywords: distribution theory, spectra analysis, actuator fault, sensor fault, FDI

published or not. The documents may come

INTRODUCTION

Fault diagnosis methods include some actions implemented in order to detect, isolate and identify any abnormal phenomenon on a system. One can find in [START_REF] Isermann | Fault-Diagnosis System[END_REF]), [START_REF] Patton | Issues of Fault Diagnosis for Dynamic Systems[END_REF]) and the references therein the classical approaches using analytical information and which allow to do robust fault diagnosis in the presence of unknown entries and parametric uncertainties, namely parity space, observer and parameter estimation. These methods depend not only on structural knowledge of the system, but also require knowledge of system parameters that can be more or less accurate. The fault diagnosis approach proposed in this paper, contrary to the usual approaches, does not need a process identification procedure. In particular, the only assumptions that we will make concern:

• the process structure,

• the fault structure.

It means that:

-for the first assumption, we consider only known the degree of the differential equation modelling the system although its coefficients are assumed to be unknown and constant.

-for the second assumption, we consider structured additive actuator and sensor faults [START_REF] Fliess | An algebraic framework for linear identification[END_REF]) which take a large variety of signals into account.

For the design of the fault diagnosis algorithm, we adopt a distributional formulation (also used in [START_REF] Fliess | An algebraic framework for linear identification[END_REF]) and (Belkoura and Richard [2006])). This formulation allows us to transform the diagnosis problem into spectrum analysis problem from which we can detect and isolate faults without estimating system parameters. The paper is organized as follows. Section 2 is devoted to the recall of mathematical tools which will be used along the paper. In section 3, different assumptions are made on the system structure and the fault signal structures before defining and solving the fault diagnosis problem. Section 4 illustrates our approach with an academic example.

MATHEMATIC TOOLS

Distribution theory

We recall in this section some definitions and results from the distribution theory [START_REF] Schwartz | Théorie des distributions[END_REF]) and fix the notations we will use later.

Let K be an open set of R, the support of a function f defined on set K is the closure of the set of points where the function is not zero and is written Supp(f ). The space of smooth functions (infinitely differentiable) with compact support in K is denoted by D(K), D ′ (K) is the space of distributions on K, i.e. the space of linear and continue functionals on D(K). Distributions extend the concept of derivative to all locally integrable functions. If a function f is locally measurable on K, we define the regular distribution T f in D ′ (K) for all φ ∈ D(K) by < T f , φ >= f (s)φ(s)ds. In addition, if function f is continuous except at point x with a finite jump s x , the associated distribution derivative is given by Ṫf = ḟ -s x δ x , where ḟ is the usual derivative of function f . The complement of the largest open subset of K in which the distribution T vanishes 1 is called the support of T ; we denote it also by Supp(T ). A subspace of distributions of great interest in control theory is the subspace D ′ + of distributions with left bounded support (contained in [0, +∞[). Equipped with convolution and addition operations, this subspace becomes a convolution algebra with Dirac delta distribution δ as convolution identity. The Dirac distribution with delay τ is written as δ τ . The distribution T is said to be finite order if there exists a natural number m so that, for all compact subsets Ω ⊂ K with no empty interior

∃C > 0 : ∀φ ∈ D(Ω), |T (φ)| ≤ C sup 0≤i≤m ||φ (i) || ∞
The smallest number m satisfying the above inequality is called the order of the distribution T and is written as Ord(T ). Example 2.1. Supp(δ τ ) = {τ } and Ord(δ τ ) = 0. Supp(T f ) = Supp(f ) and Ord(T f ) = 0.

Derivation, integration and translation can be formed from the convolution product ẏ = δ (1) * y, t 0 y = H * y, y(t -τ ) = δ τ * y, where H is the unit step function (Heaviside distribution). The next theorem is the main result from which we will design our fault diagnosis algorithm. Theorem 2.1. [START_REF] Schwartz | Théorie des distributions[END_REF]) If a distribution T has a compact support Supp(T ) and a finite order m then αT = 0 whatever smooth function α and all its derivatives of order≤ m vanish on Supp(T ).

According to the above theorem, it follows

t k δ (n) =    0 if k > n (-1) k n! (n -k)! δ (n-k) τ else (1)
Finally, if the convolution product of distributions S and T is well defined, we have the following equality

t n (S * T ) = n k=0 C k n (t k S) * (t n-k T ), ∀t ∈ R, n ∈ N (2)

Spectral analysis

Given matrices Ã0 , Ã1 , .., Ãk ∈ R m×n , we assume that there exists at least one vector of reals β which reduces the rank of the pencils Ã0 -β1 Ã1 -β2 Ã2 -... -βk Ãk , i.e., there exists X ∈ R n /{0} so that ( Ã0 -β1 Ã1 -β2 Ã2 -... -βk Ãk ) X = 0 (3) Consider now matrices A 0 , A 1 , .., A k originated from matrices Ã0 , Ã1 , .., Ãk by perturbation (additive noises). The problem we want to solve in this section is the following: Given perturbed matrices A 0 , A 1 , .., A k (without knowing matrices Ã0 , Ã1 , .., Ãk ), we search vector β = (β 1 , ..., β k ) solution of the problem

min ||X||=1,λ ||(A 0 -λ 1 A 1 -λ 2 A 2 -... -λ k A k )X|| 2 ( Ã0 -λ 1 Ã1 -λ 2 Ã2 -... -λ k Ãk )X = 0 (4) Note that if || • || is the 2-norm, then min ||X||=1 ||(A 0 -λ 1 A 1 -... -λ k A k )X|| = σ min (A 0 -λ 1 A 1 ... -λ k A k )
and β can be characterized by

β = argmin λ σ min (A 0 -λ 1 A 1 ... -λ k A k )
where σ min (•) denotes the smallest singular value of •.

Let, for each ǫ ≥ 0, Λ ǫ (A 0 , ..., A k ) be the set

{λ ∈ R k : σ min (A 0 -λ 1 A 1 -.. -λ k A k ) ≤ ǫ} Sets Λ ǫ (A 0 , ..., A k )
are most likely to be empty for very small values of ǫ. As this value grows, we obtain a family of nonempty sets Λ ǫ (A 0 , ..., A k ). Then, β can be estimated as local minimum of this famly of sets. The case k = 1 which has been the subject of a wide literature ( [START_REF] Boutry | The generalized eigenvalue problem for non-square pencils using a minimal perturbation approach[END_REF]) and [START_REF] Wright | Pseudospectra of rectangular matrices[END_REF])) is called generalized eigenvalue problem. If k > 1, we will show that, with a suitable choice of dimensions m and n, one can always bring back to the case k = 1. Definition 2.1. (cokernel) The cokernel of matrix M , denoted by coker(M ), is a matrix N such that matrix N M is the zero. A matrix can admit more than one matrix as cokernel.

Let m, n ∈ N with m > n and A a matrix in R m×n . There exists an orthogonal matrix Q ∈ R m×m and an upper triangular matrix R ∈ R m×n so that A = QR. This decomposition is called qr decomposition. It also gives a cokernel of dimension ((m-n)×m) of matrix A, it's matrix N which rows are the transpose of the m -n last columns of matrix Q (i.e. N = Q(:, n + 1 : m) T ). Indeed Q being orthogonal, we have

N Q = (0 m-n,n I m-n ) and N A = N QR = [0 m-n,n I m-n ]R = 0 m-n,n R(1 : n, :) + I m-n R(n + 1 : m, :) = 0 m-n,n
We give now a solution of the above problem. Let m > (k -1)n and consider N k the cokernel of matrix A k with dimension (m -n) × m computed from qr decomposition. Then, we can write

σ min (A 0 -β 1 A 1 -.. -β k A k ) = σ min (N k A 0 -β 1 N k A 1 -... -β k-1 N k A k-1 ) because σ min is invariant under orthogonal transforma- tions. Now, let N k-1 ∈ R (m-2n)×(m-n) be the cokernel of matrix N k A k-1 . We have σ min (N k A 0 -β 1 N k A 1 -... -β k-1 N k A k-1 ) = σ min (N k-1 N k A 0 -... -β k-2 N k-1 N k A k-2 )
and so on, this process gives the following equality

σ min (A 0 -β 1 A 1 -... -β k A k ) = σ min (P 1,k A 0 -β 1 P 1,k A 1 )
where

P 1,k = N 2 N 3 ...N k and N j is a cokernel of di- mension (m -(k -j + 1)n) × (m -(k -j)n) of matrix N j+1 N j+2 ...N k A j , 2 ≤ j ≤ k -1 (that justifies the choice m > (k -1)n).
These operations allow to determinate β 1 . The estimations of β 2 ,..., β k can be achieved by applying the same steps.

The obvious algorithm to calculate argmin λ σ min (A -λB)

where A and B are matrices of same dimension is to evaluate σ min (A -λ i B) for values of λ i on a grid in the complex plane and then to generate a contour plotted from this data. Finally, we can take the possible minima on the intersection of the real axis and the contour. The number of possible local minima depends on the matrix rank. For example, in the above study, for each element β i , the number of distinct solutions varies from 1 to inf (n, m-(k -1)n). Additional constraints on β i permit to search particular solution of the problem.

In the next section, we show the interest of concepts presented in this section. The design of the fault diagnosis algorithm leads to solve a problem expressed as (4), where matrices A i are composed of process signal measurements and vector β contains the useful parameter to accomplish the fault diagnosis task.

FAULT DIAGNOSIS

Before detailing the algorithm stages, we make in the next subsections some assumptions about the system structure and the fault signals.

Model structure

Consider system ( ) where control signal u r (t) and output signal y r (t) satisfy a differential equation described by

a n y (n) r + a n-1 y (n-1) r + ... + a 0 y r = bu r y (i) r (0), i = 0, ..., n -1 (5)
Parameter a i , b and the initial conditions are supposed unknown and constant. In the distributional formulation (by kipping the same notation for the signals and their associated distributions), the system (5) becomes

a n y (n) r + a n-1 y (n-1) r + ... + a 0 y r = bu r + φ 0 (6
) where φ 0 is a linear combination of derivatives Dirac distribution of order less or equal to n -1 and contains the contributions of the initial conditions.

When an actuator fault denoted f a and a sensor fault denoted f c occur, then the system is controlled by u(t)+f a and generates output y(t)-f c where u(t) is the true control computed by the controller and y(t) is the measured output, as illustrated in figure (1). Then, signals u and y satisfy an equation of the form:

a n y (n) + ... + a 0 y = φ 0 + bu + bf a + a n f c (n) + ... + a 0 f c (7) u u r y r f a f c y Figure 1
. System with actuator and sensor faults.

Fault signal structure

We deal with fault signals modelled by structured signals [START_REF] Fliess | An algebraic framework for linear identification[END_REF]). A structured signal can be defined, in an informal way, as a solution of a linear differential equation. The main fault signals found in literature (abrupt, ramp, intermittent faults) can be modelled as structured signals [START_REF] Fliess | Algebraic changepoint detection[END_REF]). For this presentation, we deal with the case of abrupt faults f a and f c modelled as delayed step signals, i.e.

f a = l a H(t -τ a ) (8) f c = l c H(t -τ c ) (9)
where τ a , l a and τ c , l c are respectively the occurrence time and the magnitude of faults f a and f c . This assumption concerning the structure of faults f a and f c is not a restriction since the proposed algorithm can be applied for all types of faults modelled by structured signals.

System with actuator fault

Under the assumption of the occurrence of an abrupt actuator fault, equation ( 6) becomes

a n y (n) + .. + a 0 y = φ 0 + bu + bl a H(t -τ a ) (10)
From this model structure, a fault diagnosis algorithm is designed by assuming known only order n, control u(t) and output y(t). Thus a i , b and φ 0 are supposed to be unknown. This algorithm is composed of two stages, the generation of redundancy relations and the analysis of these relations.

In that follows, steps needed to generate redundancy equations are described.

First, let us derive the two members of equation ( 10) to obtain a n y (n+1) + .. + a 0 y (1) -bu (1) = φ

(1) 0 + bl a δ τa (11) where δ τa and φ

(1) 0 (singularities) have the following properties:

(1) Supp(φ To eliminate these singularities, let us just multiply equation (11) by a suitable function (theorem 2.1). Set

α 1 (t) = t n+1 (t -τ a )
(12) then α 1 (a n y (n+1) + .. + a 0 y (1) -bu (1) ) = 0 (13) because t n+1 φ

(1) 0 = 0 and (t -τ a )δ τa = 0. It follows, from (13) that for t ≥ 0

t n+2 [a n y (n+1) + ... + a 0 y (1) -bu] - τ a t n+1 [a n y (n+1) + ... + a 0 y (1) -bu] = 0 (14)
or, by setting for i = 1, 2 and j = 0, .., n W i,j = t n+i y (j+1) and

V i = t n+i u (1) n j=0 a j W 2,j -bV 2 -τ a n j=0 a j W 1,j -bV 1 = 0 (15)
W i,j and V i can be expressed, according to (1) and (2), like a linear finite combination of the derivatives of products t k y and t k u respectively. Indeed, we have

W i,j = inf (n+i,j+1) k=0 (-1) k C n+i k (j + 1)! (j + 1 -k)! δ (j+1-k) * (t n+i-k y) V i = δ (1) * (t n+i u) -(n + i)(t n+i-1 u)
Equality ( 15) is satisfied in spite of τ a (τ a is not identifiable before fault occurence) because, for t < τ a t n [ n i=0 a i y (i+1) -bu (1) ] = 0 However, for t ≥ τ a , the redundancy relations obtained by applying successive convolutions on equation ( 15) by H lead to the following spectral formulation

[A 2 -τ a A 1 ]X = 0 (16)
by setting, for i = 1, 2

A i =        H * n+1 * W i,n ... H * n+1 * V i H * n+2 * W i,n ... H * n+2 * V i . . . . . . . . . H * n+m * W i,n ... H * n+m * V i        X = ( a n a n-1 . . . a n-p . . . a 0 -b ) T (17) m > 0 (number of rows), H * p = H * H * ... * H p times .
The initial successive convolution (H * n+1 ) ensures the elimination of all derivatives (numerically less robust than integrations).

In general, for any structured fault, the generation of redundancy relations is the same than the previous one and consists of three steps (1) Derivation (apparition of singularities) (2) Multiplication (elimination of singularities) (3) Integration (spectral formulation) From expression ( 16), the estimation of τ a consists in the computation of generalized eigenvalues of A 2 and A 1 . In practice, signals u and y being corrupted by noises, we rather have τ a ∈ Λ ǫ (A 2 , A 1 ), ǫ > 0 (see section 2.2). Matrices A 2 and A 1 being time-varying, thus the estimation of τ a is also. At each time t, the value of τ a estimated (chosen among the possible values) consists of the generalized eigenvalue associated with vector X wich is the closest possible to the vector associated with the earlier value.

The decision rule is based on the time behaviour of τ a and l a estimations (characteristics of the fault).

For an actuator fault (here delayed step), the estimation of τ a becomes stationary just after the time of the fault occurrence.

After having estimated τ a , a similar reasoning leads to the estimation of l a . Indeed, setting α 2 = t n+1 yields α 2 (a n y (n+1) + .. + a 0 y (1) -bu (1) ) = bl a α 2 δ τa ⇒ ∀ t ≥ τ a and k ≥ n + 1:

l a = H * k * [α 2 (q n+2 y (n+1) + .. + q 2 y (1) + q 1 u (1) )] H * k * [-q 1 α 2 δ τa ]
where q = (q n+2 , ..., q 1 ) T is the "smallest singular vector" 2 of matrix A 2 -τ a A 1 .

2 associated to the smallest singular value

The successive convolutions H * p (successive integrations) proceed like a low-pass filter and contribute to attenuate high frequency noises.

For the generation of matrices A i , we have taken the polynomial function α 1 in order to simplify computations. Damped functions, like α = e -at t n+1 (t -τ a ), a > 0, can improve the approach robustness to measurement noises. Indeed, these functions contribute to attenuate noises, in opposition to polynomial functions which accentuate the noise effects.

System with sensor faults

In this part, it is provided the algorithm ensuring system fault diagnosis when ( ) is corrupted by a sensor fault. This fault is modelled as a delayed step. We keep the same assumptions as previously.

So, under the assumption of the occurrence of an abrupt sensor fault at time τ c with a magnitude l c , the faulty system equation becomes

a n y (n) + ... + a 0 y = φ 0 + bu +l c (a n δ (n-1) τc + ... + a 0 H(t -τ c )) (18)
The derivation yields to

a n y (n+1) + ... + a 0 y (1) -bu (1) = φ (1) 0 + l c (a 1 δ (n) τc + ... + a 0 δ τc ) (19)
To eliminate the singularities in the second member of equation ( 19), we can consider the function

α 3 (t) = t n+1 (t -τ c ) n+1 (
20) which yields (according to theorem (2.1)) α 3 (a n y (n+1) + .. + a 0 y (1) -bu (1) ) = 0 (21)

By setting for i = 1, .., n + 2 and j = 0, .., n W i,j = t n+i y (j+1) and V i = t n+i u (1) we obtain the spectral formulation

[A n+2 -λ n+1 A n+1 -... -λ 2 A 2 -λ 1 A 1 ]X = 0 (22)
where A i , for i = 1, ..., n + 1 and X are defined like in (17) and

λ i = (-1) n-i+1 C n+1 i-1 τ n+2-i c , 1 ≤ i ≤ n + 1 is the coefficient associated to term t i-1 in the polynomial function -(t -τ c ) n+1 .
From ( 22), τ c can be estimated with the estimation of any λ i . Since only λ n+1 gives a linear estimation of τ c (λ n+1 = (n+1)τ c ), then the estimation is obtained through the computation of Λ ǫ (P A n+2 , (n + 1)P A n+1 ) where P = coker(coker(...coker(A 1 )A 2 ...)A n-1 )A n ) The decision rule is made like in the previous subsection.

As well, the amplitude l c is estimated after the estimation of τ c by setting α 4 = t n+1 (t -τ c ) n : ∀ t ≥ τ c and k ≥ n + 1 l c = (-1) n H * k * [α 4 (q n+2 y (n+1) + .. + q 2 y (1) + q 1 u (1) )]

n!H * k * δ τc where q = (q n+2 , ..., q 1 ) T is the smallest singular vector of matrix P A n+2 -(n + 1)τ c P A n+1 .

System with both actuator and sensor faults

In the case of multiple faults (actuator and sensor faults), the use of the previous approaches do not ensure the system fault diagnosis.

Assuming for example the presence of an actuator fault f a (t) = l a H(t -τ a ) and a sensor fault f c (t) = l c H(t -τ c ), the corresponding model is

a n y (n) + .. + a 0 y = φ 0 + bu (23) +bl a H(t -τ a ) + l c (a n δ (n-1) τc + .. + a 0 H(t -τ c ))
The multiplication of the previous equation by

α 1 = t n+1 (t -τ a ) or α 3 (t) = t n+1 (t -τ c ) n+1
does not eliminate all singularities obtained after derivation. For this, we can consider the function

α 5 = t n+1 (t -τ a )(t -τ c ) n+1 The spectral expression obtained thanks to α 5 is [A n+2 -β n+1 A n+1 -... -β 1 A 1 -β 1 A 1 ]X = 0 (
24) where A i and X are obtained according to (17) and β i is the coefficient associated to term t i-1 in the polynomial function -(t -τ a )(t -τ c ) n+1 . From ( 24), τ a and τ c can be estimated, as for sensor fault case, by the means of the estimation of any β i . The first fault can be detected like in subsections 3.3 and 3.4. The detection of the second fault can be made through the estimation of β n+1 = τ a + (n + 1)τ c which becomes stationary after time t ≥ max(τ a , τ c ).

Finally, the second fault magnitude can be also estimated through τ a and τ c estimations.

APPLICATION

Consider the first order system (n = 1) modelled by the input-output equation a 1 ẏ + a 0 y = bu This system is controlled with a PI-controller. The proposed approach is applied to diagnose the system above, simulated in first with a sensor drift and then with actuator and sensor bias. Recall that the proposed appoach assumes that parameters a 1 , a 0 and b are unknown and will not be explicitely estimated. 4) show the graphical results of the estimations (of τ c and l c ) given by the proposed algorithms. According to these observations, we can conclude to the occurence of a sensor drift, modelled as f c (t) = 0.3(t -1)H(t -1). According to these two observations, the detection of an actuator fault, modelled by f a (t) = 0.8H(t -1), can be accomplished in the time interval [1,2]. From t > 2, the estimation of τ a and l a fails due to the occurrence of a sensor fault. It is important to notice that the later cannot The algebraic approach for fault diagnosis presented in this paper is only based on the structural knowledge of the system and the faults. On the one hand, the parameters of the system can be unknown and we do not need to estimate them explicitly. On the other hand, the decision rule is not based on the outputs of some residual generators, but rather on the time behaviour of some fault characteristic estimations. Simulation results show that the proposed method gives good results for detecting and isolating actuator and sensor faults even in the case of multiple faults. In future works, robustness of the proposed method with respect to noise signals will be investigated.

  Supp(δ τa ) = {τ a } and Ord(δ τa ) = 0

Figure ( 2 Figure 2 .

 22 Figure (2) shows input u and output y simulated with a drift (of slope l c = 0.3 at instant τ c = 1) of the output sensor.

Figure 4 .

 4 Figure 3. evolution of τ c estimation

Figure 5 .

 5 Figure 5. Input and output signals in the case of multiple faults Graphical results of the estimations given by the proposed algorithms are represented on figures (6)-(10). On figure (6), we can easily see the stationary value of τ a estimation around value τa = 1 in the time interval [1, 2]. The associated magnitude estimation is equal to 0 before time t = 1 and takes a value around 0.8 in the time interval [1, 2] (figure 7).

Figure 7

 7 Figure 6. time evolution of τ a estimation

Figure 8 .

 8 Figure 8. time evolution of τ c estimation

  Figure (10) shows the time evolution of the estimation of the sensor fault magnitude obtained by means of β 2 and τ a estimations. Thanks to τ a and β 2 estimations, the detection of a sensor bias can be accomplished, it is about f c (t) = 0.3H(t -2).

  Figure 9. time evolution of β 2 estimation

We said that a distribution T vanishes on Ω ⊂ K if T (φ) = 0 for all φ ∈ D(Ω)


Mirjam Dür

Jean-Baptiste Hiriart-Urruty

Keywords: Copositive matrices, difference-of-convex (d.c.), Legendre-Fenchel transforms, nonconvex duality, subgradient algorithms for d.c. functions

We consider the problem of minimizing an indefinite quadratic form over the nonnegative orthant, or equivalently, the problem of deciding whether a symmetric matrix is copositive. We formulate the problem as a difference of convex functions (d.c.) problem. Using conjugate duality, we show that there is a one-to-one correspondence between their respective critical points and minima. We then apply a subgradient algorithm to approximate those critical points and obtain an efficient heuristic to verify non-copositivity of a matrix.

Testing copositivity with the help of difference-of-convex optimization 1 Introduction

Consider the problem of minimizing a (homogeneous) quadratic form over the nonnegative orthant:

(P) min 1 2 x T Ax s. t. x ∈ R n + . As we do not assume the matrix A to be positive semidefinite, (P) is a nonconvex constrained minimization problem. It is equivalent to the problem of determining whether A is copositive: A matrix A is called copositive, if x T Ax ≥ 0 holds for all x ∈ R n + . Copositivity has recently attracted quite an amount of interest in mathematical optimization, see [START_REF] Bomze | Copositive optimization -recent developments and applications[END_REF][START_REF] Dür | Copositive Programming -a Survey[END_REF][START_REF] Hiriart-Urruty | A variational approach to copositive matrices[END_REF] for surveys. The problem of determining copositivity, and consequently also problem (P), are known to be co-NP-complete [START_REF] Murty | Some NP-complete problems in quadratic and nonlinear programming[END_REF].

In this note, we transfer the copositivity problem into an unconstrained problem in the same number of variables. A similar approach has been pursued by Bomze and Palagi [START_REF] Bomze | Quartic Formulation of Standard Quadratic Optimization Problems[END_REF], who use an exact penalization method to arrive at an objective function which is a polynomial of degree six. Our approach differs from theirs in that we formulate (P) as an unconstrained difference-of-convex (d.c. in short) optimization problem which allows us to use specific duality schemes from convex analysis and nonconvex optimization. As opposed to the objective function of Bomze and Palagi [START_REF] Bomze | Quartic Formulation of Standard Quadratic Optimization Problems[END_REF] which is C ∞ , we derive a nonsmooth objective which is a combination of a quadratic function and the indicator function of the convex feasible set.

To be more specific, let r > max{λ max (A), 0}, where λ max denotes the maximal eigenvalue of A. Then the matrix rI -A is positive definite. We now decompose the objective function as follows: For all x ∈ R n , let

f (x) = g(x) -h(x), with g(x) = r 2 x 2 + Ψ R n + (x), h(x) = 1 2 x T (rI -A)x, where Ψ R n + (x) is the indicator function of R n + , that is: Ψ R n + (x) = 0 if x ∈ R n + and Ψ R n + (x) = +∞ if x / ∈ R n + .
Note that both g and h are convex functions, whence g -h is a description of f as a d.c. function.

On R n + , the feasible set of (P), f coincides with the objective function of (P). Outside R n + , f ≡ +∞. Hence, (P) translates to the unconstrained d.c. minimization problem

(P) min x∈R n f (x) = g(x) -h(x).
The copositivity property of A is now that inf x∈R n f (x) = 0, (actually, the minimal value 0 of f on R n is attained since f (0) = 0). In dealing with d.c. minimization problems, there is a duality scheme specially tailored for these problems which has proved useful in various contexts, starting with calculus of mechanics, cf. [START_REF] Toland | A duality principle for non-convex optimization and the calculus of variations[END_REF][START_REF] Toland | On the stability of rotating heavy chains[END_REF][START_REF] Toland | Duality in non-convex optimization[END_REF]: the adjoint or dual problem associated with (P) is

(P ) min f (y) := h * (y) -g * (y) s. t. y ∈ R n ,
where ϕ * stands for the Legendre-Fenchel conjugate of ϕ (see [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms I and II[END_REF]Chapter X]).

In our situation, for all y ∈ R n we have

h * (y) = 1 2 y T (rI -A) -1 y, g * (y) = 1 2r y + 2 ,
where y + means the vector (y + 1 , . . . , y + n ) with y + i = max{y i , 0}. The new problem (P ) is again a d.c. minimization problem, now unconstrained, with a C 1 (but not C 2 ) objective function. Note that with the help of the Legendre-Fenchel transformation, the nonnegativity constraint in (P) has been integrated in the expression y + 2 appearing in (P ).

The construction of (P ) from (P) depends the decomposition f = g -h more than on f itself. This flexibility is indeed an advantage, for example adding a convex function ϕ to both g and h would not change the initial problem since f = (g + ϕ) -(h + ϕ) = g -h, but would give rise to a different adjoint problem (P ). Most of the arguments in this paper hold for d.c. decompositions other than the one described above, but the numerical performance of d.c. optimization algorithms will depend on the choice of the decomposition. Bomze and Locatelli [START_REF] Bomze | Undominated d.c.decompositions of quadratic functions and applications to branch-and-bound approaches[END_REF] discuss this problem. The authors give evidence that in case of a quadratic objective x T Ax, the so called spectral d.c. decomposition is preferable. This decomposition consists of writing A as A = A + -A -, where both A + and A -are (close to) nonsingular and collect the positive (negative) eigenvalues of A, respectively. In our approach, however, we prefer to work with nonsingular decompositions, as we will frequently use the one-to-one relation between problems (P) and (P ).

The two problems (P) and (P ) are related in an "involutive" manner, i.e., (P ) := (P ) = (P), since ϕ * * := (ϕ * ) * = ϕ holds for any closed convex function ϕ.

The next results relate infima, global minimizers and critical points in problems (P) and (P ). They are adaptations to our specific context of general results in d.c. dualization ( [START_REF] Toland | A duality principle for non-convex optimization and the calculus of variations[END_REF][START_REF] Toland | Duality in non-convex optimization[END_REF], see also [START_REF] Singer | A Fenchel-Rockafellar type duality theorem for maximization[END_REF]).

Comparing infima and critical points in (P) and (P )

In minimizing f over R n + (our problem (P)), there are only two possibilities: either f is not bounded below, i.e., inf x∈R n + f (x) = -∞, or f is bounded from below. In that case, necessarily inf x∈R n + f (x) = 0 and x = 0 is among the global minimizers of f on R n + . Like for quadratic objective functions, this peculiar property is due to the positive homogeneity of degree 2 of f : if there were x ∈ R n + with f (x) < 0, then f (αx) = α 2 f (x) for all α > 0 would imply inf

x∈R n + f (x) = -∞. Theorem 1 We have inf x∈R n + x T Ax = inf x∈R n f (x) = inf y∈R n f (y),
that is, the infima of (P) and (P ) (which are either 0 or -∞) are equal.

This relationship has been proved directly in [START_REF] Hiriart-Urruty | A variational approach to copositive matrices[END_REF]Section 7.8] where it gave rise to the following characterization of copositivity: A matrix A is copositive if and only if y T (rI -A) -1 y ≥ 1 r y + 2 for all y ∈ R n . We next define what we mean by a critical point of f on R n + , or critical point of (P).

Definition 1 A point x ∈ R n + is called a critical point of (P) if -Ax ∈ N R n + (x), (1)
where N R n + (x) stands for the normal cone to R n + at x.

As a motivation for this definition, note hat ∇h(x) = (rI -A)x and that the generalized subdifferential of g at x is ∂g(x) = rx + N R n + (x), whence writing the optimality condition 0 ∈ ∂f (x) = ∂g(x) -∇h(x) directly leads to [START_REF] Bomze | Copositive optimization -recent developments and applications[END_REF].

The set

R n + is a closed convex cone whose polar cone is R n -:= -R n + . Hence, the normal cone N R n + (x) can be given explicitly (see [7, Chapter III]): ȳ = (ȳ 1 , . . . , ȳm ) ∈ N R n + (x) ⇔ ȳ ∈ R n -and ȳi = 0 whenever xi > 0.
As a consequence, we have the following.

Proposition 1 A point x ∈ R n + is a critical point of (P) if and only if Ax ≥ 0 and (Ax) i = 0 for all i such that xi > 0. (2)
Observe that x = 0 is always a critical point of (P), and so is any local minimizer of f on R n + . Due to the positive homogeneity of degree 2 of the objective function f in (P) (i.e., f (αx) = α 2 f (x) for all α > 0 and x ∈ R n ), positive homogeneity also holds true for the set of critical points in (P): if x is a critical point of (P) and α > 0, then αx is also a critical point of (P). Therefore, we have "critical rays".

In [START_REF] Seeger | Local minima of quadratic forms on convex cones[END_REF], the authors considered the problem of minimizing the function

1 2 x T Ax over S = R n + ∩ {x ∈ R n | x = 1}.
The resulting definition of critical point they adopted ([15, Theorem 1]),

Ax -(x T Ax)x ∈ R n + is less stringent than ours.
In problem (P ), the objective function f is differentiable, so that defining critical points does not pose any difficulty.

Definition 2 A point ȳ ∈ R n is a critical point of (P ) if (rI -A)ȳ + = rȳ. (3) Indeed, ∇f (ȳ) = ∇h * (ȳ) -∇g * (ȳ) = (rI -A) -1 ȳ -1 r
ȳ+ , so that the relation ∇f (x) = 0 directly leads to [START_REF] Bomze | Undominated d.c.decompositions of quadratic functions and applications to branch-and-bound approaches[END_REF].

Observe that ȳ = 0 as well as any local minimizer of f on R n are always critical points. Like for (P), the positive homogeneity of degree 2 of the objective function f in (P ) yields the positive homogeneity of the set of critical points of (P ): if ȳ is a critical point of (P ), then αȳ is a critical point of (P ) for any α > 0.

There are some specific properties of the critical points in (P ) which can be derived from [START_REF] Bomze | Undominated d.c.decompositions of quadratic functions and applications to branch-and-bound approaches[END_REF]. Here are a few examples: Proposition 2 The critical points of (P ) have the following properties:

1. The critical points located in R n + are the points in R n + ∩ Ker A. Hence, if A is nonsingular, there is no critical point in R n + , except for ȳ = 0. 2.
There is no critical point in R n -, except for ȳ = 0. 3. All the critical points belong to the cone

K := R n + -A(R n + ). ( 4
) Proof 1. The characterization (3) of critical points ȳ located in R n + is (rI -A)ȳ = rȳ,
that is to say Aȳ = 0.

If ȳ lies in R n

-, then ȳ+ = 0, so (3) reads rȳ = 0 whence the result follows. 3. According to (3), a critical point ȳ equals ȳ+ -1 r Aȳ + which belongs to the cone

R n + -A(R n + ). 2
With the help of ȳ-= (y - 1 , . . . , y - n ) where y - i = max{-y i , 0} so that y = y + -y -, we can rewrite characterization (3) as follows: y is a critical point in (P ) if and only if

Ay + = ry -.

While the notion of a critical point in (P) does not depend on r (still chosen greater than max{λ max (A), 0}), the critical points in (P ) depend a priori on r. However, the closed convex cone K defined in (4) does not depend on r. It gives a rough estimate of the location of the critical points in (P ) for all the admissible r. The main point in this section is the correspondence between critical points (as well as critical values) in the original problem (P) and in the adjoint problem (P ).

Theorem 2 1. If x is a critical point in (P), then ȳ := (rI -A)x is a critical point in (P ), and the corresponding critical values are equal:

f (x) = f (ȳ). 2. If ȳ is a critical point in (P ), then x := (rI -A) -1 ȳ = 1
Since (g * ) * = g and (h * ) * = h, we deduce from ( 6):

(g * ) * (x) + g * (ȳ) = xT ȳ (h * ) * (x) + h * (ȳ) = xT ȳ (7)
as well as

g(x) -h(x) = h * (ȳ) -g * (ȳ). (8)
Equations ( 7) state that x ∈ ∂g * (ȳ) ∩ ∂h * (ȳ). Since both g * and h * are differentiable at ȳ, this means that ∇g * (ȳ) = ∇h * (ȳ) or equivalently, 0 = ∇h * (ȳ)-∇g * (ȳ) = ∇f (ȳ). So ȳ is indeed a critical point of (P ). Equality of the critical values has been observed in (8). 2. We play the same game as above. Let ȳ be a critical point of (P ). We have:

(rI -A) -1 ȳ = ∇h * (ȳ) = ∇g(ȳ) = 1 r ȳ+ . (9)
Define x := (rI -A) -1 ȳ. Similar to (6), we have

x = ∇h * (ȳ) ⇔ (h * ) * (x) + h * (ȳ) = xT ȳ, x = ∇g * (ȳ) ⇔ (g * ) * (x) + g * (ȳ) = xT ȳ. (10)
Again, since (h * ) * = h and (g * ) * = g, we derive from [START_REF] Hall | Copositive and Completely Positive Quadratic Forms[END_REF]:

h(x) + h * (ȳ) = xT ȳ g(x) + g * (ȳ) = xT ȳ (11)
and

h * (ȳ) -g * (ȳ) = g(x) -h(x). (12)
What [START_REF] Mongeau | Computing eigenelements of real symmetric matrices via optimization[END_REF] expresses is:

ȳ = ∇h(x) and ȳ ∈ ∂g(x),
that is to say: ∇h(x) ∈ ∂g(x), a formulation equivalently saying that x is a critical point in (P). The fact that x = (rI -A) -1 ȳ ∈ R n + was a by-product of [START_REF] Thi | Solving a class of linearly constrained indefinite quadratic problems by D.C. algorithms[END_REF]. Equality of the critical values, i.e., f (ȳ) = f (x) was observed in [START_REF] Murty | Some NP-complete problems in quadratic and nonlinear programming[END_REF].

In the scheme below we summarize the correspondence between critical points in (P) and (P ) via the nonsingular linear mapping rI -A:

(P) (P ) critical point x rI-A -→ critical point ȳ = (rI -A)x critical point x = (rI -A) -1 ȳ (rI-A) -1
←-critical point ȳ

3 Comparing global and local minimizers in (P) and (P )

We have shown in Section 2 that the infima in (P) and (P ) are the same. Now, using the same techniques from convex analysis as in Section 2, we show that there is also a one-to-one correspondence between global minimizers in (P) and (P ).

Theorem 3 Let x be a global minimizer in (P).

Then (i) xT Ax = 0, (ii) ȳ := (rI -A)x is a global minimizer in (P ), (iii) ȳT (rI -A) -1 ȳ = 1 r ȳ+ 2 .
Conversely, let ȳ be a global minimizer in (P ).

Then (i') ȳT (rI -A) -1 ȳ = 1 r ȳ+ 2 (ii') x := (rI -A) -1 ȳ = 1 r
ȳ+ is a global minimizer in (P).

Proof Let x be a global minimizer in (P). Then necessarily xT Ax = 2f (x) = 0. With ȳ := (rI -A)x, we have f (ȳ) = 0 (cf. Theorem 1), and it remains to prove that f (ȳ) = h * (ȳ) -g * (ȳ) = 0. By linearity of the constraints, no qualifications are needed to ensure that x is a critical point. Then Theorem 2 gives the desired results (ii) and (iii). The proofs of (i') and (ii') are in the same vein.

2

What we have proved for global minimizers also holds true for local minimizers: the proofs have to be slightly adapted. We present below the version "from (P ) to (P)", which is the case of more interest.

Theorem 4 If ȳ is a local minimizer in (P ), then x = (rI -A) -1 ȳ = 1 r ȳ+ is a local minimizer in (P).
Proof Let ȳ be a local minimizer in (P ). Then there exists a neighborhood V of ȳ such that

h * (y) -g * (y) = f (y) ≥ f (ȳ) = h * (ȳ) -g * (ȳ) for all y ∈ V. ( 13
)
Since ∇h : x ∈ R n → ∇h(x) = (rI -A)x is a bijective linear mapping, there exists a neighborhood N of x = (rI -A) -1 ȳ such that ∇h(N ) ⊂ V . So choose x ∈ N . Then y = ∇h(x) and x satisfy h * (y) + h(x) = x T y (this is the characterization of ∇h(x) in terms of h * ), and

g * (y) + g(x) ≥ x T y
by the definition of the conjugate fuction. Because y = ∇h(x) lies in V , these two inequalities combined with (13) imply

g(x) -h(x) ≥ h * (y) -g * (y) ≥ f (ȳ). ( 14
)
We know from Theorem 2 that f (ȳ) = f (x). So, we have proved with ( 14) that f (x) ≥ f (x) for all x in the neighborhood N of x, which completes the proof. 2

DC Algorithm

Any local optimization procedure can be used to compute critical points of (P) and (P ). In our context, it seems suitable to adapt the DCA algorithm of [START_REF] Pham | Algorithms for solving a class of nonconvex optimization problems. Methods of subgradients[END_REF][START_REF] Pham | Duality in D.C. (difference of convex functions) optimization Subgradient methods[END_REF][START_REF] Thi | Solving a class of linearly constrained indefinite quadratic problems by D.C. algorithms[END_REF]. This method works as follows: given a starting point x 0 ∈ R n , we construct sequences {x k } and {y k } with

y k ∈ ∂h(x k ) and x k+1 ∈ ∂g * (y k )
such that the sequences (g -h)(x k ) and (h * -g * )(y k ) are decreasing and such that the limit points of {x k } (resp. {y k }) is a critical point of (g -h) (resp. (h * -g * )).

The condition y k ∈ ∂h(x k ) is equivalent to x k ∈ ∂h * (y k ), which in turn means by definition of the subdifferential that

h * (y) -(x k ) T y ≥ h * (y k ) -(x k ) T y k for all y ∈ dom h * ,
or equivalently, that y k can be derived as

y k ∈ Argmin{h * (y) -(x k ) T y : y ∈ dom h * }. ( 15
)
The latter is equivalent to saying that y k is a solution of the problem

(D k ) inf y∈dom h * {h * (y) -g * (y k-1 ) -(x k ) T (y -y k-1 )}.
which can be interpreted as a linearized version of (P ), where we have linearized g * at the point y k-1 . Observe that in doing so, we have obtained a convex problem (D k ). Since x k ∈ ∂g * (y k-1 ) by construction, we moreover have that inf (P ) ≤ inf (D k ).

Analogously, the condition x k+1 ∈ ∂g * (y k ) says that x k+1 can be obtained as

x k+1 ∈ Argmin{g(x) -(y k ) T x : x ∈ R n }. (16)
which in turn is equivalent to saying that x k+1 is a solution to the problem

(P k ) inf x∈R n {g(x) -h(x k ) -(y k ) T (x -x k )}
. which just as above can be interpreted as a linearized version of (P), where we have linearized h at the point x k . Since y k ∈ ∂h(x k ), we get inf (P) ≤ inf (P k ).

Note that to actually compute the sequences {x k } and {y k }, we can use ( 15) and ( 16). In our particular setting, (15) reads as:

y k ∈ Argmin{ 1 2 y T (rI -A) -1 y -(x k ) T y : y ∈ R n }.
Since this is a differentiable problem, the solution is obtained by setting the gradient to zero, which gives

y k = (rI -A)x k .
Similarly, in our context [START_REF] Singer | A Fenchel-Rockafellar type duality theorem for maximization[END_REF] translates to

x k+1 ∈ Argmin{ r 2 x 2 -(y k ) T x : x ∈ R n + }.
Here, the objective function is separable, so we can consider each component function separately. For each of them, this amounts to minimizing a strictly convex quadratic function subject to a sign constraint, and the solution is easily seen to be

x k+1 i = 1 r y k i if y k i ≥ 0, 0 else, i.e., x k+1 = 1 r (y k ) + .
Combined, this yields the recursions:

x 0 given,

x k+1 = 1 r (rI -A)x k + = (I -1 r A)x k + y 0 = (rI -A)x 0 , y k+1 = (rI -A) 1 r (y k ) + = (I -1 r A) y k + .
A detailed convergence theorem for this method is proved in [START_REF] Pham | Duality in D.C. (difference of convex functions) optimization Subgradient methods[END_REF]. We present here a reduced version of this theorem which, however, is sufficient for our purposes:

Theorem 5 Suppose that the sequences {x k } and {y k } are defined as described above. Then we have

(h * -g * )(y k+1 ) ≤ (g -h)(x k+1 ) ≤ (h * -g * )(y k ) ≤ (g -h)(x k ). If (g -h)(x k+1 ) = (g -h)(x k ), then x k , x k+1 are critical points of g -h.
Similarly, if (h * -g * )(y k+1 ) = (h * -g * )(y k ), then y k , y k+1 are critical points of h * -g * .

A heuristic for testing non-copositivity

The DC algorithm described in the last section inspires the following easy heuristic for testing whether a matrix is not copositive.

Algorithm 1: Heuristic to detect whether a matrix A is not copositive.

Input: symmetric matrix A; parameter r > max{λmax(A), 0} Output: a certificate that A is not copositive or statement "this instance is undecidable by the heuristic" 1 while stopping criterion is not fulfilled do 2 generate a starting point x 0 ;

3 k ← 0; 4 repeat 5 x k+1 = (I -1 r A)x k + ; 6 if (x k+1 ) T A(x k+1 ) < 0 then 7 return "A is not copositive" 8 end 9 k ← k + 1;
10 until x k+1 = x k or stopping criterion is fulfilled; 11 end 12 return "this instance is undecidable by the heuristic" As a stopping criterion, one may choose for example that a prescribed number of iterations has been performed. Note that this algorithm has the ability to do many multistarts, and can be parallelized without any further effort.

Examples

We first study the behavior of our heuristic for a copositive input matrix. In particular, it is interesting to see the role of the parameter r (recall that by assumption r > max{λ max (A), 0}). Consider the so-called Horn-matrix

H =       1 -1 1 1 -1 -1 1 -1 1 1 1 -1 1 -1 1 1 1 -1 1 -1 -1 1 1 -1 1      
which is known to be copositive (see e.g. [START_REF] Hall | Copositive and Completely Positive Quadratic Forms[END_REF]). It has λ max (H) = 3.2361, so we need to choose r bigger than this value. Applying our heuristic to H will obviously never yield a certificate for non-copositivity.

Running the heuristic for different values of r and 1000 starting points each gives the results displayed in Table 1 iterations gives the respective number of iterations needed by the heuristic until a critical point was reached. In column "% failure" we list the percentage of instances where the heuristic did not yield a critical point within 1000 iterations.

Observe that convergence seems to be faster for values of r that are closer to λ max (H).

As a second example, consider the behavior of the heuristic when applied to noncopositive matrices. To construct such matrices, recall that the clique number ω of a graph can be written as:

ω = min{α ∈ R : [α(E -A) -E] is copositive },
where E is the all-ones matrix and A is the adjacency matrix of the graph. Now consider the 5-cycle on a graph with five vertices whose adjacency matrix is given by

C =       0 1 0 0 1 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1 1 0 0 1 0       .
Its clique number is ω = 2. Therefore, matrices of the form

B α := α(E -C) -E
are noncopositive for α < 2 and copositive for α ≥ 2. Note that α = 2 is the smallest value for which this is true, and that we have

B 2 = 2(E -C) -E = E -2C = H, the Horn-matrix from above.
Below is a summary of the behavior of the heuristic. We choose different values of α, and for the value of r we choose r = cλ max (B α ) with c ∈ {1.1, 2, 10}. Observe that as α → 2 from below, the matrix B α approaches the boundary of the copositive cone, since for α = 2 = ω, the matrix B ω lies on the boundary of that cone. In Table 2, we list the results for various values of α and r for 1000 starting points each. The column "noncopositive" gives the percentage of runs that noncopositivity was detected. The next columns list the percentage of runs where the method converged to a critical point. We also list in which iteration convergence occurred. Finally, the right-most column indicates the percentage of runs where the method failed (i.e., neither detected noncopositivity nor converged to a critical point).

It can be seen that as α → 2, the heuristic performs poorer. As discussed above, this is due to the fact that in this case B α approaches the boundary of the copositive cone. This behavior corresponds to what has been observed in other methods, see e.g. [START_REF] Bundfuss | Algorithmic Copositivity Detection by Simplicial Partition[END_REF].

Observe that in this setting, the choice of the parameter r does not seem to have any impact.

In all cases, if the heuristic detected noncopositivity, it was in the very first iteration, which means the sampled starting point has already been a certificate for noncopositivity.

r 2 xT

 2 ȳ+ is a critical point in (P), with equality in the critical values:f (ȳ) = f (x).Proof 1. Let x be a critical point in (P). By the very first definition, we have (rI -A)x = ∇h(x) ∈ ∂g(x) (5) and f (x) = 1 Ax. Let ȳ := (rI -A)x. By an elementary property of the Legendre-Fenchel conjugate, we have ȳ ∈ ∂g(x) ⇔ g * (ȳ) + g(x) = xT ȳ (or, equivalently, ≤ xT ȳ), ȳ = ∇h(x) ⇔ h * (ȳ) + h(x) = xT ȳ.

Table 1

 1 . In this table, the columns min/max/average Results for the Horn-matrix: Running the heuristic for different values of r and 1000 starting points each.

		iterations
	r	average min max % failure
	20	219.91	2	362	0.12
	10	102.99	2	171	0.28
	6	54.50	2	87	1.20
	4	28.74	2	55	2.93
	3.5	21.77	2	45	0.46
	3.2365	18.09	2	42	0.53

Table 2

 2 Results for the matrixBα := α(E -C) -E.

			noncopositive		critical point		failure
	α	c	%	%	av.	min	max	%
	0.5 1.1	66.3	33.7	2	2	2	0
		2	67.2	32.8	2	2	2	0
		10	66.2	33.8	2	2	2	0
	1	1.1	52.7	47.3	2.782	2	3	0
		2	52.0	48.0	2.79	2	3	0
		10	52.5	47.5	2.78	2	3	0
	1.5 1.1	9	91	306.45	2	5181	0
		2	8.4	91.6 355.23	2	5181	0
		10	10	90	280.99	3	5179	0
	1.7 1.1	6.3	0	-	-	-	93.7
		2	5.3	0	-	-	-	94.7
		10	5.5	0	-	-	-	94.5
	1.9 1.1	2.5	1	2	2	2	97.4
		2	2.4	1	2	2	2	97.5
		10	2.9	0	-	-	-	97.1

Acknowledgements

This research was initiated during a research visit of Mirjam Dür to Université Paul Sabatier in Toulouse. She would like to thank the members of the optimization group for their warm hospitality and inspiring discussions. Her work was partially supported by the Netherlands Organisation for Scientific Research (NWO) through Vici grant no.639.033.907. Both authors wish to thank the referees for highly detailed and valuable comments.