Charles Consel
email: consel@labri.fr

Generative Programming from a DSL Viewpoint

Introduction

A domain-specic language (DSL) is typically created to model a program family [START_REF] Consel | From A Program Family To A Domain-Specic Language[END_REF]. The commonalities found in the target program family suggest abstractions and notations that are domain specic. In contrast with general-purpose languages (GPL), a DSL is readable for domain experts, often concise, and usually declarative. As an illustration, consider a program family aimed to communicate data in a distributed heterogeneous system. Such a layer is commonly needed in a variety of distributed applications and relies on a mechanism like the Sun Remote Procedure Call (RPC) [START_REF]NFS: Network le system protocol specication[END_REF]. The XDR layer of the Sun RPC consists of marshaling and un-marshaling both arguments and returned value to/from a machine-independent format. The program family, represented by all the possible (un-)marshaling variations, has lead to the development of a DSL that allows a programmer to concisely express the type of the remote procedure arguments and returned value, and to obtain the corresponding marshaling layer on both the client and server sides.

From a DSL viewpoint, generative programming [START_REF] Czarnecki | Generative Programming[END_REF] provides a variety of approaches and techniques to produce and optimize a DSL implementation such as the marshaling layer in the XDR case.

Outline. Section 2 discusses how generative tools can be used to compile DSL programs into GPL programs, from a DSL interpreter. When compiled into a GPL, a DSL program can be processed by existing generative tools for various purposes, including optimization, instrumentation and verication. In this context, the generative tools are driven by domain-specic information that is translated into dierent forms: declarations (Section 3), annotations (Section 4), and meta-programs (Section 5). In essence, these forms enable a DSL to be interfaced with existing generative tools. High-Level Compilation Some level of compilation can be achieved by specializing an interpreter with respect to a DSL program. Traditionally, this approach has been used to generate compilers from denotational-style language denitions [START_REF] Consel | Tutorial notes on partial evaluation[END_REF][START_REF] Jones | Partial Evaluation and Automatic Program Generation[END_REF]. This approach was later promoted by Consel and Marlet in the context of DSLs, where the language user base often forbids major compiler development. Furthermore, the high-level nature of DSLs facilitates the introduction of optimizations.

An example of such a compilation strategy was used for a DSL, named Plan-P, aimed to specify application-specic protocols (e.g., stream-specic degradation policies) to be deployed on programmable routers [START_REF] Thibault | Safe and ecient active network programming[END_REF]. Program specialization was used at run time to achieve the eect of a Just In Time compiler, by specializing the Plan-P interpreter with respect to a Plan-P program. The resulting compiled programs ran up to 50 times faster than their interpreted counterparts [START_REF] Thibault | Safe and ecient active network programming[END_REF]. Importantly, such late compilation process enabled the safety and security of programs to be checked at the source level by the programmable routers, before being deployed. The use of specialization allowed to achieve late compilation without requiring any specic development, besides writing the interpreter.

At a lower level, meta-programming provides an alternative approach to deriving a compiler from an interpreter [START_REF] Taha | A Gentle Introduction to Multi-stage Programming[END_REF][START_REF] Czarnecki | DSL Implementation in MetaOCaml, Template Haskell, and C++[END_REF]. This approach involves a careful annotation, and sometimes re-structuring, of the interpreter. A key feature of the DSL approach is to make domain-specic information an integral part of the programming paradigm. As such the programmer can be viewed as being prompted by the language to provide domain-specic information. This information may take the form of domain-specic types, syntactic constructs and notations. This language enrichment over GPLs is typically geared towards collecting sucient information to make some domain-specic properties decidable [9] and thus to enable domain-specic verications and optimizations. The collection of information may be achieved by dedicated program analyses, which are, by design of the DSL, simpler than the ones developed for GPLs.

Because the scope of computations to be expressed by a DSL is usually narrow, GPL constructs and operations are restricted or excluded. Furthermore this language narrowing may also be necessary to enable key properties to be statically determined. In fact, a DSL is commonly both a restricted and an enriched version of a GPL.

Once key properties are exhibited, the DSL program can be compiled into a GPL program. To retain domain-specic information, the generated program needs to be accompanied by some form of declarations specifying its properties. Of course, the declarations are tailored to a set of verication and/or optimization tools.

In the XDR case, an XDR description often denes xed-size RPC arguments. When this description is compiled into GPL code, it can be accompanied by declarations aimed to drive some transformation tool. This situation is illustrated by the XDR compiler that conventionally generates C code gluing calls to a generic library. We have modied the XDR compiler to generate binding-time information about the RPC argument sizes, besides the marshaling code. As a result, a declaration is attached to the marshaling code of each data item to be transmitted. This declaration denes the size parameter as static, if it corresponds to a data item that has a xed size; it is dynamic otherwise 1 . In this work, the generated declarations are targeted for a program specializer for C, named Tempo [10]. This tool performs a number of optimizations on both the marshaling code and the generic XDR library, including the removal of buer overow checks and the collapsing of function layers.

The XDR case is interesting because it demonstrates that, although a DSL introduces a new programming paradigm, it can still converge with and re-use GPL technology. Additionally, because GPL tools are often used in a narrow context, the results may be more predictable. In the XDR case for instance, since the compilation schemas, the XDR library, and the specialization contexts are xed, specialization is fully predictable. This situation obviously does not exist for an arbitrary program with an arbitrary specialization context.

Aspect-oriented declarations could also be generated from a DSL program. For example, one could imagine adding a data compression phase to marshaling methods when invoked with data greater than a given size. In this case, the pointcut language has to be expressive enough to enable any program point of interest to be used to insert the compression phase. Alternatively, annotations can be directly injected in the generated program. From a DSL Program to Annotations Instead of generating a GPL program together with declarations, annotations can be directly inserted into the generated program. This strategy allows information to be accurately placed in the program. Like declarations, these annotations are geared towards specic tools. They can either be processed at compile time or run time.

At compile time, annotations can be used to guide the compilation process towards improving code quality or code safety. In the XDR case, for example, one could imagine a library where there would be two sets of buer operations, with or without overow checks. The selection of an operation would depend on annotations inserted in the program.

At run time, annotations can trigger specic actions upon run-time values. For instance, a data of an unknown size could be tested before being transmitted and be compressed if it is larger than a given threshold. [START_REF] Jones | Partial Evaluation and Automatic Program Generation[END_REF] From a DSL Program to Meta-programming A DSL program can also benet from meta-programming technology. In this context, the compilation of a DSL program produces a code transformer and code fragments. For example, in a multi-stage language like MetaOCaml [START_REF] Taha | A Gentle Introduction to Multi-stage Programming[END_REF], the code transformer consists of language extensions that enable to concisely express program transformations.

In the XDR case, multi-stage programming would amount to generate code that expects some data size and produces code optimized for that size. Like program specialization, multi-stage programming corresponds to GPL tools that are directly used by a programmer. In the context of DSLs, these tools can implement a specic phase of an application generator.

Conclusion

We examined generative programming approaches and techniques from a DSL viewpoint. We showed that a DSL can make use of these approaches and techniques in very eective ways. In essence, the DSL approach exposes information about programs that can be mapped into the realm of generative programming and be translated into declarations or annotations, which would normally be provided by a programmer. This situation illustrates the high-level nature of the DSL approach.

2

 2

3

  From a DSL Program to Declarations

4

 4

  Logic Programming. Number 1490 in Lecture Notes in Computer Science, Pisa, Italy (1998) 170194 10. Consel, C., Lawall, J., Le Meur, A.F.: A tour of Tempo: A program specializer for the C language. Science of Computer Programming (2004) 11. Le Meur, A.F., Lawall, J., Consel, C.: Specialization scenarios: A pragmatic approach to declaring program specialization. Higher-Order and Symbolic Computation 17 (2004) 4792

Note that these declarations go beyond data sizes. A detailed description of the language of specialization declarations and its application to the XDR example can be found elsewhere[10, 11].


Keywords: distance, HFF-AVRP distance, HFF-AVRP

The Heterogeneous Fixed Fleet Asymmetric Vehicle Routing Problem (HFF-AVRP) is a N P-hard optimization problem. Instances analysis and in particular, fitness landscape analysis, may help problem solving. Such analysis require the definition of a distance between feasible solutions. Such a distance does not exist for the HFF-AVRP and this report aims at proposing a new distance measure defined from the exchange operator. In order to compute the exchange-distance between two solutions, four algorithms are suggested and then experimented. One of them is proved to be robust and to give the exact distance whereas others only compute an upper bound.

Introduction

The Capacited Vehicle Routing Problem (CVRP) defines a category of problems often met by transportation companies. It consists in satisfying a set of customer demands thanks to a fleet of vehicles subject to limited capacity, trying to minimize costs due to road travelled. The Heterogenous Fixed Fleet Asymmetric Vehicle Routing Problem (HFF-AVRP) has been proposed [START_REF] Taillard | A heuristic column generation method for the heterogeneous vrp[END_REF] as a realistic extension of the basic Vehicle Routing Problem where a fleet with a fixed number of vehicles of various types with different capacities and variable unit running costs is considered.

The HFF-AVRP may be defined as follows: let G = (V, A) be a complete directed graph, where V = (v 0 , v 1 , v 2 , . . . , v nc ) is the node set and A = {(v i , v j ) : v i , v j ∈ V } is the arc set. Vertex v 0 corresponds to the depot and vertices v 1 to v nc correspond to the nc customers. A fleet of heterogeneous (different types) vehicles (total of nv vehicles) uses the depot as a starting base. Let us denote by c t the variable cost per distance unit of a vehicle of type t and by Q t its capacity. The number of vehicles of each type is limited (fixed). Each vehicle can be allocated at most to one route. A non-negative demand q i is associated with vertex v i (i = 1..nc) corresponding to the customers. A non-negative distance d ij (d ii = 0 for all i = 1..nc) is associated with each arc (v i , v j ) ∈ A, v i = v j , of the graph G and there exists at least two different i, j in 1..nc such that d ij = d ji (the distance matrix is asymmetric). This problem aims at determining a set of vehicle routes, each starting and ending at the depot, such that a single vehicle supplies each customers demand, the total demand of customers assigned to each route does not exceed the capacity of the vehicle assigned to it, and the total cost is minimized.

