Sophie Jacquin

Lucien Mousin

Igor Machado

El-Ghazali Talbi

Laetitia Jourdan

A Comparison of Decoding Strategies for the 0/1 Multi-objective Unit Commitment Problem

Keywords: UCP, metaheuristics, heuristic, multi-objective optimization

In the single objective Unit Commitment Problem (UCP) the problem is usually separated in two sub-problems : the commitment problem which aims to fix the on/off scheduling of each unit and the dispatching problem which goal is to schedule the production of each turned on unit. The dispatching problem is a continuous convex problem that can easily be solved exactly. For the first sub-problem genetic algorithms (GA) are often applied and usually handle binary vectors representing the solutions of the commitment problem.Then the solutions are decoded in solving the dispatching problem with an exact method to obtain the precise production of each unit. In this paper a multiobjective version of the UCP taking the emission of gas into account is presented. In this multi objective UCP the dispatching problem remains easy to solve whereas considering it separatly remains interesting. A multi-objective GA handling binary vectors is applied. However for a binary representation there is a set of solutions of the dispatching problem that are pareto equivalent. Three decoding strategies are proposed and compared. The main contribution of this paper is the third decoding strategy which attaches an approximation of the Pareto front from the associated dispatching problem to each genotypic solution. It is shown that this decoding strategy leads to better results in comparison to the other ones.

Introduction

The UCP is used to find the scheduling of commissioning and production of generating units that minimizes the production cost. However, environmental protection has become a major issue. In response to the growth of the negative impacts on the environment and due to the growing importance of environmental interest in society, governments have developed and implemented laws or technical standards in order to reduce the negative impacts of human activity on the environment. For this reason new modeling of the UCP taking into account some limitation constraints on gas emissions [START_REF] Catalao | A practical approach for profitbased unit commitment with emission limitations[END_REF] has been proposed. A bi-objective model has also been proposed, but is usually solved by reducing the two objectives to one, using the weighted sum approach [START_REF] Zhang | Multi-objective unit commitment fuzzy modeling and optimization for energy-saving and emission reduction [j[END_REF][START_REF] De Moura Gomes | Metaheuristics for the Unit Commitment Problem The Constraint Oriented Neighbourhoods Search Strategy[END_REF]. In [START_REF] Srinivasan | An evolutionary algorithm for evaluation of emission compliance options in view of the clean air act amendments[END_REF], a multiobjective version of the UCP is solved with a genetic algorithm, but it does not directly exploit the concept of Pareto dominance. The authors applied a classical genetic algorithm, with the exception that the selection process is based on a specific version of the tournament selection. Two individuals are randomly selected from the population and a stochastic competition of the objective that are chosen randomly, is performed to determine the winner that will survive into the next generation.

The UCP can be separated into two sub-problems. The commitment problem which is to give the on/off scheduling of each unit and then the dispatching problem which is to give the exact production for each turned on unit. In the single objective case, the dispatching problem is a quadratic continuous problem which is easy to solve exactly. For this reason evolutionary algorithms generally handle binary vectors giving the on/off scheduling of each unit to solve the single-objective UCP [START_REF] Damousis | A solution to the unitcommitment problem using integer-coded genetic algorithm[END_REF][START_REF] Jeong | A thermal unit commitment approach using an improved quantum evolutionary algorithm[END_REF][START_REF] Kazarlis | A genetic algorithm solution to the unit commitment problem[END_REF][START_REF] Lau | Quantum-inspired evolutionary algorithm approach for unit commitment[END_REF][START_REF] Swarup | Unit commitment solution methodology using genetic algorithm[END_REF]. Then to obtain a complete description of the solution (the phenotypic solution), the production of each turned on unit is determined optimally using a λ-iteration method [START_REF] Saramourtsis | Genetic algorithm solution to the economic dispatch problem-application to the electrical power grid of crete island[END_REF]. This binary representation takes advantage of a real one because the search space is considerably reduced. For the multi-objective the dispatching problem remains interesting to exploit separately because it has still good properties. All the solutions of the dispatching problem are supported and can be found by the scalarizing method. This method is to transform the problem with two objectives f 1 and f 2 into a single objective problem optimizing the function λf 1 + (1λ)f 2 , where λ ∈ [0, 1]. The quality of the results found in the single objective case and the fact that the multiobjective version of the dispatching problem has good properties lead us to a specific interest in considering a similar two level method of resolution for the multi-objective UCP.

In this paper a multi-objective GA based on NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: Nsga-II[END_REF] is proposed. This GA handle binary vectors representing the solutions of the commitment problem. Then the phenotypic solutions corresponding to the production of each unit has to be find by a decoding method. Since the dispatching problem is also a multi-objective problem there are many candidates of phenotypic solutions for one genotypic solution. Therefore three different decoding methods are proposed and compared. The first one is to construct the phenotypic solution by solving the dispatching problem by optimizing the function f 1 +f 2 . The second is to add a real λ into the genotypic representation of a solution. The phenotypic solution is decoded by solving the dispatching problem minimizing λf 1 + (1λ)f 2 . The third is to associate a set of Pareto equivalent solutions of the corresponding multi-objective dispatching problem to each genotypic solution. For the last decoding system, the process of fitness and diversity assignment of NSGA-II has to be adapted.The fact that many phenotypic solutions are attached to a single genotypic solution must be taken into account. The main contribution of this paper is to show that the multi decoding embedded approach has an advantage over the two other less complex decoding systems that are proposed.

The paper is organized as follows. In the next section the multi-objective UCP is described in details. Then the three solving methods corresponding to the three binary genetic algorithms using the different decoding strategies are presented. Finally the experimental process and the obtained results are presented and discussed before the conclusion section.

Multi-objective Unit Commitment Problem -MO-UCP

In this section, the MO-UCP is presented in details. This problem is the same as the classical UCP but an objective is added to take into account the gas emission of SO 2 and CO 2 .

Unit Commitment Problem

Unit Commitment Problem (UCP) [START_REF] Kazarlis | A genetic algorithm solution to the unit commitment problem[END_REF] is to schedule generating units online or offline over a scheduling horizon. The goal is to minimize the power production cost while satisfying a set of operational constraints. The production cost includes the fuel and start-up costs. Constraints are capacity of production of each units, minimum up/down time and spinning reserve. UCP is usualy modeled as as a mixed integer non-linear problem. It consists of binary variables u i,t that takes value 1 if a unit i is turned on at time t and 0 otherwise, continuous variables p i,t that denotes their prodduction amounts. It is a very complex problem to solve because of its enormous dimension, a non-linear objective function, and time-dependent constraints. Indeed, the UCP is a NP-complete problem [START_REF] Guan | Optimization based methods for unit commitment: Lagrangian relaxation versus general mixed integer programming[END_REF].

Multi-objective UCP

The first objective is the same as in single-objective UCP. It minimizes the cost of production. This production cost is divided into two components, the fuel cost and the start up cost. For a system of N units and a time horizon of T periods, the objective function can be described as follows:

f 1 (u, p) = T t=1 N i=1 F C i (p i,t ) × u i,t + CS i (T of f i,t-1 ) × (1 -u i,t-1 )u i,t
, where :

-F C i is the fuel cost function of the unit i, which is modeled by a quadratic function:

F C(p i,t ) = a 1,i + a 2,i × p i,t + a 3,i × p 2 i,t
, where a 1,i , a 2,i and a 3,i are real cost coefficients for the unit i.

-CS i is the start-up cost for unit i, which depends on time T of f i,t-1 the unit i has been turned off at time t -1:

CS i (T of f i,t-1 ) = CS cold if T of f i,min + T cs,i ≤ T of f i,t-1 CS hot else ,
where T of f i,min + T cs,i is the time it takes the unit i to become cold. The second objective function measures the SO 2 and CO 2 emissions:

f 2 (p, u) = T t=1 i,ui,t=1 b 0,i + b 1,i p i,t + b 2,i p 2 i,t
The coefficients b 0,i , b 1,i , b 2,i used in this paper are the ones proposed in [START_REF] De Moura Gomes | Metaheuristics for the Unit Commitment Problem The Constraint Oriented Neighbourhoods Search Strategy[END_REF].

The minimization of the objectives of the UCP is subject to the following system and unit constraints:

1. Power balance constraints:

N i=1 p i,t u i,t = D t ∀t ,
where D t is a real number giving the load demand at time t.

Spinning reserve constraints:

N i=1 p i,max u i,t ≥ D t + R t ∀t ,
where R t is a real giving the minimal reserve at time t.

Unit output constraints:

p i,min ≤ p i,t ≤ p i,max ∀t,
where p i,min and p i,max are the lower and upper bounds on the energy production of unit i respectively. 4. Minimum up time limit:

T on i,t-1 ≥ T on i,min × (1 -u i,t )u i,t-1 ∀t,
where T on i,t-1 is the time from which the unit i is turned on at time t -1 and T on i,min is the minimal time during which unit i has to stay turned on. 5. Minimum down time limit:

T of f i,t-1 ≥ T of f i,min × (1 -u i,t-1 )u i,t ∀t,
where T of f i,t-1 is the time from which the unit i is turned off at time t -1 and T of f i,min is the minimal time during which unit i has to stay off.

The feasible outcome vectors of the objective space are compared using the Pareto dominance ≻. In this minimization context, a solution x ∈ Ω is said to dominate a solution y ∈ Ω, denoted by x ≻ y, if they satisfy relation [START_REF] Bleuler | Pisa: A platform and programming language independent interface for search algorithms[END_REF].

∀i ∈ {1, 2}, f i (x) ≤ f i (y) ∃i ∈ {1, 2}, f i (x) < f i (y) (1)
Many evolutionary algorithms involving a decoding system have been proposed to solve the single objective version of the UCP [START_REF] Damousis | A solution to the unitcommitment problem using integer-coded genetic algorithm[END_REF][START_REF] Jeong | A thermal unit commitment approach using an improved quantum evolutionary algorithm[END_REF][START_REF] Kazarlis | A genetic algorithm solution to the unit commitment problem[END_REF][START_REF] Lau | Quantum-inspired evolutionary algorithm approach for unit commitment[END_REF][START_REF] Swarup | Unit commitment solution methodology using genetic algorithm[END_REF]. A binary vector u of size T × N is used to represent the solutions. In this representation each u i,t gives the state of a unit i (on or off) at a given time period t of the scheduling.

Then the exact production of each unit is decoded by solving the dispatching sub-problem (D(u)) using the lambda-iteration method [START_REF] Saramourtsis | Genetic algorithm solution to the economic dispatch problem-application to the electrical power grid of crete island[END_REF]:

min p T t=1 i s.t ui=1 f 1 ((p i,t ) i ) such that: i s.t ui=1 p i,t = D t (2)
p i,min ≤p i,t ≤ p i,max ∀i s.t u i = 1 (3)
The advantage of the binary vector representation over a real vector representation where the productions are directly given is obvious: the search space is considerably reduced. In the multi-objective case, the dispatching problem is also a multi objective one as f 1 and f 2 have to be minimized. But this subproblem has still very good properties that make it easy to solve. It is a convex and continuous bi-objective problem. As the objective functions are convex and the decision variables are defined in a convex set all the solutions are supported [START_REF] Ehrgott | Multicriteria optimization[END_REF] . Then the Pareto front solution is convex and totally defined by the set:

{f 1,λ , f 2,λ |λ ∈ [0, 1]},
where

f 1,λ = f 1 (p * λ ) and f 1,λ = f 1 (p * λ )
, with p * λ solution of the dispatching problem D(u, λ) defined as follows:

p * λ = arg(min p T t=1 i s.t ui=1 λf 1 ((p i,t ) i )) + (1 -λ)f 2 ((p i,t ) i ))
such that (2) and (3) are met.