HFF-AVRP is a N P-hard combinatorial optimization problem where it is very difficult for any resolution method to find a global optimum. With each HFF-AVRP instance, a structure, named landscape [START_REF] Wright | The roles of mutation, inbreeding, crossbreeding and selection in evolution[END_REF], can be associated. It corresponds to a triplet (Ω, N , f ) where Ω is a search space, N is a neighborhood definition and f a fitness evaluation. Landscape analysis helps to find particularities of instances which may have influences on the resolution and so give appropriate information to the decision maker.

Some landscape analysis studies [START_REF] Bachelet | Métaheuristiques parallèles hybrides : application au problème d'affectation quadratique[END_REF][START_REF] Schiavinotto | A review of metrics on permutations for search landscape analysis[END_REF][START_REF] Verel | The connectivity of NK landscapes' basins: a network analysis[END_REF] have been made for different optimization problems, and they use several indicators to find interesting properties. A lot of indicators require to define a distance between feasible solutions. This distance should evaluate the number of elementary moves required to move from one solution to another. Hence, the distance may be different regarding the operator used. So, in order to make a landscape analysis of HFF-AVRP instances, a distance between solutions is required. As far we know, up to now, no exact distances have been proposed for such problems. Our study aims at filling this gap giving a reliable distance between HFF-AVRP solutions as well as an algorithm to compute it efficiently. Section 2 proposes a representation and a neighborhood operator to HFF-AVRP solutions. Section 3 presents the distance and several algorithms to compute it. In section 4, experimentations of the different proposed algorithms are lead. Finally, section 5 draws a conclusion and exposes perspectives for this work.

An HFF-AVRP structure

The HFF-AVRP is less studied than the symmetric CVRP with a homogeneous fleet. The dedicated literature (for example [START_REF] Laporte | An exact algorithm for the asymmetrical capacitated vehicle routing problem[END_REF][START_REF] Toth | An overview of vehicle routing problems[END_REF]) does not deal with neither solution representation nor neighborhood definition. As studies on CVRP are numerous, they give a reliable basis for an HFF-AVRP study. Figure 1 presents an HFF-AVRP solution with 7 customers and 4 vehicles. This solution contains four routes: the first one, allocated to vehicle 1, visits customer 1, then customer 3 and finishes with customer 2. The second route, allocated to vehicle 2, visits customer 4, then visits customer 7 to finish with customer 5. The thirs one, allocated to vehicle 3, only visits customer 6. The last one is empty. This figure illustrates the key concepts of the current section.

State of the art on solution representation

Solving a combinatorial optimization problem highly depends on how it has been modelled. For CVRP, different representations may be found in the literature. In 1983, Beasley [START_REF] Beasley | Route-first cluster-second methods for vehicle routing[END_REF] introduces the big tour representation. Customers are listed in the order of their visit, and then a split procedure forms the routes : it is route first / cluster second principle. Thus, a solution corresponds to the route list in which customers are organised by vehicle visit order. For example, the HFF-AVRP solution of Figure 1 could be first represented as 1 3 2 4 7 5 6 and then the split procedure would give the representation 1 3 2, 4 7 5, 6, ∅.

Kubiak [START_REF] Kubiak | Distance measures and fitness-distance analysis for the capacitated vehicle routing problem[END_REF] prefers to consider a set of routes where each route is itself an ordered set of customers. This representation is the most used because it allows to easily access any route and then any customer. For example, the HFF-AVRP solution of Figure 1 would be represented as

[[1 3 2], [4 7 5], [6], []].
Another representation [START_REF] Michalewicz | How to Solve It: Modern Heuristics[END_REF], named adjacency representation, put the customers in a vector whose size corresponds to the number of customers. Each index (numbered from 1 to nc) represents a customer. To each index corresponds its successor. When a customer finishes a route, the vector is filled with a zero at the correponding index. For example, the HFF-AVRP solution of

Distance and neighborhood operator

In the literature, three different ways to compute the distance between CVRP solutions exist (these distances are detailed in part 3.1.):

1. distance with structural comparisons [START_REF] Kubiak | Distance measures and fitness-distance analysis for the capacitated vehicle routing problem[END_REF], 2. approximative distances such as Hamming distance [START_REF] Ho | Path relinking for the vehicle routing problem[END_REF] and pseudo-distance [START_REF] Prins | Solving the capacitated locationrouting problem by a GRASP complemented by a learning process and a path relinking[END_REF],

3. distance defined by an operator [START_REF] Sörensen | Path relinking for the vehicle routing problem using the edit distance[END_REF].

In a context of landscape analysis, the third type of distance is particularly interesting as the concepts of neighborhood and distance have to be linked.

The exchange operator (also called interchange operator) is a classical operator for permutation problems (see Schiavinotto and Stützle [START_REF] Schiavinotto | A review of metrics on permutations for search landscape analysis[END_REF]). It consists in switching two different clients. For a permutation of size N , this operator (i, j) exchanges element i with element j. In the traveling salesman problem (TSP), this operator is often used. Since a VRP can be viewed as a multiple TSP, we choose to build the neighborhood structure using the exchange operator.

A solution representation adapted to a neighborhood operator

Even if representations of the literature are widely used for different approaches of CVRP resolution, they are maladaptive to work with the HFF-AVRP and the exchange operator. Therefore, the chosen representation for this study consists in a list of customers (numbered between 1 and nc, the number of customers) ranked following their visit order and separated by zeros meaning the end of a route and the start of another. For readability, a solution starts and finishes with a zero. In the following, indexation is from 0 to nc + nv -1. For example, the HFF-AVRP solution of Figure 1 could be expressed as:

(0 1 3 2 0 4 7 5 0 6 0 0)

With this representation, it is easy to know in which route a customer is and its visiting order.

This representation suits to the exchange operator which consists in exchanging either two customers or one customer and a zero (exchanging two zeros is forbidden otherwise the solution does not change). So, with the previous example, the exchange operator (1, 7) acts on customers located at indexes 1 and 7 and gives the solution:

(0 5 3 2 0 4 7 1 0 6 0 0)

Notice that this solution can become unfeasible after an exchange if the capacity constraints of vehicles are not respected. For this representation and its associated operator, we aim to define an associated distance.

3 An adaptated distance

State of the art

Only few studies exist on HFF-AVRP, so we have studied distances proposed for the CVRP. These distances can be divided into three groups:

RR n°7263

Comparative distance Kubiak, in a fitness-distance analysis [START_REF] Kubiak | Distance measures and fitness-distance analysis for the capacitated vehicle routing problem[END_REF] uses three distances that can be described as comparative. Indeed, he measures the similarities between two solutions by examining either if customers have the same predecessor and successor in the two solutions, or if each route of the two solutions visits the same customers. Even if these distances are not based on a specific operator, they can help to know how dissimilar or close two CVRP solutions are.

Approximate distance Hamming distance is a well-known distance in combinatorial optimization. a path-relinking study, Ho and Gendreau [START_REF] Ho | Path relinking for the vehicle routing problem[END_REF] adapt this distance to CVRP solutions by numbering all edges (an edge is the road between customer/vehicle to its successor) and creating a bit-vector. Hamming distance corresponds to the number of different bits between two solutions. Another approximate distance is proposed by Prins and al. [START_REF] Prins | Solving the capacitated locationrouting problem by a GRASP complemented by a learning process and a path relinking[END_REF] in another pathrelinking study. This distance consists in giving a weight to each edge of the first solution according to the similarity or the difference in the second solution.

Operator related distance In a study on path relinking for the CVRP, Sörensen [START_REF] Sörensen | Path relinking for the vehicle routing problem using the edit distance[END_REF] presents the edit-distance based on three elementary operators: add, remove and substitute. The distance between two feasible solutions is the minimal number of operator applications required to move from the first solution to the second one. Then, he tranforms this distance keeping only add and remove operators. Since solutions are represented as a big-tour, this distance corresponds to the insert distance for permutation, defined in [START_REF] Schiavinotto | A review of metrics on permutations for search landscape analysis[END_REF]. This distance is an interesting proposition of a distance between CVRP solutions defined from an operator.

Regarding all these distances, we see that comparative and approximate distances give information about closeness or dissimilarities between solutions, their edges and routes, whereas distances based on operators give a measure that can be considered through the search space with a neighborhood structure. As explained before, to make a landscape analysis, an operator based distance is needed.

The proposed exchange-distance

The interest of having a distance associated to an operator is that the number of modifications required to move from a solution to another can be measured and evaluated. In this study, we aim at providing a distance associated to the exchange operator. Our proposed distance, called exchange-distance d e , is defined by the minimal numbers of exchange operators to be applied to move from an HFF-AVRP solution to another. This definition verifies the mathematical conditions:

• positive definiteness is clear,

d e ≥ 0 ∀(s 1 , s 2 ), d e (s 1 , s 2 ) = 0 ⇔ s 1 = s 2 RR n°7263
• the symmetry property comes from the symmetry property of the exchange operator ∀(s 1 , s 2 ), d e (s 1 , s 2 ) = d e (s 2 , s 1 )

• and the triangle inequality is given by the minimal caracteristic of the definition ∀(s 1 , s 2 , s 3 ),

d e (s 1 , s 2 ) ≤ d e (s 1 , s 3 ) + d e (s 3 , s 2 )
Indeed, if p exchange operator applications are needed to move from s 1 to s 3 and q exchange operator applications are needed to move from s 3 to s 2 , then it is possible to move from s 1 to s 2 by at most p+ q exchange operator application.

With this distance, the neighborhood of an HFF-AVRP solution defined with the exchange operator can be described as the set of all the HFF-AVRP solutions which are at a distance equal to one. In the same way, the n-neighborhood of an HFF-AVRP solution is the set of HFF-AVRP solutions which are at a distance equal to n from the first one. In the case of an HFF-AVRP instance with nc customers and nv vehicles, we can proove that the maximal neighborhood is nc-order (it does not depend on the number of vehicles).

Algorithms to compute the exchange-distance

Defining and computing a distance is the aim of this study. So, an algorithm which computes the exchange-distance between HFF-AVRP solutions s 1 and s 2 is required. Obviously, this algorithm depends on the HFF-AVRP representation. In the following, we present four algorithms of different complexity and accuracy.

Permutation-adapted algorithm As the representation looks like a permutation, we first try to use an existing algorithm which calculates the exchangedistance in the case of permutations. Bachelet [START_REF] Bachelet | Métaheuristiques parallèles hybrides : application au problème d'affectation quadratique[END_REF] proposed an O( N ) algorithm for such purpose (N is the permutation size). This algorithm, based on the inverse operator of the exchange operator for a permutation, has to be adapted to our representation. Since zeros appear several times in our representation, we decide to number vehicles from nc + 1 to nc + nv in order to be able to use Bachelet algorithm. Thus, we obtain a real permutation. This adapted algorithm, called permutation-adapted, is fast. Unfortunately having differenciated vehicles leads to wrong computed distance. Even if it gives an upper bound of the real value of the distance, a robust calculation is needed. We target to find another one. Three algorithms emerge from this work.

Simple greedy algorithm This second strategy is greedy i.e. solution s 1 is scanned from left to right in order to transform it into s 2 with the minimal number of applications of the exchange operator. At the end of the p-step, the p first elements (customers and vehicles) in the two considered solutions are the same. So, at each step, if the values in s 1 and s 2 are different, the value of s 2 is sought through s 1 from index p + 1. And then, the exchange operator is called. This algorithm computes an upper bound as a feasible sequence of exchange operators is produced.

With-choice greedy algorithm The simple greedy algorithm has a major drawback: it does not optimize the replacement of zeros. A better strategy needs to be defined. If a zero has to be replaced, instead of choosing the first one in s 1 , the idea is to try to find a zero which puts the element p at its position . This new algorithm is more efficient than the previous one when the exchange operator puts both vehicle and customer in their right place. The calculation of the distance becomes more precise i.e. the value of the upper bound is ameliorated.

On-hold greedy algorithm Ordering solution at each step introduces errors which increase the distance. We decided to break this approach by ordering all the customers first. We realize all the zeros are then well placed back. In order to have an optimized algorithm, a memory with a list of indexes is created. When a zero has to be replaced and when any exchange allows the current customer to be replaced at his right position, the index of the step is put in the memory. Hence, the memory contains all the indexes ever scanned where s 1 values are not equal to s 2 ones. This list has to be scanned at each step because a customer refered by it can be found. Using a classical exchange procedure, the on-hold greedy algorithm is detailed in Algorithm 1.

The following example illustrates the algorithm. To move from s 1 to s 2 , four exchange operations are needed. Let us see the detailed execution of the algorithm:

s 1 : (0 1 3 2 7 0 6 0 4 5 0) s 2 : (0 1 3 2 0 4 7 5 0 6 0)

• i = 1 to 3 s 1i = s 2i : the values are equal in both solutions. distance = 0, list = ().

• i = 4 s 1i = s 2i : 0 is needed. list is empty. s 2i = 0: no 0 at the right position of 7. Index 4 is put in the memory list. distance = 0, list = (4).

• i = 5 s 1i = s 2i : 4 is needed. It is not refered by the list.

4 is found at index 8. Exchange (5, 8) gives (0 1 3 2 7 4 6 0 0 5 0). distance = 1, list = (4).

• i = 6 s 1i = s 2i : 7 is needed. It is refered by the list. s 1i = 0: the list remain unexchange. Exchange (4, 6) gives (0 1 3 2 6 4 7 0 0 5 0). distance = 2, list = (4).

• i = 7 s 1i = s 2i : 5 is needed. It is not refered by the list.

5 is found at index 9. Exchange (7, 9) gives (0 1 3 2 6 4 7 5 0 0 0). distance = 3, list = (4).

• i = 8 s 1i = s 2i : the values are equal in both solutions. distance = 3, list = (4).

• i = 9 s 1i = s 2i : 6 is needed. It is refered by the list. s 1i = 0: the index which refers 6 is removed from the list. Exchange (4, 9) gives (0 1 3 2 0 4 7 5 0 6 0). distance = 4, list = ().

This algorithm needs more conditions than the other two greedy algorithms. Experiments reported later show it always finds the exact distance. On the algorithmic point of view, the three proposed greedy algorithms are more and more difficult to implement. Indeed, more and more conditions are RR n°7263 defined in order to be closer to the exact distance. On the complexity point of view, they are equivalent. Indeed, they have on average an O( N 2 ) complexity (N is the length of the solution). Therefore, the on-hold greedy algorithm has to be considered to compute the distance between two HFF-AVRP solutions, because it gives the exact distance and it is not more time consuming than the other two.

Experimentations

The permutation-adapted algorithm and the three greedy algorithms are tested on the proposed HFF-AVRP representation in order to find cases where the calculation of the distance between two solutions is exact. Here, a solution is represented by a vector of N elements with nc customers and nv + 1 zeros (as previously said, nv is the number of vehicles). Note, the number of vehicles taken into account is smaller than half the number of customers, to avoid as possible, routes with only one customer to visit.

Implementation

Tests consist in randomly creating solutions at a d chosen distance from a random initial solution. Then each algorithm computes the distance between the initial and the final solutions to compare it with the exact value (d).

Initialization To generate HFF-AVRP solutions, every customers are ranked in the numerical order (from 1 to nc) in a vector. Then, two zeros are added at the beginning and at the end of the vector. They cannot be exchanged because the representation has to start and to finish by zero. To have nv routes, nv -1 zeros are randomly inserted in the vector in such a way that no route is empty (zero can not be closed).

Successive neighborhood From this initial generated solution, the exchange operator is applied d times to get a solution which belongs to its d-neighborhood. Remember the maximal neighborhood of a solution is the number of customers (nc). Two cases can be distinguished to create a d-neighborhood. Indeed, if d is equal to nc, as in permutation problems, a cycle is sufficient to have a solution in the d-neighborhood (order d). Since no route is empty, it is impossible to have a zero which replaces another. This is why we are sure of the existence of a ncneighbor. For all d-neighborhood with d strictly lower than nc, a list is created with all the indexes between 2 and (N -1) corresponding to the elements which could be exchanged. Two indexes are chosen in this list (corresponding to at most one zero). If the indexes correspond to two customers then the first index is removed from the list. Thus, during next steps, the first element will not be put in its first position because it does not belong to the list anymore and so for the second element, because it cannot move anymore to its first position. If the two chosen indexes corresponds to a vehicle and a customer, then the two indexes are removed from the list: for sure any zero will be able to substitute another. At the end of step d = nc -1, the list is empty.

Experiments The four algorithms are tested using several parameters:

• number of customers,

• number of vehicles,

• distance of the neighborhood. Indeed, the larger is the neighborhood and the more the number of customers and/or vehicles are/is huge, the bigger is the probability to have a wrong computed distance. Algorithms are tested with 10, 100 and 1000 customers. For each case the number of vehicles equals: 1 vehicle (case of the traveling salesman problem), then 5%, 20% and 50% of the number of customers. The nc/8, nb/4, nc/2, nc -1 and nc -neighborhood are computed for each case.

Analysis

To realize the analysis over the thousand computed values, two indicators are used. The first one is the error percent (1) which gives information about the difference between a theoretical value and its computed one. The error percent for each case could be computed as:

Err% = d -D D ( 1
)
where d is the computed distance and D is the theoretical distance. As the four algorithms compute an upper bound, d is positive and so d -D, the average of the difference, is meaningful. The second indicator is the coefficient of variation of a data series (2) which is dimensionless and gives information about dispersion of a distribution in comparison to its mean. The coefficient of variation is computed as:

∆ = σ d ( 2
)
where σ is the standard deviation. Table 2 presents the error percents for each algorithm and the computed distances. The values are given in three different vertical groups to show the efficiency of each algorithm depending on the number of customers. Moreover, the values are split into three horizontal groups in order to show the influence of the number of vehicles.

nc 0 0 0 0 3 0 0 0 4.8 0 0 0 (0) (0) (0) (0) (0) (0) (0) (0) (0) (0) (0) (0)
(0) (0) (0) (0) (0) (0) (0) (0) (0) (0) (0) (0)
Table 2: Error percents of the distance and coefficients of variation for each measure of the experimentation in brackets. Algorithm 1 stands for permutationadapted algorithm, 2 for simple greedy algorithm, 3 for with-choice greedy algorithm and 4 for on-hold greedy algorithm.

RR n°7263

Undoubtedly, the on-hold greedy algorithm can be considered as robust. Indeed, the error between the exact and the computed values is always null. Besides, Table 2 presents in brackets the coefficients of variation which give information about the behaviour of the other three algorithms. Indeed, regarding only error percents, the on-hold greedy algorithm is the best and the with-choice greedy algorithm is better than the simple greedy algorithm which is better than the permutation-adapted algorithm.

For these three last algorithms, the larger the number of customers and/or the number of vehicles are, the higher the error percent is. The permutationadapted algorithm and the simple greedy algorithm tend to be confused as soon as these numbers increase whereas the with-choice greedy algorithm seems not to be perturbed too much when the number of customers increases because values remain low for 100 or 1000 customers. However, the coefficient of variation gives an interesting information about how different or not the thousand values for each case of the experimentation are. Even if the results of the permutationadapted algorithm and the simple greedy algorithm are closed and often bad, the coefficient of variation is always greater for the last one. The chance to give a better value is huge. Also, we can see that the with-choice greedy algorithm has got, in the case of 100 and 1000 customers, a coefficient of variation greater than the permutation-adapted algorithm and the simple greedy algorithm ones; it shows that this algorithm, even if is better, the chance to get a wrong value is sizeable.

Let us remark, for 1-vehicle (i.e. travelling salesman problem), all algorithms are robust, and so, permutation-adapted algorithm is the fastest one.

Conclusion and perspective

This study aims to fill a gap by defining a distance between HFF-AVRP solutions and an algorithm to compute it. Therefore, a solution representation and a neighborhood operator are needed. A representation close to the permutation representation and the exchange operator were chosen. Then a distance, called the exchange-distance, was defined by the exchange operator as the minimal number of exchanges necessary to move from an HFF-AVRP solution to another. With this theoretical definition of the exchange-distance, a robust algorithm was searched to compute it. Four algorithms were suggested and tested. Experimentations showed that the on-hold greedy algorithm computes the exact distance, and so it could be used in order to work on distance between the HFF-AVRP solutions. Since this study gives a reliable distance based on the exchange operator between feasible HFF-AVRP solutions, the next step is to realize a landscape analysis using it.

Figure 1 :

 1 Figure 1: Example of an HFF-AVRP solution with 7 customers and 3 vehicles.

5 ]

 5 (indexes are added to help the understanding).

  end if else {It exists j such that s 1j = s 2i i.e. s 2i = 0 and j ∈ list}

	temp ← j
	end if
	end if
	if temp exists then
	Exchange(s 1 ,i,temp)
	distance ← distance + 1
	end if
	end if
	end for