It seems interesting to consider a two level based method in the multiobjective version which will be similar to those proposed for the single-objective UCP. The method proposed is a GA, NSGA-II, handling binary vectors representing the solutions of commitment problem. Then the production of the units are obtained by using a decoding method. However, since the dispatching problem is a multi-objective one, there are many pareto equivalent solution to fix the production values. It becomes difficult to choose a method to associate phenotypic solutions with genotypic ones. Then three approaches of decoding are proposed and compared. In the first approach the solution associated with a binary representation is the one obtained by solving the dispatching problem for λ fixed to 0.5. This approach is a naive one because it might miss some possibly good solutions. The second approach consists of adding λ to the representation of an individual. Then a binary vector associated with an on/off scheduling can be present many times in the population with different value of λ and all the solutions are reachable. In the last approach, an approximation of the Pareto front of the dispatching sub-problem is associated with each individual. Therefore an adapted version of NSGA-II is proposed to manage the association of many phenotypic solutions to a single genotypic solution.

In the next subsection, the common components of the three methods, which are essentially evolutionary operators, are presented. Then each method are explained in details.

Common components

Each multi-objective GA proposed hereafter are based on NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: Nsga-II[END_REF]. In each case the following operators are used:

Crossovers: Two crossovers are used. The first one is the classical one-point crossover and the second one is an intelligent two-points crossover. It is to randomly choose a window size and if the window size is smaller than the remaining portion of the vector, a new individual is created with the window portion of the worst parent and the remaining portion from the best parent. The reverse is done if the window size is larger than the number of bits in the remaining vector.

Mutations: Two mutations are used. The first one is the standard 1-bit-flipmutation. This operator randomly flips a bit of the vector with a low probability. The second one randomly chooses a window whose size is randomly chosen and flips all the bits of this window.

Repair operator: The aim of this operator is to correct a solution if it does not meet the constraints of demand or of minimum on/off time. Firstly, it corrects the violations of time constraint hour by hour in modifying the states of unities if necessary. Then corrections are done on the constraints of power balance. This is done hour-by-hour. If the maximal capacity of production of the turned on units is lower than the power balance then a unit is turned on. The unit to turn on is chosen randomly among the units that can be turned on depending on the past hour (i.e. if this unit was turned on the previous hour or if it was turned off for a time long enough). The reverse process is done if minimal capacity of the turned on units exceeds the load demand. Naturally this correction process does not guarantee obtaining a feasible solution but this operator speeds up the algorithm and hence increases the possibility of finding such solutions.

Objective function:

The objective values are computed by using the exact production p i,t at each time period t and for each unit i which is a phenotypic solution. The way to obtain phenotypic solutions from genotypic ones is different for each approach. The objective functions correspond to f 1 , the cost production and f 2 , the quantity of SO 2 and CO 2 emission. Nevertheless some penalties have to be added on the violation of time and load constraints. Therefore the objectives are:

       obj 1 = f 1 + c p T t=1 ( N i=1 p i,t -D t ) + c p (T t-1 i,of f (u i,t -x i,t-1 ) + T t-1 i,on (u i,t-1 -u i,t ))
obj 2 = f 2 + c p T t=1 ( N i=1 p i,t -D t ) + c p (T t-1 i,of f (u i,t -u i,t-1 ) + T t-1 i,on (u i,t-1 -u i,t )), (4)
where c p is a constant positive number.

Naive approach

In this approach a genotypic solution is decoded in solving a dispatching problem that is reduced to a single objective one by scalarization. Then it is possible that some pareto optimal phenotypic solutions cannot be reached because they are not solution of the chosen scalarized sub-problem. The decoding process is explained in detail hereafter.

Decoding process: The phenotypic solution (p i,t ) i,t associated to a genotypic is computed by solving the dispatching problem D(u, 0.5) thank to the λ-iteration method.

Scalarized decoding

In this approach an integer (λ u ∈ 0, 100 ) is added to the genotypic representation of a solution. This value is chosen randomly during the initialization process and then can be modified by the evolutionary operators. It is used to define the coefficients of scalarization during the decoding process. Hence, the main difficulty of the previous approach becomes possible to overcome, which is the inaccessibility of some pareto optimal solutions. Actually as the pareto front solution of the dispatching problem will always be convex all the pareto optimal solutions are reatchable by the scalarization method.

Representation: A value λ u is added at the end of the representation vector with a binary encoding. λ u is an integer between 0 and 100. The representation is shown Fig. 1.

Decoding process: The phenotypic solution (p i,t ) i,t associated to a genotypic is computed in solving the dispatching problem D(u, λu 100 ) thank to the λ-iteration method.