Algorithm 1 On-hold greedy algorithm Input: s 1 , s 2 , N Output: ExchangeDist(s 1 ,s 2 ) distance ← 0 list ← () temp {variable to keep the index if exchange is possible} for i = 1 to N do {where N is the length of the solutions s 1 and s 2 } if s 1i = s 2i then if list is not empty and it exists j such that s 1listj = s 2i then {s 2i is refered in list} temp ← list j if s 1i = 0 then {0 finds its positions, index has to be removed from the list} list ← list \ list j end if else if s 2i = 0 then if It exists j such that s 1j = 0 and s 1i = s 2j then {s 1i can be directly well repositioned} temp ← j else list ← list :: i

Table 1

 1 reports the chosen values for the tests.

	Number of	Number of	Neighborhood
	Customers	Vehicles	order
	10	1, 2, 5	1, 2, 5, 9, 10
	100	1, 5, 20, 50	12, 25, 50, 99, 100
	1000	1, 50, 200, 500 125, 250, 500, 999, 1000

Table 1 :

 1 Tested values.

	Each possibility has been tested 1000 times.

RR n°7263


D Basic
email: duro.basic@schneider-electric.com

F Malrait

Pierre Rouchon
email: pierre.rouchon@mines-paristech.fr

Euler-Lagrange models with complex currents of three-phase electrical machines and observability issues

Keywords: Lagrangian with complex coordinates, space-harmonics, magnetic saturation, induction machine, permanent-magnet machine, sensorless control

A new Lagrangian formulation with complex currents is developed and yields a direct and simple method for modeling three-phase permanentmagnet and induction machines. The Lagrangian is the sum a mechanical one and of a magnetic one. This magnetic Lagrangian is expressed in terms of rotor angle, complex stator and rotor currents. A complexification procedure widely used in quantum electrodynamic is applied here in order to derive the Euler-Lagrange equations with complex stator and rotor currents. Such complexification process avoids the usual separation into real and imaginary parts and simplifies notably the calculations. Via simple modifications of such magnetic Lagrangians we derive new dynamical models describing permanent-magnet machines with both saturation and saliency, and induction machines with both magnetic saturation and space harmonics. For each model we also provide its Hamiltonian thus its magnetic energy. This energy is also expressed with complex currents and can be directly used in Lyapunov and/or passivity based control. Further, we briefly investigate the observability of this class of Euler-Lagrange models, in the so-called sensorless case when the measured output is the stator current and the load torque is constant but unknown. For all the dynamical models obtained via such variational principles, we prove that their linear tangent systems are unobservable around a one-dimensional family of steady-states attached to the same constant stator voltage and current. This negative result explains why sensorless control of three-phase electrical machines around zero stator frequency remains yet a difficult control problem.

Introduction

Modeling electrical machines with magnetic-saturation and space-harmonics effects is not a straightforward task and could lead to complicated developments when a detailed physical description is included (see, e.g., [START_REF] Chiasson | Modeling and High Performance Control of Electric Machines[END_REF][START_REF] Boldea | The Induction Machine Handbook[END_REF]). Even if such effects are not dominant they could play an important role for sensorless control (no rotor position or velocity sensor). For a permanent-magnet machine, the rotor position will be unobservable without saliency. For standard models of induction machines (no magnetic saturation, no space-harmonics) the rotor velocity is always unobservable at zero stator frequency [START_REF] Canudas | Observability conditions of induction motors at low frequencies[END_REF][START_REF] Holtz | Sensorless control of induction motor drives[END_REF]. A global observability analysis based on such standard models is given in [START_REF] Ibarra-Rojas | Global observability analysis of sensorless induction motors[END_REF]. In this note we develop a systematic method to include into such standard control models magnetic-saturation, saliency and/or space-harmonics effects. Our initial motivation was to see whether such non-observability is destroyed by such modelling changes or not. It appears that any physical-consistent model admits the same kind of observability deficiency at zero stator frequency.

By physically consistent models we mean Lagrangian-based models. We propose here an extension of Lagrangian modeling of three-phase machines with real variables (see, e.g. [START_REF] Ortega | Passivity-Based Control of Euler-Lagrange Systems[END_REF]) to complex electrical variables. It is directly inspired from quantum electro-dynamics where Lagrangian with complex generalized positions and velocities are widely used (see, e.g., [START_REF] Cohen-Tannoudji | Photons and Atoms: Introduction to Quantum Electrodynamics[END_REF], page 87). We obtain, from such Lagrangian functions, physically consistent and synthetic Euler-Lagrange models directly expressed with complex stator and rotor currents. Such modeling method by-passes the usual detailed physical descriptions that are not easily accessible to the control community. Here we propose a much more direct way: it just consists in modifying the magnetic part of the Lagrangian directly expressed with complex currents and then in deriving the dynamic equations from the Euler-Lagrange with complex variables. We obtain automatically the dynamics of the electrical part as a set of complex differential equations. We suggest here simple Lagrangians modeling simultaneously magnetic-saturation, saliency and space-harmonic effects. The obtained dynamics extend directly the ones used in almost all control-theoretic papers and include also more elaborate ones that can be found in specialized books such as [START_REF] Boldea | The Induction Machine Handbook[END_REF].

For permanent-magnet three-phase machines, the general structure of any physically consistent model including magnetic-saturation, saliency and other conservative effects is given by equations ( 4) with magnetic Lagrangian L m and Hamiltonian H m related by [START_REF] Fudeh | Modeling and analysis of induction machines containing space harmonics. Part-I: modeling and transformation[END_REF]. For induction three-phase machines the physically consistent models are given by (15) where the magnetic Lagrangian L m is related to the magnetic energy H m by (16). Such synthetic formulation of the dynamical equations is new and constitutes the first contribution of this note. We propose here natural modifications of the standard Lagrangian to include magnetic-saturation, saliency and space-harmonics, derive the corresponding dynamical equations and magnetic energies that could be used in the future to construct controlled Lyapunov functions and/or storage functions for passivitybased feedback laws.

From a control theoretic point of view we just prove here that the severe ob-servability difficulties encountered in sensorless control and well explained in [START_REF] Ibarra-Rojas | Global observability analysis of sensorless induction motors[END_REF] for the standard model resulting from the quadratic Lagrangian (13), remain present for models (15) where the magnetic Lagrangian is any function of the rotor angle, stator and rotor currents. Consequently, addition of magnetic saturations, saliency and harmonics effects, do not remove observability issues at zero stator frequency in the sensorless case (see proposition 1). This observability obstruction has nether been stated for models with magnetic-saturation and space-harmonics of three-phase machines and constitutes the second contribution of this notes. Contrarily to observability, non-observability is not a generic property and could be destroyed by generic and small changes in the equations.

Since proposition 1 is based on the class of models derived from ( 1) or (15) with arbitrary magnetic Lagrangian L m , we prove here that any physically consistent model of three-phase machines where the non-conservative effects result only from voltage supply and Ohmic losses, such non-observability holds true around zero stator frequency. This means that non-observability around zero stator frequency is robust to generic and physically consistent modifications of the equations. As far as we know this negative and physically robust result is new. It indicates that sensorless control of three-phase electrical machines around zero stator frequency cannot be just addressed via refined physical models but also requires advanced and nonlinear control techniques.

In section 2 we recall the simplest model of a permanent-magnet machine and its Euler-Lagrange formulation based on the two scalar components of the complex stator current. Using the complexification procedure detailed in appendix, we show how to use complex representation of stator-current in Lagrangian formulation of the dynamics. This leads us to the general form of physically consistent models (equation ( 4)). Finally we obtain, just by simple modifications of the magnetic Lagrangian, physically consistent models with magnetic saturation and saliency effects (equation [START_REF] Yourgrau | Variational Principles in Dynamics and Quantum Theory[END_REF]). Section 3 deals with induction machines and admits the same progression as the previous one: we start with the usual (α, β) model, describe its complex Lagrangian formulation, derive physically consistent models (equations (15)) and specialize them to saturation and space-harmonics effects (equation (21) ). In section 4, we prove proposition 1 that states the main observability issues of these Euler-Lagrangian models at zero stator frequency. In conclusion we show how to transpose this modelling based on complex currents associated to a Lagrangian formulation to complex fluxes associated to a Hamiltionian formulation with complex generalized positions and momentums. The appendix details the complexification procedure. It explains how to derive the Euler-Lagrange equations when some generalized positions and velocities are treated as complex quantities. Throughout the paper, we define models in (α, β) frame, using the standard transformation from three phases frame (see, e.g., [START_REF] Leonhard | Control of Electrical Drives[END_REF]).

2 Permanent-magnet three-phase machines 2.1 The usual model and its magnetic energy

In the (α, β) frame (total power invariant transformation), the dynamic equations read (see, e.g., [START_REF] Leonhard | Control of Electrical Drives[END_REF]):

     d dt J θ = n p ℑ φe npθ * ı s -τ L d dt λı s + φe npθ = u s -R s ı s (1)
where

• * stands for complex-conjugation, ℑ means imaginary part,  = √ -1 and n p is the number of pairs of poles.

• θ is the rotor mechanical angle, J and τ L are the inertia and load torque, respectively.

• ı s ∈ C is the stator current, u s ∈ C the stator voltage.

• λ = (L d + L q )/2 with inductances L d = L q > 0 (no saliency here).

• The stator flux is φ s = λı s + φe npθ with the constant φ > 0 representing to the rotor flux due to permanent magnets.

The Lagrangian associated to this system is the sum of the mechanical one L c and magnetic one L m defined as follows:

L c = J 2 θ2 , L m = λ 2 ı s + īe npθ 2 (2)
where ī = φ/λ > 0 is the permanent magnetizing current.

It is well known that (1) derives from a variational principle (see, e.g., [START_REF] Ortega | Passivity-Based Control of Euler-Lagrange Systems[END_REF]) and thus can be written as Euler-Lagrange equations with source terms corresponding to energy exchange with the environment. Consider the additional complex variable q s ∈ C defined by d dt q s = ı s . Take the Lagrangian L = L c +L m as a real function of the generalized coordinates q = (θ, q sα , q sβ ) and generalized velocities q = ( θ, ı sα , ı sβ ):

L(q, q) = J 2 θ2 + λ 2 `( qsα + ī cos npθ) 2 + ( qsβ + ī sin npθ) 2 ´(3)
Then the mechanical equation in (1) reads

d dt ∂L ∂ θ - ∂L ∂θ = -τ L
where -τ L corresponds to the energy exchange through the mechanical load torque. Similarly, the real part of complex and electrical equation in (1) reads

d dt ∂L ∂ qsα - ∂L ∂q sα = u sα -R s ı sα
and its imaginary part

d dt ∂L ∂ qsβ - ∂L ∂q sβ = u sβ -R s ı sβ
since ∂L ∂qsα = ∂L ∂q sβ = 0 and qs = ı s . The energy exchanges here are due to the power supply through the voltage u s and also to dissipation and irreversible phenomena due to stator resistance represented by the Ohm law -R s ı s .

Euler-Lagrange equation with complex current

The drawback of such Lagrangian formulation is that we have to split into real and imaginary parts the generalized complex coordinates with q s = q sα + q sβ , (q sα and q sβ real) and velocities qs = ı s = qsα +  qsβ , ( qsα and qsβ real). We do not preserve the elegant formulation of the electrical part through complex variables and equations.

Let us apply the complexification procedure detailed in appendix to the Lagrangian L(θ, q sα , q sβ , θ, qsα , qsβ ) defined in [START_REF] Chiasson | Modeling and High Performance Control of Electric Machines[END_REF]. The complexification process only focuses on q s and qs = ı s by considering L as a function of (θ, q s , q * s , θ, ı s , ı * s ):

L(θ, θ, ı s , ı * s ) = J 2 θ2 + λ 2 ı s + īe npθ ı * s + īe -npθ .
Following the notations in appendix, n c = 1 with q c = q s , n r = 1 with q r = θ, S r = -τ L and S c = u s -R s ı s . Then according to (25) the usual equations ( 1)

read d dt ∂L ∂ θ = ∂L ∂θ -τ L , 2 d dt ∂L ∂ı * s = u s -R s ı s since ∂L ∂q * s = 0 and ∂L ∂ q * s = ∂L ∂ı * s .
More generally, for any magnetic Lagrangian L m that is a real value function of θ, ı s and ı * s and that is 2π np periodic versus θ, we get the general model (with saliency, saturation, space-harmonics, ...) of three-phase permanent-magnet machine:

d dt J θ = ∂L m ∂θ -τ L , 2 d dt ∂L m ∂ı * s = u s -R s ı s (4)
We recover the usual equation with φ s = 2 ∂Lm ∂ı * s corresponding to the stator flux and to the conjugate momentum p c of q c as shown in appendix. According to (28), the Hamiltonian H is the sum of two energy: H = H c + H m . The mechanical kinetic energy H c = J 2 θ2 and the magnetic energy

H m (θ, ı s , ı * s ) = ∂L m ∂ı s ı s + ∂L m ∂ı * s ı * s -L m . (5)
The standard model (1) derives from a magnetic Lagrangian of the form L m = λ 2 ı s + īe npθ 2 with λ and ī are two positive parameters. Its corresponding magnetic energy reads H m = λ 2 |ı s | 2 -ī2 . We recover the usual magnetic energy λ 2 |ı s | 2 up to the constant magnetizing energy λ 2 ī2 . In (4), many other formulations of L m are possible and depend on particular modeling issues. Usually, the dominant part of L m will be of the form λ 2 ı s + īe npθ 2 ( λ, ī positive constants) to which correction terms that are "small" scalar functions of (θ, ı s , ı * s ) are added.

Saliency models

Adding to L m the correction -µ 2 ℜ ı 2 s e -2npθ with |µ| < λ (ℜ means real part) provides a simple way to represent saliency phenomena while the dominant part of the magnetic Lagrangian (and thus of the dynamics) remains attached to λ 2 ı s + īe npθ 2 . With a magnetic Lagrangian of the form

L m = λ 2 ı s + īe npθ ı * s + īe -npθ - µ 4 ı * s e npθ 2 + ı s e -npθ 2 (6)
where

λ = (L d + L q )/2 and µ = (L q -L d )/2 (inductances L d > 0 and L q > 0), equations (4) become (λī = φ)      d dt J θ = n p ℑ (λı * s + φe -npθ -µı s e -2npθ )ı s -τ L d dt λı s + φe npθ -µı * s e 2npθ = u s -R s ı s (7)
and we recover the usual model with saliency effect. In this case the magnetic energy is given by:

H m = λ 2 |ı s | 2 -ī2 - µ 4 ı * s e npθ 2 + ı s e -npθ 2 . (8)

Saturation and saliency models

We can also take into account magnetic saturation effects, i.e., the fact that inductances depend on the currents. Let us assume first that only the inductances λ and µ in (6) depend on the modulus ρ = ı s + īe npθ . The magnetic Lagrangian now reads

L m = λ ı s + īe npθ 2 ı s + īe npθ 2 - µ ı s + īe npθ 4 ı * s e npθ 2 + ı s e -npθ 2 . ( 9
)
The dynamics is given by ( 4

d dt " J θ" = npℑ "" Λ " ı * s + īe -npθ " -M ıse -2npθ " ıs " -τL d dt " Λ " ıs + īe npθ " -M ı * s e 2npθ " = us -Rsıs (10)
with Λ = λ+

|ıs+īe np θ | 2 λ ′ and M = µ+ |ıs+īe np θ | 2 µ ′ .
It is interesting to compute the magnetic energy H m from general formula (5):

Hm = λ + ˛ıs + īe npθ ˛λ′ 2 |ıs| 2 - λ 2 ī2 + ˛ıs + īe npθ ˛λ′ 4 ī " ıse -npθ + ı * s e npθ " - µ + µ ′ ı * s (ıs+īe np θ )+ıs(ı * s +īe -npθ ) 2|ıs+īe np θ | 4 " " ı * s e npθ " 2 + " ıse -npθ " 2 « . (11)
Such magnetic energy formulae are not straightforward. They are not obtained by replacing λ and µ in the standard magnetic energy (8) by Λ and M respectively.

Induction three-phase machines

We will now proceed as for permanent-magnet machines. Let us recall first the usual dynamical equations of an induction machine with complex stator and rotor currents. They admit the following form

             d dt J θ = n p ℑ L m ı * r e -npθ ı s -τ L d dt L m ı r + ı s e -npθ + L f r ı r = -R r ı r d dt L m ı s + ı r e npθ + L f s ı s = u s -R s ı s (12)
where

• n p is the number of pairs of poles, θ is the rotor mechanical angle, J and τ L are the inertia and load torque, respectively.

• ı r ∈ C is the rotor current (in the rotor frame, different from the (d, q) frame) , ı s ∈ C the stator current (in the stator frame, identical to the (α, β) frame) and u s ∈ C the stator voltage (in the stator frame). The stator and rotor resistances are R s > 0 and R r > 0.

• The inductances L m , L f r and L f s are positive parameters with L f r , L f s ≪ L m .

• The stator (resp. rotor) flux is φ s = L m ı s + ı r e npθ + L f s ı s (resp. φ r = L m ı r + ı s e -npθ + L f r ı r ).

Euler-Lagrange equation with complex current

With notations of appendix, n c = 2 with q c = (ı r , ı s ), n r = 1 with q r = θ, S c = (-R r ı r , u s -R s ı s ) and S r = -τ L . The Lagrangian associated to (12), expressed with complex currents ı r and ı s , reads:

L(θ, θ, ı r , ı * r , ı s , ı * s ) = J 2 θ2 + L m 2 ı s + ı r e npθ 2 + L f r 2 |ı r | 2 + L f s 2 |ı s | 2 .
The first term J 2 θ2 represents the mechanical Lagrangian and the remaining sum the magnetic Lagrangian L m . The dynamics (12) read

d dt ∂L ∂ θ - ∂L ∂θ = -τ L , 2 d dt ∂L ∂ı * r = -R r ı r , 2 d dt ∂L ∂ı * s = u s -R s ı s ,
The magnetic Lagrangian L m has the following form that coincides here with the magnetic energy H m :

L m = L m 2 ı s + ı r e npθ 2 + L f r 2 |ı r | 2 + L f s 2 |ı s | 2 . (13)
More generally physically consistent model should be obtained with a Lagrangian of the form

L IM = J 2 θ2 + L m (θ, ı r , ı * r , ı s , ı * s ) ( 14
)
where L m is the magnetic Lagrangian expressed with the rotor angle and currents. It is 2π np periodic versus θ. Any physically admissible model of a threephase induction machine reads

d dt J θ = ∂L m ∂θ -τ L , 2 d dt ∂L m ∂ı * r = -R r ı r , 2 d dt ∂L m ∂ı * s = u s -R s ı s , (15)
where the rotor and stator flux are given by

φ r = 2 ∂L m ∂ı * r , φ s = 2 ∂L m ∂ı * s .
In general the magnetic energy does not coincide with L m . It is given by (28) that yields:

H m (θ, ı r , ı * r , ı s , ı * s ) = ∂L m ∂ı r ı r + ∂L m ∂ı * r ı * r + ∂L m ∂ı s ı s + ∂L m ∂ı * s ı * s -L m . (16)

Saturation models

A simple way to include saturation effects is to consider that the main inductances L m appearing in (13) depends on the modulus ρ = ı s + ı r e npθ . Thus we consider the following magnetic-saturation Lagrangian:

L m = L m ı s + ı r e npθ 2 ı s + ı r e npθ ı * s + ı * r e -npθ + L f r 2 ı r ı * r + L f s 2 ı s ı * s . (17) Since (L ′ m = dLm dρ ) ∂Lm ∂θ = n p ℑ(ı * r e -npθ ıs) |ıs+ıre np θ | L ′ m and ∂L m ∂ı * r = ı s e -npθ + ı r 2 |ı s + ı r e npθ | L ′ m , ∂L m ∂ı * s = ı s + ı r e npθ 2 |ı s + ı r e npθ | L ′ m ,
the saturation model (formula (15) with L m given by ( 17)) reads

             d dt J θ = n p ℑ Λ m ı * r e -npθ ı s -τ L d dt Λ m ı r + ı s e -npθ + L f r ı r = -R r ı r d dt Λ m ı s + ı r e npθ + L f s ı s = u s -R s ı s (18) with Λ m = L m + |ıs+ıre np θ | 2
L ′ m function of ı s + ı r e npθ . We recover here usual saturation models (see, e.g., [START_REF] Boldea | The Induction Machine Handbook[END_REF], page 428). Notice the similarity with permanent-magnet machines and [START_REF] Yourgrau | Variational Principles in Dynamics and Quantum Theory[END_REF]. Following (16), the associated magnetic energy reads then

H m = L m + ı s + ı r e npθ L ′ m 2 ı s + ı r e npθ 2 + L f r 2 |ı r | 2 + L f s 2 |ı s | 2 .
Notice also that such similarity is not complete and could be misleading: contrarily to (18) that can be derived intuitively form (12) by replacing L m by Λ m , the above magnetic energy is not provided by ( 13) with L m replaced by Λ m .

Space-harmonics with saturation

We can take into account space-harmonic effects by adding their contribution to the magnetic Lagrangian (17). According to [START_REF] Boldea | The Induction Machine Handbook[END_REF], page 298, the iron path in general gets shorter as the harmonic order gets higher. Thus saturation effect has relatively smaller influence on the spatial harmonics. Following [START_REF] Fudeh | Modeling and analysis of induction machines containing space harmonics. Part-I: modeling and transformation[END_REF], the Lagrangian L ν of harmonic ν is

L ν = L ν 2 ı s ı * r e -σν νnpθ + ı * s ı r e σν νnpθ ( 19
)
with L ν a small parameter (|L ν | ≪ L m ) and with σ ν = ±1 depending on arithmetic conditions on ν (see [START_REF] Fudeh | Modeling and analysis of induction machines containing space harmonics. Part-I: modeling and transformation[END_REF], equations ( 25) to ( 29)). The total magnetic Lagrangian now reads

L m = L m ı s + ı r e npθ 2 ı s + ı r e npθ ı * s + ı * r e -npθ + L f r 2 ı r ı * r + L f s 2 ı s ı * s + L ν 2 ı s ı * r e -σν νnpθ + ı * s ı r e σν νnpθ . ( 20
)
Now the saturation model ( 18) is changed as follows:

d dt " J θ" = npℑ "" Λme -npθ + Lν σν νe -σν νnp θ " ı * r ıs " -τL d dt " Λm " ır + ıse -npθ " + L f r ır + Lνıse -σν νnp θ " = -Rrır d dt " Λm " ıs + ıre npθ " + L f s ıs + Lν ıre σν νnp θ " = us -Rsıs (21)
Following ( 16), the associated magnetic energy reads then

H m = L m + ı s + ı r e npθ L ′ m 2 ı s + ı r e npθ 2 + L f r 2 |ı r | 2 + L f s 2 |ı s | 2 + L ν 2 ı s ı * r e -σν νnpθ + ı * s ı r e σν νnpθ .
Several space-harmonics can be included in a similar way. Moreover saturation of space-harmonics can be also tackled just by choosing L ν as a function of |ı s + ı r e -npθ |. As far as we know such explicit models including magnetic saturation and space harmonics have never been given.

Observability issues at zero stator frequency

The sensorless control case is characterized by a load torque τ L constant but unknown, control inputs u s and measured outputs ı s . Models derived from (4) for permanent-magnet machines (resp. from (15) for inductions machines) can be always written in state-space form

d dt X = f (X, U ), Y = h(X) (22)
where

X = (τ L , θ, θ, ℜ(ı s ), ℑ(ı s )) (resp. X = (τ L , θ, θ, ℜ(ı r ), ℑ(ı r ), ℜ(ı s ), ℑ(ı s ))) with U = (ℜ(u s ), ℑ(u s )), Y = (ℜ(ı s ), ℑ(ı s )) and d dt τ L = 0.
A stationary regime at zero stator frequency corresponds then to a steady state ( X, Ū , Ȳ ) of (22) satisfying f ( X, Ū ) = 0 and Ȳ = h( X). The tangent linear system around this steady state is then

d dt x = Ax + Bu, y = Cx (23)
where A = ∂f ∂X ( X, Ū ), B = ∂f ∂U ( X, Ū ) and C = ∂h ∂X ( X). If we assume that the linearized system (23) is observable, the Kalman criteria implies that the rank of the matrix C A must be equal to dim(X). If it is the case, the mapping X → (f (X, Ū ), h(X)) is maximum rank around X. This maximum rank condition just means that the set of algebraic equations characterizing the steady-state from the knowledge of Ū and Ȳ , f (X, Ū) = 0 and h(X) = Ȳ admits around X the maximum rank dim(X). Such rank is not changed by any invertible manipulations of this set of equations characterizing the steady-state from the knowledge of the input and output values, Ū and Ȳ . Putting the implicit Euler-Lagrange equations ( 4) and (15) into their explicit state-space forms (22) involves such invertible manipulations.

For permanent-magnet machines described by ( 4), this set of equations yields the following mapping

(τ L , θ, θ, ı s ) → (0, θ, ∂L m ∂θ -τ L , Ū -R s ı s , ı s )
where L m depends on θ and ı s . Its rank should be maximum, i.e., equal to 5. This is not the case since its rank is obviously equal to 4. For induction machines described by (15), the mapping is

(τ L , θ, θ, ı r , ı s ) → (0, θ, ∂L m ∂θ -τ L , -R r ı r , Ū -R s ı s , ı s )
where L m depends on θ, ı s and ı r . Its rank is equal to 6 whereas the maximum rank is 7. The above arguments yield following proposition:

Proposition 1. Any dynamical model of permanent-magnet machines (4) (resp. induction machines (15)) is unobservable around zero stator frequency regime when the measured output is the stator current ı s and the load torque is constant but unknown. By unobservable we mean that:

• to any constant input and output ūs and īs satisfying ūs = R s īs correspond a one dimensional family of steady states parameterized by the scalar variable ξ with τ L = ∂Lm ∂θ (ξ, īs , ī * s ), θ = ξ, ı s = īs for the permanent-magnet machines, -τ L = ∂Lm ∂θ (ξ, 0, 0, īs , ī * s ), θ = ξ, ı r = 0, ı s = īs for induction machines.

• the linear tangent systems around such steady-states are not observable;

Conclusion

The models proposed in this note, [START_REF] Yourgrau | Variational Principles in Dynamics and Quantum Theory[END_REF] for permanent-magnet machine and (21) for induction drives, are based on variational principles and Lagrangian formulation of the dynamics. Such formulations are particularly efficient to preserve the physical insight while maintaining a synthetic view without describing all the technological and material details (see [START_REF] Yourgrau | Variational Principles in Dynamics and Quantum Theory[END_REF] for an excellent and tutorial overview of variational principles in physics). Extensions to network of machines and generators connected via long lines can also be developed with similar variational principles and Euler-Lagrange equations with complex currents and voltages.

In this note we have put the emphasis on currents and thus Lagrangian modelling. Since flux variables are conjugated to current variables, Hamiltonian modelling is also possible when fluxes are used instead of currents. For permanent-magnet machines, the Hamiltonian counterpart of Lagrangian models (4) reads:

d dt J θ = - ∂H m ∂θ -τ L , d dt φ s = u s -2R s ∂H m ∂φ * s
where the magnetic energy H m is considered as a function of the rotor angle θ, the stator flux φ s and its complex conjugate φ * s . The stator current ı s corresponds then to 2 ∂Hm ∂φ * s . For induction machines, the Hamiltonian formulation associated to (15) becomes

d dt J θ = - ∂H m ∂θ -τ L , d dt φ r = -2R r ∂H m ∂φ * r , d dt φ s = u s -2R s ∂H m ∂φ * s .
The magnetic energy H m now depends on θ, the rotor flux φ r and its complex conjugate φ * r , the stator flux φ s and its complex conjugate φ * s . The rotor (resp. stator) current is then given by 2 ∂Hm ∂φ * r (resp. 2 ∂Hm ∂φ * s ). As for Lagrangian modelling, one can modify the magnetic energies of the standard models (1) and (12) to include, for example, magnetic-saturation or space-harmonics effects. This yields new formula expressing H m as function of angle and fluxes. The corresponding flux-based models are then given by the above equations.

[Lagrangian and Hamiltonian with complex variables] It is explained in [4, page 87] how to use complex coordinates for Lagrangian and Hamiltonian systems. Here we propose a straightforward extension where the complexification procedure is only partial. Such extension cannot be found directly in text-book. In the context of electrical drives, such complexification applies only on electrical quantities whereas mechanical ones remain untouched.

Assume we have a Lagrangian system with generalized positions q ∈ R n , n ≥ 3 and an analytic Lagrangian L(q, q). Let us decompose q into two set of components:

• the first set q c = (q 1 , . . . , q 2n c ) with 0 < 2n c ≤ n will be identified with n c complex numbers q c k = q 2k-1 + q 2k , k = 1, . . . , n c ; • the second set q r = (q 2n c +1 , . . . , q n ) gathers the n r = n -2n c components that will remain untouched and real.

Thus we can identified q with (q c , q r ) where q c ∈ C n c and q r ∈ R n r . Since the Lagrangian L is a real-value and analytic function, it can be seen as an analytic function of the complex variables q c , qc , q c * , qc * and of the real variables q r and qr (q c * corresponds to the complex conjugate of q c ). This function with be denoted by L(q c , q c * , q r , qc , qc * , qr ) and is equal to L(q, q) where q

= " q c 1 + q c * 1 2 , q c 1 -q c * 1 2 , . . . , q c n c + q c * n c 2 , q c n c -q c * n c 2 , q r 1 , . . . , q r n r « q = " qc 1 + qc * 1 2 , qc 1 -qc * 1 2 , . . . , qc n c + qc * n c 2 , qc n c -qc * n c 2 , qr 1 , . . . , qr n r « .
Let us consider the Euler-Lagrange equations

d dt ∂L ∂ qk = ∂L ∂q k + S k , k = 1, . . . , n (24)
where the S k -terms correspond to non conservative energy exchanges with the environment. Similarly to q, we decompose S = (S 1 , . . . , S n ) into S c ∈ C n c and S r ∈ R n r : S is identified with (S c , S r ). We will reformulate these equations with L and its partial derivatives.

k ′ = ∂L ∂ q2k ′ -1 - ∂L ∂ q2k ′ , 2 ∂ L ∂ qc * k ′ = ∂L ∂ q2k ′ -1 +  ∂L ∂ q2k ′ .
This provides the following complex formulation of the real Euler-Lagrange equations (24)

d dt 2 ∂ L ∂ qc * k = 2 ∂ L ∂q c * k + S c k , k = 1, . . . , n c (25)
d dt ∂ L ∂ qr k = ∂ L ∂q r k + S r k , k = 1, . . . , n r . ( 26
)
In the usual complexification procedure ([4, page 87]) the coefficient 2 appearing in the above equations does not appear. This is due to our special choice q c k = q 2k-1 + q 2k instead of the usual choice q c k = q 2k-1 +q 2k √ 2

. This special choice preserves the correspondence, commonly used in electrical engineering, between complex and real electrical quantities Let us assume that, for each q, the mapping q → ∂L ∂ q is a smooth bijection. Then the Hamiltonian formulation of (24) reads

d dt q k = ∂H ∂p k , d dt p k = - ∂H ∂q k + S k , k = 1, . . . , n (27)
with H = ∂L ∂ q • q -L and p = ∂L ∂ q . Let us decompose p into p c ∈ C n c and p r ∈ R

H = ∂ L ∂ qc qc + ∂ L ∂ qc * qc * + ∂ L ∂ qr qr - L (28)
where H denotes the Hamiltonian H when is a considered as a function of (q c , q c * , q r , p c , p c * , p r ). Then (27) becomes

d dt q c k = 2

  ) with such L m . Denote λ ′ = dλ dρ . With ∂λ ∂θ = n p ℑ(īe -np θ ıs) |ıs+īe np θ | λ ′ and ∂λ ∂ı * s = ıs+īe np θ 2|ıs+īe np θ | λ ′ , we get the following dynamical model with both saliency and magnetic-saturation effects:

  For k = 2n c + k ′ with k ′ = 1, . . . , n r , they remain unchanged k ′ . The two scalar equations corresponding to k = 2k ′ -1 and k = 2k ′ with k ′ = 1, . . . , n c yield to the following single complex equation d dt 2

  n r . Then p r = ∂ e L ∂ qr and p c = 2 ∂ e L ∂ qc * . Simple computations yield another derivation of the Hamiltonian from L directly:

  k , k = 1, . . . , n r


K Abdennadher

P Venet

G Rojat

J M Rétif

C Rosset

Karim

Abdennadher Karim
email: karim.abdennadher@mgeups.com

Online Monitoring Method and Electrical Parameter Ageing Laws of Aluminium Electrolytic Capacitors Used in UPS

Keywords: <<Industrial application>>, <<Uninterruptible Power Supplies (UPS)>>, <<Estimation technique>>, <<Maintenance>>, <<Real time processing>>

This paper presents a work used for a low cost method to create a real time condition monitoring and a predictive maintenance system of an electrolytic capacitor in Uninterrupted Power Supplies (UPS). It consists of evaluating the evolution laws for equivalent series resistance (ESR) and the capacitance C values of the electrolytic capacitors versus operating time. Experimental results of parameter evolutions and the law estimations are presented. The laws funded are used in UPS where waveforms are continuously varying in amplitude, frequency and versus temperature. The estimated laws are used for the on line failure prediction method which has the merits of using only the existent resources in UPS.

Introduction

Electrolytic capacitors have been widely used in power electronic systems because they can achieve high capacitance and voltage ratings with volumetric efficiency and low cost. This type of capacitors have been traditionally used for filtering, coupling, timing networks, by-pass and other many applications in power electronics requiring a cost effective and volumetric efficiency components. We know that the common faults in electrolytic capacitor include initial catastrophic failures due to manufacturing or misapplication defects, and wear-out faults which can't be avoid. Unfortunately, electrolytic capacitors are some of the weakest components in power electronic converters [START_REF] Abdennadher | A real time predictive maintenance system of electrolytic capacitors used in uninterruptible power supplies[END_REF][START_REF] Venet | Dectection of faults of filter capacitors in a converter. Application to predictive maintenance[END_REF][START_REF] Venet | Realization of a smart electrolytic capacitor circuit[END_REF][START_REF] Lahyani | Failure prediction of electrolytic capacitors during operation of a switchmode power supply[END_REF][START_REF] Harada | Use of ESR for deterioration diagnosis of electrolytic capacitor[END_REF][START_REF] Imam | Real-Time Condition Monitoring of the Electrolytic Capacitors for Power Electronics Applications[END_REF][START_REF] Gasperi | Life prediction modeling of bus capacitors in ac variable-frequency drives[END_REF][START_REF] Sankaran | Electrolytic capacitor life testing and prediction[END_REF][START_REF] Ondel | A decision system for electrolytic capacitors diagnosis[END_REF][START_REF] Lee | Online Capacitance Estimation of DC-link Electrolytic Capacitors for Three-phase AC/DC/AC PWM Converters using Recursive Least Squares Method[END_REF][START_REF] Myungchul | Condition Monitoring of DC Link Electrolytic Capacitors in Adjustable Speed Drives[END_REF][START_REF] Perisse | Simple Model of Electrolytic Capacitor Taking into Account the Temperature and Aging Time[END_REF]. Many papers have proposed different methods or algorithms to determinate the ESR and/or capacitance C of the electrolytic capacitor [START_REF] Venet | Dectection of faults of filter capacitors in a converter. Application to predictive maintenance[END_REF][START_REF] Venet | Realization of a smart electrolytic capacitor circuit[END_REF][START_REF] Lahyani | Failure prediction of electrolytic capacitors during operation of a switchmode power supply[END_REF][START_REF] Harada | Use of ESR for deterioration diagnosis of electrolytic capacitor[END_REF][START_REF] Imam | Real-Time Condition Monitoring of the Electrolytic Capacitors for Power Electronics Applications[END_REF][START_REF] Gasperi | Life prediction modeling of bus capacitors in ac variable-frequency drives[END_REF][START_REF] Sankaran | Electrolytic capacitor life testing and prediction[END_REF][START_REF] Ondel | A decision system for electrolytic capacitors diagnosis[END_REF][START_REF] Lee | Online Capacitance Estimation of DC-link Electrolytic Capacitors for Three-phase AC/DC/AC PWM Converters using Recursive Least Squares Method[END_REF][START_REF] Myungchul | Condition Monitoring of DC Link Electrolytic Capacitors in Adjustable Speed Drives[END_REF][START_REF] Perisse | Simple Model of Electrolytic Capacitor Taking into Account the Temperature and Aging Time[END_REF]. However, many parameters such as additional measurements and many computations are required, which makes it complicated, difficult, expensive and impractical for actual application. Nowadays, a UPS disposes of many accurate sensors, numerical treatment systems and powerful computation resources which are used to control and regulate the UPS in order to improve its performance and efficiency. The proposed method has the merits of making a real time predictive maintenance system of electrolytic capacitors with using existing resources in the UPS. This predictive maintenance system works in background task and without perturbing operating system. Also this document presents studying the ageing of electrolytic capacitors to determine the evolution laws of the internal parameters: equivalent series resistance (ESR) and capacitance (C). These evolution laws will be used for the ageing algorithm that will be integrated on DSP as a part study of the real time predictive maintenance system of aluminium electrolytic capacitors used in UPS [START_REF] Abdennadher | A real time predictive maintenance system of electrolytic capacitors used in uninterruptible power supplies[END_REF]. This part of work is divided into two main steps to determine all the evolution laws versus temperature, frequency and operation time.

Electrolytic capacitors monitoring

Estimation of ESR and C

The electrolytic capacitors can be modelled as a series combination of a capacitance, inductance, and resistance as showing in Fig. 1 As it can be observed from bode-plot in Fig. 2, there exists three widely separated frequency bands. The resonant frequency of the capacitor is given by this equation.

f r = 1 2.π. √ ESL.C (1)
Therefore, generally converters in UPS works at low frequency band compared to the resonant frequency's one, so ESL is usually neglected and the equivalent model of the capacitor is RC series equivalent circuit.

The transfer function of this model is given by:

H(p) = Uc(s) Ic(s) = ESR.C.s + 1 C.s (2)
However, converters are controlled with using digital treatments; we consider the z-transform corresponding of Equation (2) to represent the discrete time domain with using the Tustin bilinear method which is given by the following equation [START_REF] Venet | Realization of a smart electrolytic capacitor circuit[END_REF].

H z -1 = b 0 + b 1 .z -1 1 -z -1 = ESR + T s 2.C + T s 2.C -ESR .z -1 1 -z -1 (3)
Ts is the sampling period. The bode-plot in Fig. 2 shows that the continuous and discrete time domains are equivalent if we use the bilinear method discrete time. Also, we can see the temperature change effect on the different parameters of the capacitor.

To identify b0 and b1 (cf. equation ( 3)), so the two parameters ESR and C, we use a least square algorithm. Compared to [START_REF] Abdennadher | A real time predictive maintenance system of electrolytic capacitors used in uninterruptible power supplies[END_REF], this paper presents the implementation of the predictive maintenance system on an AC / DC converter. Also we discuss on this paper how to improve the precision identification of capacitor parameters (ESR and C). In fact, the identification of these two parameters is made separately by using simple frequency adjusting channel measurements needed to identify each of these two parameters.

A Matlab simulation was performed using dynamic equation from equivalent circuit shown in Fig. 3 of an AC/DC converter with 230V AC voltage input, for a switching frequency of 5 kHz and which is regulated to have 400V DC voltage output. A Recursive Least Square algorithm which is a particular case of Kalman filters [START_REF] Foulard | Commande et régulation par calculateur numérique[END_REF] was used. In practice, we don't have a capacitor current sensor; because it's difficult and expansive to implement it in UPS. To perform and regulate the output voltage, we have only one current sensor used to measure the input inductance L current of the boost converter. However, we can identify the capacitor current [START_REF] Abdennadher | A real time predictive maintenance system of electrolytic capacitors used in uninterruptible power supplies[END_REF] with using Kirchhoff laws referencing to the schematic of a boost converter given in Fig. 3 and the PWM generated to control the converter.

i c = i L .PW M -avg i L .PW M (4)
Simulation results presented for the parameters identification are summarized in Table I.

The %Error presented on this table are the ones calculated between ESR and C values fixed on simulation and the estimated values given by the least square algorithm.

Experimental study for online monitoring method Experimental setup

To verify the validity of the proposed method, an experimental study was performed on an AC/DC converter with input AC voltage of 230 V, 2200 µF -450 V filtering capacitor, an R-load of 100 Ω, a switching frequency of 5 kHz. Then a fixed point DSP controller (TMS320F2812) was used for i L , u c , i c and capacitor ripple voltage ΔUc acquisition measurements and switching operations required to produce a 400 V DC output voltage. Moreover, a recursive least square algorithm (a particular case of Kalman filter) had been implemented to work in background task and to identify the two parameters ESR and C. As existing literature shows [START_REF] Landau | Digital Control Systems Design, Identification and Implementation[END_REF], we know that using different cases of Kalman filter provides centred input and output system measurements which are respectively capacitor current and capacitor voltage. Then, the capacitor voltage and the capacitor current are passed through a high-pass filter with cut-off frequency of around f c = 200 Hz to suppress DC component of voltage measurements and without including phase shift. Consequently, we have only the capacitor ripple voltage which is in phase with capacitor current.

An FFT analysis respectively on the capacitor ripple voltage and capacitor current is made and we conclude that the bandwidth of the system is about 40 kHz, therefore and to respect Shannon theorem we choose a sampling frequency f s = 80 kHz.

To eliminate signal distortion anti-aliasing filters with a cut-off frequency of f c = f s 2 ≈ 40kHz are used like showing in the global capacitor condition monitoring scheme Fig. 4.

Experimental results

Experimental results for online estimated parameters are summarized in Table II. The least square algorithm implemented on DSP to work in background task without perturbing operat- ing system use about 1K x 16bits of DSP memory and take about two seconds to converge and gives the ESR and C estimated values. This time is too small compared to the life time of a capacitor. It can be seen that there is a small and less error between estimated and actual values of ESR and C compared to experimental results presented in [START_REF] Abdennadher | A real time predictive maintenance system of electrolytic capacitors used in uninterruptible power supplies[END_REF] due to the fact that we identify the two parameters separately. These results are consistent with the measured value using an impedancemeter.

Ageing algorithm

Ageing laws

To deduce the number of hours before failure with knowing the ambient temperature measurement and the ESR and C values found in real time with the least square algorithm, we implement an ageing algorithm (Fig. 4). This algorithm uses the limit values to consider that the capacitor is at the end of life and which are fixed at:

• An increase of 100% of the ESR value compared to the initial one: ESR limit =2.ESR(0)

• A decrease of 20% of the capacitance value compared to the initial one: C limit = 0.8.C(0)

The parameters ESR and C depend on the inverter switching frequency and the ambient temperature.

The experimental tests carried out off-line and presented latter allow us to highlight the following laws: Parameters evolutions versus temperature T' by knowing the ambient temperature: these laws allow us to determine the values of ESR (0) and C (0) using the following laws [START_REF] Perisse | Simple Model of Electrolytic Capacitor Taking into Account the Temperature and Aging Time[END_REF]:

ESR(0) = α + β. exp -T γ (5)
C(0) = χ + λ.T (6)
with α, β, χ and λ coefficients depend of type of capacitor.

The evolution curve parameters of experimental results allow us to determine the parameter evolution laws as a function of operating time (t op ) [START_REF] Perisse | Simple Model of Electrolytic Capacitor Taking into Account the Temperature and Aging Time[END_REF].

ESR (t op ESR ) = (ESR (0) + A 1 ) . exp (B 1 .t op ESR ) (7)
C (t op C ) = (E.C(0) + F.t op C ) (8)
with A 1 , B , E and F coefficients depend of type of capacitor. From these evolution laws, we determine the ageing time (ta) corresponding to the two ageing temperatures (85 • C and 90 • C) using the following equations:

ta ESR = log ESR limit ESR(0)+A 1 B 1 ( 9
)
ta C = C limit -E.C (0) F ( 10
)
After having calculated the ageing time for two temperatures (85 • C and 90 • C for example), we calculate the activation energies (Ea ESR and Ea C ) used in the Arrhenius law as follows: [START_REF] Perisse | Simple Model of Electrolytic Capacitor Taking into Account the Temperature and Aging Time[END_REF] ta ESR : Ageing time corresponding to the ESR limit at T a . ta C : Ageing time corresponding to the C limit at T a . k = Boltzmann constant (8,617.10 -5 eV/ • K). Ea ESR : Activation energy with taking into account ESR as an indicator of ageing. Ea C : Activation energy taking into account C as an indicator of ageing. The procedure to determinate different laws and different parameters off-line is illustrated in Fig. 5 below: In real time and online, we have three values which are the equivalent series resistance (ESR), the capacitance (C) emerged by the identification algorithm and the ambient operating temperature T (Fig. 4). The value of temperature T allows us also to determine the time needed for the ageing with the use of following equations:13

Ea ESR = k. log ta ESR(85 • C) ta ESR(90 • C) (90-85) (273+90).(273+85) (11) Ea ESR = k. log ta ESR(85 • C) ta ESR(90 • C) (90-85) (273+90).(273+85)
t ESR = ta ESR . exp Ea ESR k . T a -T (T a + 273) . (T + 273) ( 13
)
t C = ta C . exp Ea C k . T a -T (T a + 273) . (T + 273) (14)
T a : Ageing temperature (85 • C for example). T : Ambient temperature (25 • C for example). The ESR and C values identified in real time and also the ESR (0) and C (0) will be used in the equations ( 9) and ( 10) to determine the respective operating time as shown in following equations:

t op ESR = log( ESR ESR(0)+A 1 ) B 1 ( 15
)
t op C = C -E.C (0) F (16)
The time until failure t f ailure (ESR) with the use of equations ( 13) and ( 14) is given by:

t f ailure (ESR) = ⎛ ⎝ ta ESR . exp Ea ESR k . T a -T (T a + 273) . (T + 273) - log ESR ESR(0)+A 1 B 1 ⎞ ⎠ (17)
The figure below explains the different times given by equations ( 13), ( 15) and ( 17 14) and ( 16) is given by:

t f ailure (C) = ta C . exp Ea C k . T a -T (T a + 273) . (T + 273) - C -E. (0) F ( 18
)
The capacitor time until failure is the less time between the two ones:

t f ailure = Min(t f ailure (ESR),t f ailure (C)) (19)
The

Experimental results used for ageing algorithm

We have for these work different types of capacitors according to four different manufacturers. For each one, we have 4 to 5 components indicated by the letters A, B ... Q and for each series we have 3 components of the same electrical characteristics with a capacitance tolerance of about ±20%. Each component has an index of 1, 2 or 3 (e.g. A1, A2 and A3) since they have the same nominal capacitance and voltage. All capacitors are first characterized versus temperature to determine the parameter evolution laws versus temperature (cf. figure 8 as example). The parameters of the equation ( 4) are determined by the non-linear least squares method of Levenberg Marquardt which allows a very close approximation of the model. For the example presented in Fig. 8, the method gives the following values:

Conclusion

Due to the large capacity and low cost, electrolytic capacitors capable of energy storage and voltage regulation are used for almost all types of power electronics system. Electrolytic capacitor, which is usually affected by wear-out faults, plays a very important role for the power-electronic systems quality and reliability. Therefore it is important to monitor the condition of an electrolytic capacitor in real-time to predict the failure. A new method is proposed to detect in real time the changes in ESR and capacitance C value of an electrolytic capacitor in order to realize a real-time failure prediction of an electrolytic capacitor. For the proposed method, capacitor current and capacitor ripple voltage measurements using cheap and simple analog circuits are required. Simulation results and hardware experimental implementation show that the proposed electrolytic capacitor failure-prediction technique can be applied to power-electronic systems successfully. Also, in continuation of our research and to improve the parameter estimation, we decide to implement adjustable gains used for the capacitor current and capacitor voltage measurements and which allow us to minimize the dynamic numerical errors. In this paper, we have also presented the experimental results achieved on different electrolytic capacitors to find the different evolution laws of electrical parameters ESR and C as a function of temperature. Then we discussed the results of the ageing of these standard components to 90 • C after about 5 months of ageing. These results enabled us to compare the capacitors from different manufacturers and conclude on their quality. As perspective to our work, we will finish the ageing standard 85 • C to determine all the laws of ageing capacitors over time. The final objective is to determine these laws directly from well detailed data sheets provided by manufacturers without having to make these experimental tests used for each capacitor. These laws of evolution are used in the ageing algorithm implemented in real time to determine the remaining time before failure of electrolytic capacitors used in UPS.

  First, the capacitors are characterized at different temperatures (-20 • C, 0 • C, 10 • C, 25 • C, 40 • C, 55 • C and 85 • C). This allows us to determine the evolution laws of ESR and C versus temperature. After, two standards ageing are launched at two different temperatures (85 • C and 90 • C). This type of ageing consists of supplying the components under nominal voltage (normally used on UPS) and the maximum operating temperature allowed with the use of a climatic chamber. This allows us to analyze the capacitors internal parameter behaviors versus operating time and determine the ageing laws of parameters ESR and C of each component. This paper presents the experimental results of temperature characterization and standard ageing at 90C after 3390 hours (about 5 months) of ageing for different types and manufacturers of electrolytic capacitors.

Figure 1 :

 1 Figure 1: Electrical equivalent circuit of an aluminium electrolytic capacitor

Figure 2 :

 2 Figure 2: Bode-plot of electrolytic capacitor in continuous and discrete time domain

Figure 3 :

 3 Figure 3: Schematic of an AC/DC converter circuit with a DC capacitor circuit

Figure 4 :

 4 Figure 4: Real time electrolytic capacitor condition monitoring in UPS

Figure 5 :

 5 Figure 5: Off-line procedure to determinate the parameter evolution laws

  T ESR : Extrapolated time at with ESR as an ageing indicator. T C : Extrapolated time at with C as an ageing indicator.

Figure 6 :

 6 Figure 6: ESR evolution versus operating time

Figure 7 :

 7 Figure 7: Ageing algorithm

Table I

 I

		: Simulation results for fixed parameters
		fixed	estimated % error
	ESR 32 mΩ 32.6 mΩ	0.2
	C	2200 µF 2184 µF	0.7

Table II :

 II Experimental results for estimated parameters

		actual	estimated[min, max] % error
	ESR 31.2 mΩ [29.6 mΩ,32.5 mΩ] +/ -5
	C	2200 µF	[2289µF, 2585µF]	+/ -6

EPE 2009 -Barcelona P.3

The parameters of the equation [START_REF] Imam | Real-Time Condition Monitoring of the Electrolytic Capacitors for Power Electronics Applications[END_REF] are determined in the same way as those identified in ESR. For the example presented in Fig. 8, we have the following values: χ=0.000812F λ=4,600790.10 -0.007 F/ • C The index 1 capacitors are used for the first ageing at 85 • C and those of index 3 are used for the ageing at 90 • C. Some experimental results show the parameter evolutions of ESR and C after 3390 hours of ageing (about 5 months) for capacitors with index 3 (ageing at 90 • C). Also we use a Matlab algorithm to make a linear interpolation and extrapolation with the use of experimental data. This allows us to find the parameters evolution laws and also the time until failure corresponding of ESR and the capacitance C limit like presented in Table III. The ESR and C variations versus temperature for a capacitor 1000µF / 450 V are shown in Fig. 8.