{ u 010 • • • 1 { unit 1 010 • • • 0 { unit 2 • • • 110 • • • 0 { unit N 101 • • • 1 { λu Fig. 1.
Representation genotypic of a solution: the scalar λ is included in the representation.

Adaptation on the crossover: The one point crossover is transformed in a two-points crossover with the first crossover point in a locus of the definition of u and the second point in a locus of definition of λ u . It is done in order that the crossover can have a significant impact on λ u .

Adaptation on the mutation: The 1-bit-flip-mutation is applied only on the bits corresponding to u. λ u is mutated by being replaced by a value chosen randomly between 0 and 100 with a normal distribution centered on its original value.

Multi decoding embedded approach

In this approach a genotypic solution is associated with a set of phenotypic solutions. This set of solutions is from the optimal pareto front solution of the dispatching problem associated with the genotypic solution. Fig. 2 helps to un- derstand how genotypic solutions are represented in the objective space. The continuous line is the optimal pareto front of the entire MO-UCP. This front can be non-convex. Each one of the convex pareto fronts composed of the round, square or diamond points is derived from a single genotypic solution. They are composed of the phenotypic solutions found in solving the dispatching problem defined by the corresponding genotypic solution. The phenotypic solutions attached to a genotypic solution are pareto equivalent and form a convex front.

Adaptation of the fitness assignment process: The fitness value assigned to a solution u is the best fitness value among the fitness values of the phenotypic solutions p u ∈ P u :

f it(u) = opt pu∈Pu (f it(p u )) (6)
In NSGA-II the fitness is the rank of the solution, then opt is the minimization operator. This process ensures that the genotypic solution from which a pareto optimal solution can be generated is not be discarded.

Adaptation of the diversity assignment process: Let F i be the set of the phenotypic solutions of rank i. In NSGA-II, the diversity measure used is the crowding distance between a solution x and the set of the other solutions having the same fitness, i.e. the same rank:

d c (F f it(x) -{x}, x).
The adapted diversity assignment is the maximal diversity measurement among the ones of the individuals p u ∈ P u ∩ F f it(u) computed without considering the elements of

P u : d c (F f it(u) , u) = max pu∈Pu∩F f it(u) d c (F f it(u) -{P u }, p u ) (7)

Experiments and Discussion

The aim of this section is to compare the three proposed methods.

Experimental Protocol

Instances Experiments will be realised on instances of 10, 40 and 100-unit data that are generated by duplicating the unit characteristics of the ten-unit system and the demand given in Tables 1 and2. The load demands are adjusted in proportion to the size system. In all cases it is supposed that the reserve is 10% of the demand. Performance Assessment The different methods of performance assessment that can be chosen to compare multi-objectives algorithms are explained in details in [START_REF] Knowles | A Tutorial on the Performance Assessment of Stochastic Multiobjective Optimizers[END_REF]. In our case the ε-indicator and the hypervolume difference indicator are selected as they are complementary. Let Z all be the set of objective vectors from all the Pareto set approximations we obtained during all our experiments. Then, a reference set R contains the non-dominated points of Z all .

-3 -3 -1 -1 -1 ε-indicator I 1 ε+
The unary version of this indicator is computed using the binary version given by ( 8) and the reference set R, with I 1 ε+ (A) = I ε+ (A, R).

I ε+ (A, B) = inf ε∈R {∀z 1 ∈ B, ∃z 2 ∈ A, ∀i ∈ 1 . . . n, z 1 i ≤ ε + z 2 i } (8)
Hypervolume difference indicator I - H The hypervolume indicator I H is computed by the measure of the hypervolume between a set of solutions and the point z = (z 1 , . . . , z n ) where z k is the upper bound of the k th objective regarding all the solutions of Z all . The hypervolume difference indicator

I - H is then computed with I - H (A) = I H (R) -I H (A).
Experimental design All the implementations are realized under the Par-adisEO 2.0 [START_REF] Liefooghe | A software framework based on a conceptual unified model for evolutionary multiobjective optimization: Paradiseo-moeo[END_REF] software framework. A sensitivity analysis is carried out for each algorithm to determine the effect of the crossover rate and of the mutation rate.

It is done thanks to the R statistical package Irace [START_REF] López-Ibáñez | The irace package, iterated race for automatic algorithm configuration[END_REF]. The population size is fixed to 100 individuals. A convergence criteria of 100 generations without improvement of the hypervolume is used as stopping criteria. For the multi decoding embedded approach the parameter n λ is fixed to 10. For each case 20 runs are launched for each decoding system using the same seeds. Most of the performance assessment procedures are next achieved using PISA [START_REF] Bleuler | Pisa: A platform and programming language independent interface for search algorithms[END_REF] platform and its performance assessment module. The existence of a significant difference between the result obtained by the different decoding systems is verified with the Friedman statistical test. Then a post-hoc test is carried out to compare the decoders by pairs. A p-value lower than 0.005 is used as a criterion for rejecting the null hypothesis.

Experimental Results

Results for the 10-units case: In Table 5 the results of the statistical tests of comparisons are summarized. The results do not differ from one indicator to the other one. In this table and on all the following the column "Naive" indicates results of the first approach, the column "Scalarized" those of the second one and the column "multi" those of the multi decoding embedded approach. From this table we can see that for the 10-unit based case the difference between the naive approach and the scalarized approach is statistically significant. However the multi decoding embedded approach gives results significantly better than those obtained by the other approaches. Table 6 gives the best and average values obtained for each indicator and decoder over the 20 runs. The decoding embedded approach improves the ε-indicator value of 1.36% in comparison to the other method. It also improves the hypervolume indicator of 2.54%.

Results for the 40-units case: In Table 7 the results of the statistical tests of comparisons are summarized. It can be observed that the multi decoding embedded approach gives significantly better results than the one obtained with the two other approaches. The scalarized decoding approach is better than the naive one. In Tab. 8 the best and average values obtained for each indicator and decoder over the 20 runs are shown. It can be seen that that for the ε-indicator the scalarizing approach improves the results of more than 7% in comparison with the first naive approach. Then the multi decoding embedded approach improves the result of the scalarizing approach by more than 99%. For the hypervolume indicator there is an improvement of 31% when adding the λ value in the representation. Then the multi decoding embedded approach improves the result of the scalarizing approach by more than 99.5%. Results for the 100-units case: In Table 9 the results of the statistical tests of comparisons are summarized. Again, the multi decoding embedded approach gives results significantly better than the ones obtained with the two other approaches. The scalarized decoding approach is better than the naive one. In Table 10 the best and average values obtained for each indicator and decoder over the 20 runs are shown. In comparison with the naive approach, the scalarazing approach improves the results by 94.7% for the ε-indicator and by 97% for the hypervolume. Then the multi decoding embedded approach improves the result of the scalarizing approach by 75.69% for the ε-indicator and by 99% for the hypervolume-indicator .

This results and the statistical tests lead to conclude that the choice of the decoder system has a significant impact on the result. Results obtained with the last decoder are drastically better than the ones obtained with the other decoders.

Conclusion

In this article a binary genetic algorithm has been proposed to solve a multi -objectives UCP. The main difficulty is that for one genotypic solution many phenotypic solutions could be attached. These phenotypic solutions are those of the pareto front solution of the dispatching problem. Three original decoding systems associating phenotypic solutions with the genotypic ones have been presented and compared.

The multi decoding embedded system is the main contribution of this paper. The efficiency of this method has been shown on three data set of different size. In each case the results obtained are significantly better than those obtained by the two other strategies. The bigger the data, the better the improvement. This decoding system is then the one selected. In a future work, the objective will be to compare the proposed binary GA using this decoder to the GA proposed in [START_REF] Srinivasan | An evolutionary algorithm for evaluation of emission compliance options in view of the clean air act amendments[END_REF] and to a more classical multi-objective GA using a real vector to encode the solutions. First of all, it will be interesting to make an analyse study of the impact of the choice of the n λ parameter. An important advantage of this method is that it could be reused to any multi-objective problem that can be written as follow:

            
such that:

x ∈ X (10)

y ∈ Y (x) (11)
And such that for a fixed x the sub-problem P(x) finding the optimal pareto front of solutions with y ∈ Y (x) is easy to solve. In this case the genetic algorithm will handle the x variables and the phenotypic solutions are found in solving P(x). The methodology chosen to solve P(x) does not matter. This is another advantage of this decoding system. Then, to test the multi decoding embedded approach on some other problem is one of our perspectives. We also believe that the multi decoding embedded approach can be generalized to any multiobjectives genetic algorithm. Then we plan to develop a generalized version of this approach that is suitable to any genetic algotrithm.

Fig. 2 .

 2 Fig. 2. Repressentation of the genotypic solutions in the objective space.

Fig. 9 .

 9 Fig. 9. Results obtained from statistical comparisons with the ε-indicator and the hypervolume indicator.

Fig. 10 .

 10 Fig. 10. Best value obtained for each indicator and decoder over the 20 runs.

  f n (x, y)

Table 2 .

 2 Demand data with 24h time horizon.

	hour	1	2	3	4	5	6	7	8	9	10 11	12
	demand (MW) 700 750 850 950 1000 1100 1150 1200 1300 1400 1450 1500
	hour	13 14 15 16 17 18 19 20 21 22 23	24
	demand (MW) 1400 1300 1200 1050 1000 1100 1200 1400 1300 1100 900 800

Table 1 .

 1 Generating unit data for the ten-unit base system.

	unit	unit1 unit2 unit3 unit4 unit5 unit6 unit7 unit8 unit9 unit10
	Pmax(MW)	455 455 130 130 162 80	85	55	55	55
	Pmin(MW)	150 150 20	20	25	20	25	10	10	10
	a1	1000 970 700 680 450 370 480 660 665	670
	a2	16.19 17.26 16.6 16.5 19.7 22.26 27.74 25.92 27.27 27.79
	a3 (×10 -5 )	48	31 200 211 398 712 79 413 222	173
	b1	712 570 700 860 350 370 480 660 665	670
	b2	12.9 10.26 10.60 15.50 7.70 9.26 3.74 5.92 7.27	7.79
	b3(×10 -4 )	4	3	22	11	10	22	30	40	13	23
	T up min	8	8	5	5	6	3	3	1	1	1
	T down min	8	8	5	5	6	3	3	1	1	1
	CS hot	4500 50000 550 560 900 170 260 30	30	30
	CS cold	9000 10000 1100 1120 1800 340 520 60	60	60
	TCS	5	5	4	4	4	2	2	0	0	0
	initial status (h) 8	8	-5	-5	-6

  Best value obtained for each indicator and decoder over the 20 runs.

	Decoder Naive Scalarized Multi
	Naive	-	=	<
	Scalarized =	-	<
	Multi	>	>	-
	Fig. 5. Results optained from statistical comparisons with the ε-indicator and the
	hypervolume indicator.
	Indicator	I 1 ε+	I -H
	Decoder	best mean best mean
	Naive	0.736 0.738 0.451 0.455
	Scalarized 0.719 0.738 0.433 0.454
	Multi	0.709 0.712 0.422 0.425
	Fig. 6.

  Best value obtained for each indicator and decoder over the 20 runs.

	Decoder Naive Scalarized Multi
	Naive	-	<	<
	Scalirized >	-	<
	Multi	>	>	-
	Fig. 7. Results obtained from statistical comparisons with the ε-indicator and the
	hypervolume indicator.
	Indicator	I 1 ε+	I -H
	Decoder	best	mean	best	mean
	Naive	0.195	0.333	0.346	0.508
	Scalarized 0.181	0.233	0.208	0.377
	Multi	0.00129 0.0840 0.000880 0.100
	Fig. 8. Decoder Naive Scalarized Multi
	Naive	-	<	<
	Scalarized >	-	<
	Multi	>	>	-

Decoding process: In this case a pareto set P u will be associated with a genotypic solution u. This set is:

where n λ is a fixed integer and p As many phenotypic solutions are attached to a single genotypic solution, the fitness assignment and diversity assignment methods of NSGA-II have to be adapted. This will be explained in detail in the following. The decoding process is represented in Fig. 3 and Fig. 4.


Mélisande Biet
email: melisande.biet@polytech-lille.fr

Frédéric Giraud

Betty Lemaire-Semail

The implementation of tactile feedback using squeeze film effect devices

Keywords: 85, 50-n Piezoelectric devices -77, 65-j Piezoelectricity and electromechanical effects

This paper describes the implementation and initial evaluation of active-type tactile displays in desktop environments. The first device comprises a stator of an ultrasonic motor supplied by only one channel. In this way, the motor twisted of its application do not induce movement but can create slipperiness on the stator's surface. This effect can be interpreted as squeeze film phenomenon and allow us to change indirectly the friction on the surface and to generate shear forces on the fingertip, which in turn creates interesting tactile effects. Moreover, if we add a position sensor, we can vary the signal spatially which gives the sensation of bumps/ridges. Then, based on the first device principle, another device is proposed in order to enable a free exploration of the surface following ergonomic requirements.

Introduction

Artificial tactile feedback provides an under-utilized channel which can be used in virtual reality not only to augment the reality of computer graphics but also to gain information about the state of a device, the occurrence of an event, or its content. In desktop environments, the most common approach to make tactile displays is to pack a dense array of actuated pins (an early example of which can be found in [START_REF] Talbi | Vibrotactile using micromachined electromagnetic actuators array[END_REF]) and make them move orthogonally to the skin to approximate a texture. However, to simulate finely textured surface (surface unevenness scale is around 100µm) and thus enhance tactile feedback reality, this kind of tactile device is not appropriated because there are some significant limitations including the difficulty of achieving a high density array, high power consumption, and high manufacturing cost. Due to these challenges, the resulting devices are typically complex or not portable. Another possibility for simulating finely textured surfaces is the use of the squeeze film effect principle. It has been shown that a beam excited by ultrasonic vibration can give rise to a reduction of friction [START_REF] Watanabe | A method for controlling tactile sensation of surface roughness using ultrasonic vibration[END_REF]. These squeeze bearings act as a cushion of air that lifts and separates the finger from the planar display surface. By varying the thickness of the squeeze bearing, the friction on the surface can be indirectly reduced. This device has the advantage of providing a continuous surface on the fingerpulp and avoids integration problems. In this article, we first use this principle on the stator of an ultrasonic motor (part 2). After having verified that the squeeze film effect can occur on the stator ring, we add a position sensor and we modulate this effect by imposing low frequency periodic signal over the ultrasonic vibration. This allows us to show that we can vary the signal spatially and that we can give the sensation of bumps/ridges. However, even if this tactile device is more compact, has a lower weight and simulate more sensations comparatively to the Langevin's transducers used in [START_REF] Watanabe | A method for controlling tactile sensation of surface roughness using ultrasonic vibration[END_REF], it is just a one degree-of-freedom (DOF) device and besides the finger movement is small. Then, in order to solve this problem, we use a tactile plate able to generate the squeeze film effect (part 3). This plate can bend by means of PZT piezoceramics glued on one of the two surfaces. The monomorph built up in this manner can stimulate sensations similar to the previous ones but in this case in 2 directions.

A 1 DOF tactile display

The goal of this section is to present the 1 DOF tactile display and to assess if the slippiness feeling that we can feel is really due to the squeeze film effect. In order to show that, we calculate the pressure between the fingertip and the vibrating ring (stator). Afterwards, we modulate the squeeze film effect by imposing low frequency periodic signals over the ultrasonic vibration and we qualitatively evaluate the sensation obtained.

Device

We propose to use a piezoelectric motor, the Shinsei USR 60 (Fig. 1a), to generate the slippery feeling under the fingerpulp. The test device is made of the stator, and the surface of this stator can be explored directly by users without any damage (Fig. 1b). In normal working condition, the motor is supplied by two alternative phases which operates at the resonant frequency around 40 kHz. Each phase creates its own stationary wave around the stator, and the two stationary waves, in quadrature, create together a traveling one thanks to a right position of each phase ceramics and to their initial polarization. Now, in our application, we just supply one channel so that a standing wave can appeared. In this way, the finger can not be dragged along of the stator and furthermore it is impossible to feel the traveling wave to come up if the finger is put on the vibrating stator . Nevertheless, we notice that we can add a smooth feeling to the stator's surface by applying standing wave ultrasonic vibration. It may be remarked that the normal displacement of the rotor is too weak (a few µm) and that the vibration frequency is too high (40kHz) to be detected by the skin sensors [START_REF] Loomis | Handbook of Perception and Human Performance[END_REF]. For that reason, we suggest the possibility that the squeeze film effect could be the cause of the slippery feeling. Thus, in the next subsection we are going to verify this assumption.

Changing the relative friction coefficient with the squeeze film effect

Squeeze film model

In this subsection, we are describing the squeeze effect [START_REF] Salbu | Compressible squeeze films and squeeze bearings[END_REF]. Indeed, we consider the gas film created between the fingertip and the vibrating stator (Fig. 2).

We are going to rely on a study developed in [START_REF] Wiesendanger | Squeeze Film air Bearings Using Piezoelectric Bending Elements[END_REF], but taking into account the fingerprints and the flexural mode of the stator. In our case, we can not suppose that the undulations of the fingerprints are negligible relatively to the roughness of the vibrating plate (few micrometers) since the height (2he) and pitch (L) of the epidermal ridge are around 100 µm and 350 µm, respectively [START_REF] Maeno | Relationship between the structure of human finger tissue and the location of tactile receptors[END_REF] (Fig. 3). Thus, we regard the tip of the finger as an undulated surface while the stator is assumed to bend sinusoidally in a vertical direction. The film thickness, h, is the sum of the gap created when the finger skin cannot follow the ultrasonic vibration of the stator at all, which equals the amplitude of oscillation , h vib , plus the surface unevenness (average roughness), h r . h also considers the amplitude of the undulation of the fingerprint, h e .

Therefore, the thickness of the film is given by [START_REF] Watanabe | A method for controlling tactile sensation of surface roughness using ultrasonic vibration[END_REF]:

h(x, t) = (h vib +h r )[1+cos(ω 0 t) cos( 2π λ x)]+h e [1+cos( 2π L x)]
(1) In which ω 0 /2π and λ are respectively the vibrating frequency and the wavelength of the piezoelectric motor (PEM). If we normalize h, the normalized airgap can be written with non dimensional parameters as:

H = 1 + ǫ cos(T ) cos(k 1 X) + δ cos(k 2 X) (2)
with h 0 = h vib + h r , the non-dimensional parameters are

H = h h 0 + h e , X = x l 0 , T = ω 0 t, ǫ = h 0 h 0 + h e , δ = h e h 0 + h e , k 1 = 2πl 0 λ , k 2 = 2πl 0 L
where the length in contact with the fingertip is l 0 .

We make the following assumptions [START_REF] Wiesendanger | Squeeze Film air Bearings Using Piezoelectric Bending Elements[END_REF]; 1) The fluid behavior is governed by laminar viscous flow, 2) The fluid is a compressible perfect gas, 3) The inertia effect of the flow is negligible, 4) The relative lateral motion is equal to zero. Under those conditions, the one-dimensional Navierstock equation is obtained. This equation associated with continuity and ideal gas equations allows us to find the governing Reynolds equation in the non-dimensional form:

∇(H 3 P 1/n .∇(P )) = σ ∂(P 1/n H) ∂T (3)
In which n is a polytropic constant and where

P = p p 0 , σ = 12ηω 0 l 0 p 0 (h 0 + h e ) 2 (4)
p, p 0 and η represent the pressure in the gap, the surrounding gas pressure and the dynamic viscosity of the fluid respectively. The squeeze number is given by σ, which represents a measurement of the fluid compressibility in the gap. At low squeeze numbers, the fluid is nearly incompressible, while at high squeeze numbers, the fluid is trapped in the gap and acts like a spring.

We further assume that the squeezed film is isothermal (n=1). The latter assumption is reasonable since the gas film is very thin and of low heat capacity when compared with that of the vibrating surfaces. Then, to simplify the time derivation term, we substitute P H by Ψ . For steady state conditions, the integration with respect to one period yields to [START_REF] Wiesendanger | Squeeze Film air Bearings Using Piezoelectric Bending Elements[END_REF] ∇[

1 2 H∇(Ψ 2 ∞ ) -Ψ 2 ∞ ∇(H)] = 0 (5)
where Ψ ∞ denotes the Ψ inside the airgap when σ → ∞ and H is the mean normalized film thickness given by

H = 1 + δ cos(k 2 X) (6)
It may be noticed that the analytical solution is accurate as far as σ is assumed at a very large value. On the contrary, a low value of σ would imply a much more complex solution.

Approximate analytical solution

The equation 5 gives us the relationship between the normalized film thickness and the normalized pressure for a fluid behavior close to a spring. The solution of ( 5) is given by (Appendix A) [START_REF] Wiesendanger | Squeeze Film air Bearings Using Piezoelectric Bending Elements[END_REF]:

Ψ 2 ∞ = K 1 [1 + δ cos(k 2 X)] 2 (7)
In which K 1 is a constant. The boundary conditions with respect to Ψ 2 ∞ (X), noted down as Ψ ∞B are

Ψ 2 ∞B (- 1 2 ) = Ψ 2 ∞B ( 1 2 ) = K 1 [1 + δ cos( k 2 2 )] 2 (8)
In which 1 2 is the normalized value of l0 2 . By replacing K 1 by its expression in [START_REF] Smith | Deployment of fingertip forces in tactile exploration[END_REF], we obtain:

Ψ ∞ = Ψ ∞B ( 1 2 ) |1 + δ cos(k 2 X)| |1 + δ cos( k2 2 )| (9)
Furthermore, in order to find Ψ ∞B , we will focus our attention on the boundary region. Considering that p 0 is constant in time and that Ψ ∞B equals Ψ ∞ when the boundary region meets the interior of the airgap, we obtain [START_REF] Wiesendanger | Squeeze Film air Bearings Using Piezoelectric Bending Elements[END_REF]:

Ψ 2 ∞B ( 1 2 ) = p 2 0 T +2π T H 3 B dT T +2π T H B dT = p 2 0 T +2π T [1 + ǫ cos(T ) cos( k1 2 ) + δ cos( k2 2 )] 3 dT T +2π T [1 + ǫ cos(T ) cos( k1 2 ) + δ cos( k2 2 )] dT = p 2 0 [1 + δ cos( k 2 2 )] 2 1 + 3 2 ǫ 2 [1 + δ cos( k2
2 )] 2 (10) By replacing Ψ ∞B by its expression in [START_REF] Biet | Using an ultrasonic transducer: Evidence for an anisotropic deprivation of frictional cues in microtexture perception[END_REF], we obtain:

P ∞ = Ψ ∞ H = p 0 (1 + δ cos(k 2 X)) (1 + δ cos( k2 2 )) 2 + 3 2 ǫ 2 cos 2 ( k1 2 ) (1 + δ cos( k2 2 ))(1 + ǫ cos(T ) cos( k1 2 ) + δ cos(k 2 X)) (11)
Remember that P ∞ is the normalized pressure inside the airgap for an infinite value of the squeeze number σ. Then, the mean squeeze pressure is expressed by:

P ∞ = p 0 (1 + δ cos(k 2 X)) (1 + δ cos( k2 2 )) 2 + 3 2 ǫ 2 cos 2 ( k1 2 ) (1 + δ cos( k2 2 )) (1 + δ cos(k 2 X)) 2 -ǫ 2 cos( k1 2 ) (12)

Results

First of all, let us remind you that in order to predetermine analytically the over-pressure on the fingertip, we have to satisfy the following condition : σ → ∞. However, we assume in practice that the squeeze pressure (P ∞ ) depends almost entirely on the amplitude of vibration when the squeeze number is larger than 10. Therefore, we calculate the squeeze number σ as shown in (4), with the parameters of Table 1. We obtain σ = 20.4 for h vib = 2µm. This result enables us to use the equation ( 5) and to calculate the overpressure between the finger and the stator.

Using the parameters shown in Table .I, we can compute the pressure profile along the airgap using (Eq.11).

On the upper part of figure 4, the finger at contact is modeled as a sinusoid with an amplitude of 50µm in order to take into account the fingerprints. The interval from the zero depicts mean roughness plus the vibration amplitude of the vibrating plate itself. The curve below shows the evolution of the pressure temporal mean as a function of the position on X-axis. In this figure 4, pressure peaks are localized where the finger-skin is as close as possible to the vibrating plate, i.e. at the fingerprint ridges, whereas at the fingerprint grooves, pressure is not far from the atmospheric pressure. The pressure P f is the mean finger pressure that it is used during an exploration task.

P f = F f /l 2 0 ( 13
)
where F f denotes the normal contact force applied by a person exploring the surface. For its value, we choose 0.5 N as a mean normal contact force according to [2][7]. We can notice from Fig. 4 that there are some zones of the fingerprint where the squeeze pressure is superior to the finger pressure. Those zones are not in contact with the vibrating plate. Moreover, the mean pressure P moy (Fig. 4) is higher than 1, which means that the squeeze film effect is significant for ultrasonic frequency contact of the stator and the finger, as it has been previously shown for the stator and the rotor of the USR 60 [START_REF] Maeno | Effect of the hydrodynamic bearing on rotor/stator contact in a ring-type ultrasonic motor[END_REF].

Since the squeeze force occurs, we can express the relative coefficient of friction, µ ′ µ in terms of (Eq.14).

µ ′ µ = 1 - (P ∞ -1) P f (14)
Physically, and following the assumptions previously defined, this equation means a weakening of the friction coefficient for a given surface when the squeeze force occurs.

To illustrate the friction coefficient variations, we can compute µ ′ µ . The results are given in Fig. 5. We can easily calculate the relative friction coefficient, taking into account the fingerprints and the flexural mode of the stator. Moreover, those results show that surface roughness, contact area and the vibration amplitude are important parameters of the squeeze film effect. We have also to notice that the relative coefficient of friction is strongly dependent on the force applied by the user (F f ), thus the results of (Fig. 5) have to be considered with care. The previous calculation is performed by assuming that the stator surface is flat and without teeth, which increase the contact area between the finger and the stator. However, owing to the ridge width is really larger than the groove width of the teeth, we assume that this choice has no important effect on the present study. Furthermore, considering the Figure 5, a vibration amplitude of the stator higher than 7µm would be necessary to levitate the finger, i.e. to reach a zero relative friction coefficient. This value could appear very high but it is not necessary to reach it because a vibration amplitude of 2µm is enough to feel a difference of sensation. Hence, by lubrication effect (squeeze film effect), we can succeed in changing the sensation by exploiting the friction coefficient between the finger and the plate as it was qualitatively checked in [START_REF] Biet | Using an ultrasonic transducer: Evidence for an anisotropic deprivation of frictional cues in microtexture perception[END_REF].

Space period variation for roughness simulation

Principle

Spatially varying voltage signals are effective at conveying a true sense of texture. Using position feedback (Fig. 6), it is fairly easy to make certain locations feel smoother than others. For that aim, we impose a periodic signal over the standing wave at low frequency. This comes down to a disturbance injection or a wave amplitude modulation. As a result, the squeeze force is also modulated as well as the tangential speed. Thus, controlling the spatial frequency between slippery area leads to grooved textures because of the alternation between slippery and braking feeling. Many possibilities can be suggested following this princi- ple: about the shape of the injected signals, about their frequency, and whether or not this frequency is tuned following the finger speed. For example we can use a square signal that turns off and on (pulse train), sine wave, localized pulses or velocity dependent signals.

Experimental device

The experimental device is given in figure 7. As we need the finger speed and direction, a LVDT sensor is used. The connection to the finger is rigid in order to avoid some inconsistency in of all of based effects. Only linear displacements are allowed, which implied that the finger movement is quite small (≈ 4cm). Nevertheless, ahead or reverse motion may be done. The stator of the PEM is supply by only one channel so as to create a standing wave on the stator. Its maximum vibration amplitude measured with the LASER Doppler Vibrometer is 2µm. Then, to add to the friction range, we apply a flat and braking adhesive on the surface of the stator (h r = 0.6µm instead of 0.2µm for the stator by itself [START_REF] Maeno | Effect of the hydrodynamic bearing on rotor/stator contact in a ring-type ultrasonic motor[END_REF]). According to Fig. 5, the choice of the adhesive is crucial for a good working of the tactile device since the adhesive unevenness has a strong incidence on the squeeze film effect. The device is run in a real time environment and textures are created in software as voltage patterns. In actual use, dragging with his finger, the user can feel the sensation of roughness thanks to a wave amplitude modulation. Moreover, textures with sharp changes in friction are much stronger than those with smooth transitions. For this reason, we use square pulses.

Evaluation of the device

If we change the space period, we can easily recognize whether the surface is finely or coarsely grooved. However, we need to evaluate this device to know if the space period we simulate correspond to the space period of a real gratings. Using a rectangular plate where several types of grooved surface are engraved, we carried out a match to sample experiment by comparing this plate and the 1 DOF device.

The experiment is conducted with twelve naive students aged between 18 and 25 years. Gratings are single sided copper printed circuit boards (PCBs). The PCBs measure 7 by 130mm and are etched with striated layouts as shown in figure 8. The copper thickness or height of ridge is equal to 70µm. The groove width and the ridge width of the gratings are respectively equal to 0.4, 2, 4, 8 and 10mm. This gratings are hidden in a box so as to not be influenced by vision. The subject's task is put his/her index finger in the sen-8. Gratings used for the match to sample experiment.

sor ring and to form a straight line on the stator surface with a backward and forward motion of his/her finger to discover the virtual surface. Then he/she has to compare his/her sensation with the gratings in the box (fig. 9). To improve the match to sample procedure, gratings are first randomly presented to participants. After this preliminary phase, subjects has to evaluate the "virtual surface". They are asked to say which gratings looks like the "virtual surface" tested. The test result is shown in Table . It can be seen from the results that users are able to distinguish the range of the spatial period. Difference of the displayed sensation was clear except for the last case (0.4 mm).

Limitations and improvements

This experiment is a preliminary evaluation, which shows that the 1 DOF tactile display can produce several distinct sensations. Effectively, several periodic rough textures can be simulated, with a low frequency modulation. The modulation signal has also to be tuned according to the fingers speed. However, even if this tactile device is more compact, has a lower weight and simulate more sensations comparatively to the Langevin's transducers used in [START_REF] Watanabe | A method for controlling tactile sensation of surface roughness using ultrasonic vibration[END_REF], it is just a one degree-of-freedom (DOF) device and besides the finger movement is small. Thus, due to the small one-dimensional workspace of the 1 DOF tactile display, the range of texture sensations is consequently limited. A larger, two-dimensional array able to generate the squeeze film effect has to be designed in order to solve this problem. This new device should be able to improved the sensation and to generate new modes of texture with spatial signals that can vary in two dimensions.

The limitation concerns finger orientation. In the course of testing, we found that there was a significant difference in the intensity of the perceived texture depending on whether the finger touched the plates with a horizontal or vertical orientation. Thus, it is important that the user interacts with the device using a consistent finger orientation [START_REF] Martinot | Mechanical responses of the fingerpad and distal phalanx to friction of a grooved surface: Effect of the contact angle[END_REF].

Roughness simulation on a 2 dof device

In this section we are going to design a plate which will be able to provide the same sensations but with a greater workspace, adapted to the amplitudes of movement in active touch condition. This plate can bend by means of PZT piezoceramics glued on one of the two surfaces. The monomorph dimensions are going to be chosen in order to satisfy ergonomic and squeeze film effect requirements. Then, we proceed to experimental verifications.

Design of the monomorph

We have at our disposal passive materials in Copper Beryllium, and a piezoelectric ceramics PZT, referred to as PI-91, marketed by Saint-Gobain Quartz company. The characteristics of the material are indicated in Table 4. The

Mechanical properties of the resonator

Young modulus Ei (10 9 N.m -2 ) 123

Mechanical properties of P1-91 ceramic

Piezoelectric constant e31 p (C.m -2 ) -4.9 Elastic constant c E 11p = s E 11 (s E 11 ) 2 -(s E 12 ) 2 (10 10 N.m -2 ) 6.79 Charge coefficient d31 (10 -12 m.V -1 )
-247

Mechanical quality factor

Qm 60

Table 3. List of parameters voltage supply of piezoactive ceramics is set to 15V in order to avoid a dangerous voltage. Moreover, the thickness of the piezoceramic is 1 mm, while the thickness of the substrate is set to 2 mm (this value has been chosen to guarantee mechanical holding and fabrication considerations).

An ergonomic workspace

The compatibility of our tactile device with the user's movement does not have to exclusively deserve the mechanisms of stimulation. It also has to offer perspectives on tangible interaction with computer-generated surfaces for example. Therefore, the definition of the ergonomic workspace plays an important role for the use of such devices. For finely textured surfaces, it is recommended to "allow the freedom of active exploration" [START_REF] Pasquero | Stress: A practical tactile display system with one millimeter spatial resolution and 700 hz refresh rate[END_REF]. For that purpose, the biggest exploration surface needs to be considered for the design of the tactile plate. There is no exact dimension of exploration surface in literature, but it is known, according to the observations of Klatzky and Lederman, that typically the finger quickly rubs back and forth along a small, homogeneous area of the surface, and that interior surfaces are explored rather than edges [START_REF] Lederman | Hand movements : a window into haptic object recognition[END_REF]. Moreover, useful information in measurements of movement amplitudes are reported in [START_REF] Martinot | Caractérisation du rôle de la dynamique du toucher dans la perception de textures[END_REF]. In his experiment, this author collects the position and normal force exerted by the fingertip while volunteers are actively exploring finely textured surfaces. The values for the movement amplitude along the X and Y directions are 97.27mm and 84.38mm. The mean and the biggest rectangular areas that are touched measure respectively 28.18cm 2 and 128.37 cm 2 . Within those guidelines, we can choose the width, b, of the plate in order to have a surface area higher than 28.18cm 2 . Thus, we choose b = 50mm.

Squeeze film effect requirements

Now that we know the thickness and the width of the monomorph, the length has still to be found. This choice is the most important for a good working of the device because the flexural wave propagate along the length. The choice relies on a condition on the frequency and on a condition relative to coefficient of friction at the same time.

a/ Condition on the frequency

The first condition is done by the squeeze number, σ, which has be higher than 10 as we previously mention. We calculate the squeeze number σ as shown in ( 4), with the parameters of Table 1 as a function of the frequency of the oscillation for a given amplitude h vib . Considering the most unfavorable case, i.e. the smallest slope of the straight line representing σ = f (ω 0 ) (a surface state, h r , of 1.6 µm and a vibration amplitude, h vib , of 3µm), the results show a first criterion on the frequency (f >20.6 kHz) which enables us to use the equation [START_REF] Wiesendanger | Squeeze Film air Bearings Using Piezoelectric Bending Elements[END_REF].

Furthermore, the resonance frequency, f n , is a function of the plate dimensions and can be expressed by the following equation:

f n = π λ/2 2 G b M b (15)
in which M b is the total mass per length, expressed as follow:

M b = ρ p h p + ρ i h i (16)
in which G b is the flexional rigidity of the monomorph in N.m. G b is the sum of the flexional rigidity of both materials where z here is measured from the neutral axis.

G b = c E 11p b piezo z 2 dz + E i b substrate z 2 dz (17)
In Fig. 10, the evolution of the resonant frequency obtained from (Eq. 15), f n , is represented according to the half wavelength. To match the condition on the frequency according to the squeeze number, we need f n > 20.6kHz, which implies that the half wavelength needs to be smaller than 16.1 mm. Thus, we obtain with this study an upper boundary for the wavelength.

b/ Condition on the relative coefficient of friction

According to the study on the stator of the PEM (Fig. 5), we need a relative coefficient of friction µ ′ /µ lower or equal to 0.7 in order to feel a difference of sensation. We calculate the relative coefficient of friction with the same parameters and for 4 half-wavelengths lower than 16.1 mm (Fig. 11). In this manner, we have 4 vibration amplitudes corresponding to this value. Considering the Fig. 11, it is necessary to obtain vibration amplitudes of 1.4, 1.55, 1.9 and 2.3µm if we choose half-wavelength of 7.5, 10, 12.5 and 15mm respectively. As we can see, a small wavelength will impose a high vibration amplitude and we will be limited quite rapidly because an vibration amplitude higher than 2µm is difficult to reach. For that aim, we choose λ/2 ≥ 10mm. Then, following the requirement for the design, we choose a surface length L = 83mm in order to put 7 ceramics measuring 11mm isolated from each other by a distance of 1mm. The final dimensions of the plate are: L × b = 83mm × 49mm. This gives us a surface area of 40.67cm 2 .

Numerical study on the entire monomorph

The FE model of the monomorph is represented in Figure 12. A modal analysis is performed considering the entire mono-morph. We find the (8; 0) mode which corresponds to a half-wavelength of 10.37mm in the X direction. The modal shape is presented in Fig. 13 and proves that this resonance frequency (35.6kHz) is in concordance with the squeeze analysis ( >20.6kHz). Afterward, we check if those dimensions meet amplitude requirements thanks to an harmonic analysis. The ceramics are activated by two electrical signals, with a 180 o phase shift between each which creates a standing wave (Fig. 13). The harmonic response of the actuator gives the vertical displacement as a function of the frequency (Fig. 14) and shows that the natural frequency is located at 34.77kHz. With this frequency, a vertical displacement of is obtained at the top of standing wave.

Experimental verifications

The prototype is presented in Fig. 15. The polarity of each ceramic is oriented to satisfy the wave production along the X direction. The device is supplied by one voltage source adjusted to a mechanical resonance frequency so as to generate a standing wave along the length of the plate. Some experiments with the prototype built at the laboratory have been performed to confirm analytical and numerical studies.

Vibration amplitude measurement

The vibration amplitude is measured using a single-point LASER Doppler Vibrometer (OFV505) linked to a controller (OFV-5000) that is connected to an oscilloscope.

As the results show, a deflection amplitude of about 2.3µm peak to peak is obtained by applying a voltage of 15V. The resonant frequency is 30.5kHz, which gives a squeeze number of 15.9 for h r = 0.6µm. The vibration amplitude of the prototype is smaller than the predicted amplitude. This amplitude attenuation may be due to the manufacturing of the monomorph and especially to the gluing phase. More precisely, the gluing phase is crucial for the vibration quality owing to the ideal gluing hypothesis. However, according to the relative friction study, this vibration amplitude of the plate can be enough to feel the difference of perception. This assumption will be checked in the next section.

Qualitative study of the sensation

Conclusion

This paper described a new, modular, high-performance, tactile transducer device based on the squeeze film effect principle. It had 60 actuators creating an active surface of about 10 mm2 in the shape of a cradle designed to minimize the finger pre-load when the tactile signals is quiescent. The system is compact (150 cm2, 60 g) and has integrated electronics requiring a small number of connection wires for operation. A prominent feature the device is the ease with which is can be fabricated and serviced. The results indicated that each actuator could produce a free deflection of 0.1 mm and a blocked force of 0.15 N (or 0.075 mm and In this article, we have presented a complete modeling of a planar standing wave ultrasonic motor. Firstly, the vibratory behavior has been analytically introduced with the variational formulation. Then, an experimental iden-tication method is presented and allows to identify some nonlinearities of the freely vibrating plate. The experi-mental results are compared to the analytical ones. This approach is validated for a single vibratory mode, but may be extended for all others, indiscriminately along the length or width of the plate. Secondly, the mechanical conversion is introduced. This part of modeling is at the origin of important non-linearities because of mechanical asymmetry, possibly increased by the intermittent contact phenomenon. Nevertheless, the overall behavior can be expressed from simple relations, since the inertia of the actuator lters the small vibratory transitions of the contact. By this way, an average modeling of the tangential motion is established. The tangential results are strongly dependent on the motion along the normal axis, mainly induced by the intermittent contact. This inter-mittent contact varies according to the normal pre-load or also with the wave amplitude. Some experimental mea-surements were sucient to emphasize the overall behav-ior and validate the simplied mechanical approach. At last, the complete identied modelling is discussed, thanks to causal ordering graph formalism. It will be useful for deducing a force control law, as presented in [?].

Fig. 1 .

 1 Fig. 1. a. The Piezo-electric motor(PEM) b. The stator of the PEM is the surface to touch.

Fig. 2 .

 2 Fig. 2. Mechanism of slippage.

Fig. 3 .

 3 Fig. 3. Approximate profile of a fingertip on the vibrating stator when the epidermal ridges are taken into account.

Fig. 4 .

 4 Fig. 4. Thickness of the film at a given time and spatial profile of the mean time pressure under the finger (hr = 0.6µm, h vib = 2µm)

Fig. 5 .

 5 Fig. 5. The relative coefficient of friction as a function of the amplitude of vibration for F f = 0.5N and for hr = 0.2µm and hr = 0.6µm.

Fig. 6 .

 6 Fig. 6. Modulation of the squeeze film effect.

Fig. 7 .

 7 Fig. 7. The experimental device.

Fig. 9 .

 9 Fig. 9. Experimental Setup.

Fig. 10 .

 10 Fig. 10. Natural frequency as a function of the half wavelength.

Fig. 11 .

 11 Fig. 11. The relative coefficient of friction as a function of the amplitude of vibration for F f = 0.5N , hr = 0.6µm and for 4 values of λ/2.

Fig. 12 .

 12 Fig. 12. FE model of the monomorph.

Fig. 13 .

 13 Fig. 13. Result of the modal analysis: deformed shape of the (8;0) mode at f=35.6kHz.

Fig. 14 .

 14 Fig. 14. Result of the harmonic analysis: Vertical displacement of a surface point.

Fig.

  Fig. Complete view of the below side with a piezo-ceramics matrix and above side which is the surface to touch.

Fig. 16 .

 16 Fig. 16. Laser vibrometric measurements of the (8; 0) mode.

Table 1 .

 1 List of parameters

	Amplitude of the fingerprint	he	50 [µm]
	Period of the fingerprint	L	350 [µm]
	Wavelength of the PEM	λ	18.7 [mm]
	Resonance frequency of the PEM

Table 2 .

 2 Percentage of answer

	Space period of the gratings	Space period simulated 0.4 2 4 8 10 mm mm mm mm mm
	0.4 mm
	2 mm
	4 mm
	8 mm
	10 mm

Table 4 .

 4 Percentage of answer

	Space period of the gratings	Space period simulated 0.4 2 4 8 10 mm mm mm mm mm
	0.4 mm
	2 mm
	4 mm
	8 mm
	10 mm


Frédéric Giraud

Betty Lemaire-Semail

Julien Aragones

Jacques Robineau

Jean-Thierry Audren

Stability Analysis of an Ultrasonic Motor for a new wave amplitude control

Using Piezo electric actuators can reduce the bulk size of servomechanisms; they are thus very useful in avionics applications. However, mechanical overload on the shaft of a Traveling Wave Ultrasonic Motor often results in the motor suddenly stalling. To avoid this drawback, one can increase the supply voltage or add a control loop in the rotating reference frame of the traveling wave. The consequences are extra power losses or lower dynamic performances.

The proposed method combines the advantages of classical controls and controls in a rotating frame: both stability and dynamic performances are obtained at low supply voltage levels. Experimental runs are presented.

I. INTRODUCTION

Traveling Wave Ultrasonic Motors (TWUM) exploit a piezoelectrically generated flexural wave which propagates at the surface of a stator. This wave is able to propel by contact a rotor strongly pressed on it. Friction produced by the contact mechanism is at the origin of numerous advantages, among which we can find: braking without supplying voltages, a high torque to mass ratio and high torque -low speed characteristics. Thus, while a speed reducer is often needed with an electromagnetic motor, it becomes useless in applications using TWUM: this leads to lightweight and compact applications. These features greatly increase the interest of these motors for servo applications in the avionics industry. Unfortunately, they are very difficult to control. First, since torque generation is non linear, it is difficult to precisely drive a load. However, solutions have been provided to solve that problem [START_REF] Senjyu | Position control of ultrasonic motors using mrac and dead-zone compensation with fuzzy inference[END_REF] and achieve precise position control in optoelectronic applications [START_REF] Giraud | Precise position control of a travelling wave ultrasonic motor[END_REF]. Secondly, the resonant behavior of the stator is very sensitive to different external conditions: temperature, load torque or mounting conditions. For example, the motor may stall if a mechanical overload is applied to the shaft, and this seriously decreases the reliability of servo mechanisms actuated by TWUM. Reducing the motor's sensitivity to external torque variations is thus a key issue for their use in avionic applications. In fact, bumps are typical of this kind of environment, where random disturbances are applied to the load and extremely accurate position or speed controls are needed. To solve this problem, one can increase the supply voltage. This displaces the torque limit which makes the motor stall, but does not remove the problem completely. Moreover, it increases the power losses and thus decreases the efficiency of the motor. Other solutions are proposed in [START_REF] Maas | Optimized drive control for inverterfed ultrasonic motors[END_REF] [START_REF] Giraud | Analysis and phase control of a piezo electric travelling wave ultrasonic motor for haptic stick application[END_REF]; [START_REF] Giraud | Analysis and phase control of a piezo electric travelling wave ultrasonic motor for haptic stick application[END_REF] is based on a control of the motor in a rotating reference frame fixed to the traveling wave. By controlling the phase shift between the supply voltage and the traveling wave, it is possible to ensure operation for any load conditions. However, this solution slows down the traveling wave's dynamic. When implemented in a position control, this leads to lower dynamic performances compared to a classical frequency control. The aim of this article is to describe a new control method which helps to improve the reliability of the motor's operation at low supply voltage levels while maintaining good dynamic performances during transitory operations. For that purpose, an overview of method used to model these motors is given. Then the proposed control is detailed. Experimental results are given in the final section for an industrial optoelectronic application.

II. MODELLING OVERVIEW A. Modelling in the αβ frame.

A TWUM is made up of a ring shaped stator on which piezo-electric elements are bonded. If these elements are supplied close to the resonant frequency of this mechanical set, a bending wave is propagated in the stator. This travelling wave is able to propel by friction a rotor pressed on it (figure 1).

Different authors [5][6] have described the energy conversion process for TWUM. This process can be divided into three stages. In the first stage, the supply voltages named v α and v β create forces via the inverse piezoelectric effect. In the second stage, the stator vibrates because of these forces. The deformation w along the z axis can be written as 1

w(θ, t) = w α (t) cos(θ) + w β (t) sin(θ) (1)
where the w α and w β are the stationary wave's amplitudes, which occur if v α and v β are supplied alone respectively, and k is the mode number. In figure 1 the stator is in contact with the rotor at one contact point along the wavelength. The wave crest location θ c can be deduced from w α and w β by [START_REF] Giraud | Causal modeling and identification of a traveling wave ultrasonic motor[END_REF]: The relationship between w α , w β and v α ,v β can be described using a 2 nd order equation

kθ c = tan -1 w β w α (2)
m ẅα + d s ẇα + cw α = N v α -f rα (3) m ẅβ + d s ẇβ + cw β = N v β -f rβ (4)
where m,c, d s and N are parameters of the stator; [START_REF] Giraud | Causal modeling and identification of a traveling wave ultrasonic motor[END_REF]and [START_REF] Giraud | A torque estimator for a traveling wave ultrasonic motor -application to an active claw[END_REF] describe the protocol to identify these parameters. The resonant frequency of the stator is given by ω 0 = c m . In equation 3, f rα and f rβ take into account the external forces acting on the stator which are:

• the normal force pressing the rotor on the stator F N • the load torque T Under the assumption of ideal contact (no sliding and conditions of ideal punctual contact), it is possible to deduce a relationship between f rα , f rβ , F N , and T :

f rα f rβ = R(kθ c ) F N k h b 2 T ( 5
)
where h is the thickness of the stator and b is the radius of the stator. R(kθ c ) is the well known rotational matrix, often used when modelling electromagnetic motors:

R(θ) = cos(θ) -sin(θ) sin(θ) cos(θ) (6)
The third stage of the energy conversion process is the friction between the stator and the rotor which produces the torque appearing on the shaft as a function of the wave's amplitude which is then the key variable controlling the motor.

B. Complex phasors.

It is quite common to use complex phasors in Electromagnetic machines with a rotating magnetic flux [START_REF] Vas | Vector control of AC machines[END_REF]; they are useful in obtaining simple representations of electrical variables ( voltage or current) of the machine at steady state conditions, and in finding their values. With complex phasors, we focus our attention on RMS or average values of the key variables.

In this section, we apply concept of complex phasors to TWUM. We attempt to construct a simple representation which will help to deduce the wave's amplitude as a function of the supply voltage's amplitude and frequency. Complex phasors rely on the concept of rotating reference frames. In fact, instant values of w α and w β can be respectively deduced from the imaginary part and the real part of a complex number w. And because w α and w β are sinusoidal functions of time, the edge of w appears to rotate in the complex plane. The same is true for v α and v β which can be deduced from v. This concept of a rotating reference frame is summarized in figure 2 As such, we define the following complex variables:

w = w α + jw β (7) v = v α + jv β ( 8
)
At this stage, we can deduce from 2 that:

w = w 2 α + w 2 β e jkθ c ( 9
)
and define

W = w 2 α + w 2 β ( 10
)
where W is the wave's amplitude; 10 leads to:

w = W e jkθc (11)
Moreover, defining f r , and using 5 and 6 we can write:

f r = f rα + jf rβ = (F N + jk h b 2 T )e jkθ c ( 12
)
The relationship between the complex phasors presented in this section is deduced from 3 and 4. In fact, (3) + j × (4) leads to:

m ẅ + d s ẇ + cw = N v -f r ( 13
)
To calculate the first and second derivatives of w, we use 11 which gives:

ẇ = d dt (W e jkθ c ) = Ẇ e jkθ c + jk θc W e jkθ c ( 14
)
When using complex phasors, we consider steady state only. Given this condition, Ẇ is obviously equal to 0, since it is the amplitude of the wave in the stator. Moreover, for a perfect travelling wave propagating in the stator, the velocity of the contact point is given by the frequency ω of the two voltages, and we have k θc = ω; ω is supposed to be constant for steadystate operation. Finally, 14 can be revised to take into account the steady-state assumption:

ẇ = jωW e jωt (15)
This gives ẅ = -ω 2 W e jωt (16)

Using 12, 15 and 16 to rewrite 13 gives

(-mω 2 W + jd s ωW + cW )e j ωt = N v -(F N + jk h b 2 T )e jωt (17) In order to simplify (17), we let v = (V d + jV q )e jωt ( 18
)
We can see that V 2 d + V 2 q = V , where V is the voltage's amplitude. Finally, 17 can be rewritten as: 19) is of prime importance because it allows us to find the amplitude of the vibration of the stator, as a function of external conditions (load torque T and normal preload F N ) and of the supply conditions (V and ω). It can be divided into real and imaginary parts, leading to two new equations:

-mω 2 W + jd s ωW + cW = N (V d + jV q ) -(F N + jk h b 2 T ) (19) Equation (
(c -mω 2 )W = N V d -F N ( 20
)
d s ωW = N V q -k h b 2 T (21)
This can be summarized as shown in figure 3. In this figure, the d -axis is the real axis, the q -axis is the imaginary one. Moreover, we assume that the amplitude of the supply voltages is constant, but the frequency varies. So the edge of underlinev moves along a circle of radius V .

Under no-load conditions (T = 0), and according to equation (21), W can be deduced from the projection of phasor v on the q -axis, while V d satisfies (20). As ω approaches the resonant frequency of the motor, V d decreases and V q increases: W thus increases. The maximum value for W is thus attained if the phasor V is fully projected on the q -axis. For this operating point V d = 0.

If the motor is supplied below its resonant frequency, the motor stalls. This is well described in [START_REF] Giraud | Analysis and phase control of a piezo electric travelling wave ultrasonic motor for haptic stick application[END_REF] and [START_REF] Maas | Controlled ultrasonic motor for servodrive applications[END_REF]. However, the stability of the motor under no-load conditions is not detailed in this article.

We can measure Ψ, the angle between v and W . One can see that Ψ decreases as the motor approaches resonance, reaching Ψ = π 2 for that point. It is possible to measure Ψ during experiments. In fact, as is the case for complex phasors, the angle measured in a complex plane is equivalent to a delay measured using a scope. In the case of TWUM, we have to measure the delay between the supply voltage v α and w α (or between v β and w β ). In this article, we give experimental results obtained with a Shinsei USR60 [START_REF]Internet website[END_REF] A sensor is glued to the stator of the TWUM, which provides a voltage proportional to the stator's deformation which we call v EA . Unfortunately, the sensor is not at the right place, and an additional delay of 254 o has to be taken into account between the sensor's voltage and w α .

Figure 4 shows values for v α , v β and v EA for steady state conditions, with W = 2µm. A phase shift of 98 o is measured, leading to Ψ = 156 o when taking into account the location of the sensor. The corresponding complex phasors are depicted in figure 5. Figure 6 presents other values at 3µm. We now measure a phase shift of 114 o , which gives Ψ = 140 o . These results are consistent with the analysis given in figure 3, because Ψ decreases as W decreases. They show that complex phasors are useful to describe the motor's operating point. In the next section, we use them to analyze stability of the motor when we vary the external load torque.

III. STABILITY ANALYSIS AS A FUNCTION OF EXTERNAL LOAD TORQUE

A. Effect of a load torque

If a load torque is applied on the shaft as depicted in figure 8, V q -the projection of v onto the q -axis -is divided into two parts: one part for the wave propagation and one part for the load torque T . For a given wave amplitude W , if the load torque increases, then v gradually approaches the q -axis. Once again, this can be experimentally verified as shown by the values in figure 9. For this run, W is equal to 2µm as in figure 4, but now, an external load torque of T = 0.5N m is applied to the shaft. As a consequence, the measured phase shift between v EA and v α is equal to 120 o , leading to Ψ = 126 o : v approaches the q -axis as T increases. From figure 8, we can deduce that the maximum load torque which allows a stable operation of the motor is then given in the case when v is aligned on the q -axis. At this point, if T becomes larger, the motor stalls, because it is not able to develop a higher torque: equation (21) has no solution. So, stable operations are ensured if

Ψ < π 2 (22)
is verifed.

B. Safe operating area.

The safe operating area of a TWUM is a zone delimited by three boundaries defined as follows:

• W must be lower than a mechanical limit beyond which the bond layer would be destroyed. We call this limit W M ax , • T cannot be larger than a limit defined by tribological contact conditions between the stator and the rotor. We call this limit T M ax , • The motor should not stall. This defines the maximum operating torque, called T lim . T lim is a function of V and W . In fact, starting from equation ( 21) we can deduce that the maximum available torque is given when V q = V , leading to:

k h b 2 T lim = N V -d s ωW (23)
In figure 11, we have drawn the boundary conditions of the safe operating area in the (T,W) plane, which seems to be more appropriate, for two values of the supply voltage

V (V 1 > V 2).
V = V 1 and V = V 2, V 1 > V 2.
Figure 11 explains why, during operation, the motor may stall if extra load torque is applied to the shaft of the motor. The maximum load torque depends on the operating point itself: bigger W is, the less torque we have available. Moreover, we also understand why increasing the motor's voltage is a solution to the problem: in fact, this increases the maximum available torque for each travelling wave's amplitude. However, this solution is not effective because it leads to additional power losses in the motor and the production of heat.

C. Experimental measurement of the safe operating area.

It is possible to measure T lim as a function of W and V . To do so, we increase the load torque on the shaft, by using a powder brake for example, and ω is adjusted so as to keep W constant. During this test, we also measure the phase shift between v α and v EA , so as to estimate Ψ. All the variables were recorded. We noted the values of T and Ψ when the motor suddenly stalled (W failed to zero then) even though it was previously rotating normally. Runs depicted in figure 12 Fig. 12. Safe Operating Area for a Shinsei USR60 at V = 150V and V = 110V .

show that the maximum available torque T lim actually depends on W and V . However, we have found curves which are not equal to the modelling presented in this article. Of course, errors may be due to measurement errors of the load torque, or to the protocol (it is difficult to know exactly whether the motor has stalled or not). It may however be due to the modelling process which does not accurately take into account the specific contact conditions between the stator and the rotor.

Moreover, we do not verify in figure 13 that Ψ = π 2 when the motor stalls; at present, we cannot find a reason why such a difference occurs. However, because Ψ can be easily measured for the motor, while we do not plan to measure the load torque in our application, keeping Ψ > 120 seems to be a good solution to avoid the motor stalling at any voltage. This method will be exploited in a control detailed in the next section.

IV. WAVE AMPLITUDE CONTROLS OF TWUM.

In TWUM, the key variable controlling the torque and the speed of the motor is W . This is why we should always find a Wave amplitude control loop in control schemes of TWUM. This section describes classical and new control strategies.

A. Frequency control

In classical controls, the frequency ω of the supply voltages is controlled so as to obtain a reference value of W , as shown in figure 14. This frequency control is straightforward, but is also very sensitive to the external conditions. In this control, W is set independently of the load torque, and if a load larger than T lim is applied, the motor stalls. To start up again, a long and specific procedure has to be followed: this closed loop control does not work in our case. In order to increase the torque limit, increasing V is useful, but does not completely remove the problem. Moreover, it is difficult to precisely tune the loop controler, because the relationship between W and ω is not linear. This control however is made up of a single loop, leading to excellent performances if we consider the response time obtained.

B. Control in a rotating reference frame

In addition to increasing the voltage, we could also adapt the reference value W ref shown in figure 14 when large load torques are applied to the shaft as described in figure 11: from the measurement of T , W is adapted so as to ensure stable operation at resonance. However, this solution is difficult since the measurement of T requires space, extra weight and cost. This is why, instead of measuring T , it is better to measure Ψ[4] [START_REF] Budinger | Feeding and control electronic of a piezoelectric actuator[END_REF]; stable operation of the TWUM above the resonant frequency will be ensured if ( 22) is verified. This is why [START_REF] Giraud | Analysis and phase control of a piezo electric travelling wave ultrasonic motor for haptic stick application[END_REF] advises us to control Ψ using an internal loop as shown in figure 15. This increases reliability, but at the expense of a time delay for the control of W due to the internal loop of Ψ. This is why a hybrid control is proposed, which benefits from both control strategies.

C. Hybrid control strategy

This control depicted in figure 16 uses two controllers running at the same time. One is using a frequency controller C W to control W with respect to the reference value of W ref .

The other uses the controller in the rotating reference frame C Ψ to control Ψ = π

2 . An additional strategy block switches from one controller to the other, depending on their output. Under normal operating conditions, the motor is far from its resonance. So ω W is larger than ω Ψ . The frequency of the motor is then controlled (point (1) in figure 17). If T increases, Ψ decreases because V q increases, as shown in section III. At point [START_REF] Giraud | Precise position control of a travelling wave ultrasonic motor[END_REF], ω W = ω Ψ because W = W ref and Ψ = π 2 . For a larger load torque, W ref cannot be obtained; the frequency controller fails to control W ; in this mode ω W < ω Ψ , so Ψ is controlled. Consequently, W automatically adapts to the load torque: performances are diminished, but the motor does not stall.

D. Experimental runs

The hybrid control was checked during experimental runs. The TWUM used is a Shinsei USR60, driven by a linear power amplifier. A specific electronic control was built and used for controlling the motor in its rotating reference frame; details can be found in [START_REF] Giraud | Analysis and phase control of a piezo electric travelling wave ultrasonic motor for haptic stick application[END_REF]. An inertial load (10 -4 kgm 2 ) is attached to the shaft along with a powder brake (figure 18). Moreover, the position of the motor is controlled as detailed in [START_REF] Giraud | Position control of a small travelling wave ultrasonic motor[END_REF] [START_REF] Giraud | Precise position control of a travelling wave ultrasonic motor[END_REF]. A response time of 200ms and a position error of 1, 3mrad -which is equivalent to 5000 dots/rev -are required. During the experiment, Ψ ref is not equal to π 2 but 120 o as shown in section III. We use a Dspace 1104 board to implement the controller. During the first run depicted in figure 19, the powder brake is off. Only load inertia is applying torque to the shaft. This is the normal operating mode: the measured value of Ψ is larger than the reference, the frequency of the motor is fully controlled. For the second run of figure 20, an additional Fig. 19. Position step response when powder brake off @150V. load torque is applied using the powder brake. Since this additional torque is too large, we should ensure that W is smaller than its reference value in order to maintain stability. This is automatically achieved by using the proposed method. In fact, the motor operates close to its resonance, without stalling, and we thus obtain the best possible performances for this load condition. Of course, the response time is increased, but reliability is ensured.

V. CONCLUSION

In this paper, we present a novel control method which increases reliability of servo mechanisms using TWUM. This method is a hybrid control between a classical frequency control and a control in a rotating reference frame. Experimental results confirm performances at low supply voltage levels. Moreover, the strategy adopted is included in a position control loop.

Fig. 1 .

 1 Fig. 1. Vibrating stator and the rotor pressed on it.

Fig. 2 .

 2 Fig. 2. Instant values of the stationary waves and voltages and the equivalent complex phasors w and v.

Fig. 3 .

 3 Fig. 3.Vector representation of steady-state operations in the rotating reference frame; T = 0.

Fig. 4 .

 4 Fig. 4. Values of the motor's supply voltages and the output produced by the sensor for W = 2µm.

Fig.

  Fig. Complex phasors for W = 2µm.

Fig. 6 .

 6 Fig. 6. Values of motor's supply voltages and the output produced by the sensor for W = 3µm.

Fig. 7 .

 7 Fig. 7. Complex phasors for W = 3µm.

Fig. 8 .

 8 Fig. 8. Vector representation at resonance with two load torque values.

Fig. 9 .

 9 Fig. 9. Values for motor's supply voltages and the output produced by the sensor for W = 2µm and a load torque T = 0, 5N m.

Fig. 10 .

 10 Fig. 10. Complex phasors for W = 2µm and T = 0, 5N m.

Fig. 11 .

 11 Fig. 11. Safe Operating Area for TWUM in the (T ,W ) plane forV = V 1 and V = V 2, V 1 > V 2.

Fig. 13 .

 13 Fig. 13. Value of Ψ when the motor stalls for V = 150V and V = 110V .

Fig. 14 .

 14 Fig. 14. Frequency control of a TWUM.

Fig. 15 .

 15 Fig. 15. Control in a rotating reference frame.

Fig. 16 .

 16 Fig. 16. Proposed strategy.

Fig. 17

 17 Fig. 17. W as a function of T with hybrid control.

Fig. 18 .

 18 Fig. 18. The experimental test bench.

Frédéric

  Giraud was born in France in 1973. He is graduated from Ecole Normale Supérieure de Cachan, France in 1996 and received the B.S. degree in electrical engineering from Paris-XI University, Orsay, France in 1995, the M.S. Degree in electrical engineering in 1997 from Institut National Polytechnique de Toulouse, France, and the phD degree from Lille University, in 2002. He's a member of the electrical engineering and power electronics laboratory of Lille (L2EP) as an Associate Professeur. His research deals with the modelling and the control of standing wave and travelling wave piezo-electric actuators, for positioning and force feedback applications. Betty Semail was born in France, in 1964. She received her Ph. D degree in 1990 from University of Paris XI, Orsay and habilitation degree in 1997 from University of science and technologies of Lille. She is now professor in university of Lille 1. She is member of the electrical engineering and power electronics laboratory of Lille (L2EP) and responsible of the research axis upon control of electrical systems. She has been working motors and her main field of interest now deals with the modelling and control of piezo-electric actuators, for positioning and force feedback applications. Julien Aragones was born in France in 1975, is graduated of University of Montpellier where he achived his PhD in automation and signal processing in 2002. He joined Sagem défense sécurité in 2003 where he works on helicopter sight heads. He is in charge of turret control for advanced studies in airborne optronic. Jacques P. Robineau is gratuated of "Ecole Supérieure d'Electricité" (french top engineering school) in 1990 and achieved his PhD in electromagnetism. He has been working in electromagnetism numerical simulation (RF and optics) for 10 years, dealing with many aircraft and optronic laser stealth projects. He joined Sagem défense séecurité where he is in charge of advanced studies in airborne optronics. Jean-Thierry Audren was born in Paris, France on March 1st , 1951. He was graduated Engineer from Ecole Nationale Supérieure d'Electronique et d'Electromécanique de Caen in 1974. Since 1975 he has been in charge of RD in guidance, avionics and optronics systems. Since 1995 he has been leading research in the field of high power piezo motors in SAGEM. He is currently Deputy Director of the Air Land Warfare Business Unit of SAGEM, in charge of business Development.