Lyes Tighzert
email: ltighzert@gmail.com

Cyril Fonlupt
email: cyril.fonlupt@univ-littoral.fr

HAL is

from evolution theory. In 1999, Corno et al have proposed a new algorithm called selfish gene algorithm (SGA) [START_REF] Corno | The selfish gene algorithm: a new evolutionary optimization strategy[END_REF]. This later is inspired from the modern theory of evolution, called the selfish gene theory (SGT), proposed by the biologist Richard Dawkins [START_REF] Dawkins | The selfish gene[END_REF]. In SGT, the evolution is based on genes rather than individuals. Thus, the natural selection acts on genes rather than individuals. In the point of view of SGT, the breeds' bodies is made by genes only to assure their own replication through the reproduction process. Hence, Species are just vehicles that reproduce genes [START_REF] Dawkins | The selfish gene[END_REF]. Unlike other evolutionary algorithms that are based on a populations of individuals, the SGA proposed by Corno et al. uses the same paradigms of Dawkins's theory of evolution. Therefore, the SGA is based on genes and the population is replaced by the concept of virtual population (VP) that consists of a pool of genes. At each generation, two individuals, denoted G1 and G2 respectively, are sampled from VP. The SGA associates for each gene a selection probability that measure the chance that it have to be selected. The selection probability is tuned by a tournament selection, i.e. a completion between G1 and G2. The winner's gene are rewarded by increasing their selection probabilities by 1/N and the loser genes are punished by decreasing their selection probability by 1/N. Note that N is the number of gene competing for each locus. One can see that in SGA, the genes are competing for loci. The pseudocode of SGA is shown in Fig. 1.

___________________________________________

Algorithm1: Selfish Gene Algorithm [START_REF] Corno | The selfish gene algorithm: a new evolutionary optimization strategy[END_REF] ___________________________________________ 1. VP = uniform initialization of the virtual population 2. P = initialize all probabilities pij to l/ni 3. B = select (VP) / * assume the best */

Intelligent Trajectory Planning and Control of a Humanoid Robot

Using a New Elitism-Based Selfish Gene Algorithm Lyes Tighzert, Thafsouth Aguercif, Cyril Fonlupt, and B. Mendil L. Tighzert and B. Mendil have proposed recently a new variants called the replacing and never penalizing selfish gene algorithm (RNPSGA) that replace the loser's genes by those of the winner rather than penalizing them [START_REF] Tighzert | Realization of Gymnastic Movements on the Bar by Humanoid Robot Using a New Selfish Gene Algorithm[END_REF]. The authors have demonstrated its performance under several benchmark functions and have applied it to the realization of gymnastic movements of a humanoid robot. In this paper, we propose the reconciliation of SGA with elitism strategies. Permanent elitism-based selfish gene algorithm (peSGA) and nonpermanent elitism-based selfish gene algorithm (neSGA) are proposed and analyzed. The proposed elitisms strategies make selection pressure that is good enough to allow global convergence. Thus, the SGA performance are improved. The rest of this paper is organized as follow.

In next section, we present the proposed algorithms, i.e. peSGA and neSGA. Third section presents the validation of the proposed algorithms on the IEEE CEC2014 benchmark functions. In section four, we give an application of the proposed algorithms to optimization of bipedal waking of seven links humanoid robot. Section 5 concludes this paper.

II. ELITISM-BASED SELFISH GENE ALGORITHM

In order to enhance the performance of SGA, we propose two elitism strategies. In this subsection, we propose permanent elitism selfish gene algorithm (peSGA). This strategy is inspired from the works of Ahn Chang Wook and Rudrapatna S. Ramakrishna [START_REF] Ahn | Elitism-based compact genetic algorithms[END_REF]. The idea is to apply a selection pressure to accelerate the convergence while avoiding local convergence. In peSGA, at each generation only one individual, denoted G, is generated from VP. Then, a tournament is done between the elite (best) individual and G. The winner's gene are reworded and the loser's gene are penalized. Furthermore, the sampling mechanism is done differently than in the classical SGA. An adjustment procedure is added. During the sampling proves, each generated gene can be adjusted with a probability Pa. The adjustment formula is given in [START_REF] Plestan | Stable walking of a 7-DOF biped robot[END_REF].

A. Permanent elitism

Gi K+1 = Gi K +AR K+1 × randn(-1,1) (1)
Where Gi k+1 is the value of gene of locus I at the iteration K+1, Gi k is the value of gene of locus I at the iteration K, AR K is the adjustment rayon at iteration K.

The Adjustment rayon is decreased using formula [START_REF] Hong | Evolutionary Optimization for Optimal Hopping of Humanoid Robots[END_REF].

AR K = ARmax+ ( ARmax-ARmin)*K/Maxiter. ( 2
)
Where ARmax is the maximum value AR, ARmin is the minimum value of AR, K is the currant iteration index and Maxiter in the maximum number of iteration. The pseudocode of the proposed peSGA is given in Fig. 2. GFit=evaluate(G) , K=K+1 16.

B. Nonpermanent elitism-based SGA ________________________________________________

Gm=adjust G's genes with probability Pa using (1) 17.

Reduce rayon adjustment according to (2) 18.

GmFit=evaluate(Gm) , K=K+1 One of the drawbacks of peSGA may that the applied selection pressure made by keeping the elite as long as possible. This high selection pressure may result in premature convergence. In nonpermanent elitism based selfish gene algorithm (neSGA), the elite is guarded for a given number of cycles in which is not surpassed. This means that the elite is updated after a predefined inherence length denoted η. For the rest, peSGA is similar to neSGA. The pseudocode of neSGA is given in Fig. 3. We note that the VP and the probability selection in the proposed peSGA and neSGA are matrix of size (N×M), where N in the number of genes competing for each locus and M is the dimension of the optimization problem. Hence, each gene contains one allele and represents an optimization parameter. The main difference between peSGA and neSGA consists of the fact that peSGA finds a near optimal solution (i.e., a winner) that is maintained as long as other solutions generated from the virtual population (the pool of genes) are not better. In contrast, neSGA further improves the performance of the neSGA by avoiding strong elitism that may lead to premature convergence. The nonpermanent SGA comes with all the benefits of the peSGA. In addition, it maintains genetic diversity as a bonus. In the next section we evaluate the performance of peSGA and neSGA under IEEE CEC2014 functions.

III.

In order to evaluate the performance of the proposed algorithms, we propose a benchmarking under the IEEE CEC2014 benchmark functions [START_REF] Liang | Problem definitions and evaluation criteria for the CEC 2014 special session and competition on single objective real-parameter numerical optimization[END_REF]. This black box contains 30 functions including multimodal, unimodal, separable and non-separable functions. The dimension is fixed to D=30, the maximum number of fitness evaluation is set to 300000 as recommended in [START_REF] Liang | Problem definitions and evaluation criteria for the CEC 2014 special session and competition on single objective real-parameter numerical optimization[END_REF]. Each test is repeated 25 times. The averaged errors and the standard deviations are reported in Table I.

The proposed peSGA and neSGA are compared to the original SGA [START_REF] Ariff | Selfish Gene Algorithm Vs Genetic Algorithm: A Review[END_REF][START_REF] Corno | The selfish gene algorithm: a new evolutionary optimization strategy[END_REF] and to the most recent SGA variant known as the replace and never penalize selfish gene algorithm (RNPSGA) [START_REF] Tighzert | Realization of Gymnastic Movements on the Bar by Humanoid Robot Using a New Selfish Gene Algorithm[END_REF]. The parameters setting of the compared algorithms are as follow.

ARmax=5; ARmin=10 -10 ; Pov=0.15; Pa=0.01; N=50; η=N/2.

The mean and standard deviation of the function error value (f (x) -f (x*)) were calculated over 25 independent runs for each test function, where x is the best solution returned by the algorithm after 300000 fitness evaluations and x* is the global optimal solution. As the compared algorithms are stochastic and their results differ freom one run to another, a nonparametric test should be used to compare the results. Hence, Wilcoxon's rank sum test, recommended by [START_REF] Derrac | A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms[END_REF], at a 0.05 significance level was performed to test the statistical significance of the experimental results between two algorithms. The obtained results are reported in Table II. For each test function the proposed peSGA and neSGA are compared both to peSGA and PSGA. "-", "+", and "≈" denote that the performance of the corresponding proposed algorithm, e.g. peSGA or neSGA, is worse than, better than, and similar to that of the corresponding algorithm, respectively.

In order to rank the compared algorithm, the statistical Freedman test was performed. This test gives a global ranking of the compared algorithms. The results are presented in Table III. From this experimental study, we can conclude that the proposed algorithms are very competitive. They present performance enhancement for both unimodal and multimodal problems. According to Freedman test presented in Table 3, the proposed neSGA is the best algorithm for this experimentation. Thus, the proposed peSGA and neSGA can be used for real word applications.

+ + = + + = + - = - = + - - + + = + + = = - + = = - = - = + - = = = + = = = = = = = = - - = - = = = - = - + = = = = - = + + = + + = + - - - - + - = + + = = + = = - + = - - = - - + - = = = - - = = = = - = = = - = - = = = - = - + = = = = - - + 12 2 11 1 - 7 7 10 10 =
11 21 9 19 "-", "+", and "=" denote that the performance of the corresponding proposed algorithm is worse than, better than, and similar to that of the corresponding RNPSGA and SGA, respectively.

IV. APPLICATION TO HUMANOID ROBOT

A. Robot Model

The humanoid robot considered in this study is the HYDROID robot [START_REF] Ting | Contribution à la commande de robots marcheurs bipèdes[END_REF]. The dynamic parameters of the robot are given in Table IV. The length of each robot link except the feet is denoted li; the distance between the each link mass center is denoted lci; each link mass is denoted mi, the inertia moment of each link is denoted Ii. The robot foot structure is represented in Fig. 5. The vector F1 represent the ground reaction, xZMP represent the position of the zero moment point (ZMP) [START_REF] Vukobratović | Zero-moment point-thirty five years of its life[END_REF]. The hydroid robot is represented in Fig. 5. We consider here the model in the sagittal plane. This seven links robot has 9 degrees of freedom (9 DOF) which are left foot angle (α1), left knee angle (α2), left leg angle (α3), body angle (α4), right foot angle (α5), right knee angle (α6), right leg angle (α7), x body position and the z body position. The geometric, cinematic and the dynamic model of this robot are not reported here in details because of the limited paper pages. HYDROID is described in details in the literature (the lectors can see also [START_REF] Mousavi | Mathematical simulation of a seven link biped robot on various surfaces and ZMP considerations[END_REF]).For simulation, we have created a virtual model of HYDROID robot in SimMechanics on Simscape/Matlab2014. Fig. 5 gives the obtained 3D model. It shows the robot links, their center of mass and their associated frames. The ground contact is also modeled using an approximate model known as springmass-dumper model. Other methods can be found in the literature [START_REF] Jung | Dynamically adjustable foot-ground contact model to estimate ground reaction force during walking and running[END_REF].

The relationship between the foot penetration, foot velocity and the ground reaction force (F1) is given in (3)

z C z K F      (3)
Where F represent the reaction force, K and C are constants describing the nature of impact.

is the hypothetic penetration of the foot in the ground.

B. Dynamic Robot Model

The dynamic model can be obtained by Newton-Euler formalism or by Euler-Lagrange equation [START_REF] Braun | Actuated dynamic walking in a seven-link biped robot[END_REF]. In this paper, the Euler-Lagrange equation is used [START_REF] Tighzert | Realization of Gymnastic Movements on the Bar by Humanoid Robot Using a New Selfish Gene Algorithm[END_REF].

                       i i i q L q L dt d  (4)
Where i  Corresponds to the force or torque resulting from the dynamics of the i th link. L in the langrangian of the i th link. The obtained dynamic model if given in (5).

2 2 1 1 ) ( ) , ( ) ( F A F A B X Q X X X H X X D           
(5) Where X is the state variable vector composed from the articulations angles and the x and z body position X= [α1, α2, α3, α4, α5, α6, α7, x, y] T ; D(X) is the (9×9) inertia matrix;

H(X, X  ) is the centrifuge and Coriolis matrix of dimension (9×1); Q(X) is the gravitational matrix of dimension (9×1);

B is a ponderation matrix;  is the joint actuators torques at the articulations; F1 is the ground reaction on the left foot; F2 is the ground reaction on the right foot; A1 is the Jacobian matrix of the left foot contact with the ground; and A2 is the Jacobian matrix of the right foot contact with the ground. The dynamic walking can be considered as a hybrid process constituted of 3 phases. The double support phase (two feet on the ground), the simple support phase (one foot on the ground) and the impact phase (when the flying foot bumps the ground surface). The model given in (5) describe the double support phase. In the case of simple support phase, F1 or F2 is equal to zero.

C. The Impact Model

We have chosen to consider impulse impacts. This choice simplify the resolve of contact between the ground and the foot. The detailed impact model is given in (6). )

)( ( F A F A X X X D        (6)
Where  X  and  X  are the velocities after and before the impact, respectively. The constraints associated to impact are given in [START_REF] Goodrich | Foundations and trends in human-computer interaction[END_REF] and [START_REF] Calinon | Incremental learning of gestures by imitation in a humanoid robot[END_REF]. These latter gives the constraints of non-slip and the non-detachment of the feet.

0 ) ( 1   X X A  (7) 0 ) ( 2   X X A  (8)

D. Stability and ZMP

The stability of humanoid robot is one of the most topics studied. The stability of these articulated machines is assured if and only if the zero moment point (ZMP) is in the sustention polygon. This principle was proposed by Miomir Vukobratović [START_REF] Vukobratović | Zero-moment point-thirty five years of its life[END_REF]. In our case, the ZMP position is given in [START_REF] Eiben | Introduction to evolutionary computing[END_REF]. When this position is in the sustention polygon the robot is stable.

) ( 1 1 0 1 1 x p x z ZMP F h g m s F x     (9)

E. Trajectory generation

The trajectories of each articulation angle is assumed to be spline cubic function. This type of function has been used in several works. The walking period is denoted T. We decompose each robot step into two phases. The mathematical definition of each spline cubic function during each robot forward step is given in [START_REF] Corno | The selfish gene algorithm: a new evolutionary optimization strategy[END_REF] and [START_REF] Dawkins | The selfish gene[END_REF].

k k i k i t a t    3 1 ) (  , 2 / 0 T t   (10) k k i k i t b t ) 1 ( ) ( 3 1      , T t T   2 / (11)
Additional constraints can be imposed to assure on the continuity of the movement, i.e. the velocity and the acceleration at t=T/2. The parameters of the spline functions are calculated depending on the robot kinematic, dynamic and technologic constraints of the movement and an energetic criterion.

F. Kinematic and dynamic constraints

Robot walking consists on constrained movements. During footsteps, the stability must be assured. This leads to two constraints given in ( 12) and [START_REF] Liang | Problem definitions and evaluation criteria for the CEC 2014 special session and competition on single objective real-parameter numerical optimization[END_REF]. During the movements, the position in the z-axis of toe and heel of the mobile foot must be positive which leads to ( 14) and [START_REF] Derrac | A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms[END_REF]. Anthropomorphic constraints on the knees' angles are given in ( 16) and ( 17). The bottom foot, i.e. toe and heel, z-axis velocity must be positive to assure that the fixed foot take off the ground, i.e. obtain a transition to the next footstep. This leads to constraints [START_REF] Mousavi | Mathematical simulation of a seven link biped robot on various surfaces and ZMP considerations[END_REF] and [START_REF] Jung | Dynamically adjustable foot-ground contact model to estimate ground reaction force during walking and running[END_REF]. The ground reactions are also positive which leads to [START_REF] Braun | Actuated dynamic walking in a seven-link biped robot[END_REF] and (21). Where P is an identity matrix. T is the footstep period.

G. Evolutionary optimization

In order to realize optimal stable walking that minimize (22) and respect the constraints presented above, evolutionary optimization of the walking can be formulated. According to the results exposed in previous section, neSGA is the most competitive algorithm compared to RNPSA, SGA and peSGA. The evolutionary optimization algorithm is summarized in Fig. 6.

Fig. 1 .

 1 Fig. 1. The pseudo code of SGA.

Fig. 2 .

 2 Fig. 2. The pseudo code of peSGA.

Algorithm 3 :

 3 nonpermanent elitism-based SGA (neSGA) ____________________________________________________________ 1. M= optimization problem dimension, K=0; τ=0; 2. VP = uniform initialization of the virtual population 3. P = initialize all probabilities pij to l/N 4. N= initialize the competing gene number for each locus 5. POV= initialize the overwriting probability 6. Pa = initialize adjustment rate 7. η= N/2, / * the predefined inherence length */ 8. Best = select (VP) / * assume the best */ 9. BestFit=evaluate(Best) , K=K+1; 10. Elite = select (VP, P) / * assume an elite */ 11. ElitetFit=evaluate(Elite) , K=K+1; 12. while ( termination condition ) 13. if τ< η 14. G = select-individual ( VP, P) 15.

Fig. 3 .

 3 Fig. 3. The pseudo code of neSGA.

Fig. 4 .Fig. 5 .

 45 Fig.4. Robot foot structure and the hydroid robot[START_REF] Ting | Contribution à la commande de robots marcheurs bipèdes[END_REF]

-based SGA ________________________________________________ Algorithm2: permanent elitism-based SGA (peSGA) _____________________________________________________________

	1. M= optimization problem dimension, K=0;
	2. VP = uniform initialization of the virtual population
	3. P = initialize all probabilities pij to l/N
	4. N= initialize the competing gene number for each locus
	5. POV= initialize the overwriting probability
	6. Pa = initialize adjustment rate
	7. Best = select (VP, P) / * assume the best */
	8. BestFit=evaluate(Best) , K=K+1;
	9. while ( termination condition )
	10.	G = select-individual ( VP, P)
	11.	GFit=evaluate(G) , K=K+1
	12.	Gm=adjust G's genes with probability Pa using (1)
	13.	Reduce rayon adjustment according to (2)
	14.	GmFit=evaluate(Gm) , K=K+1
	15.	if GmFit > GFit	/* ">" means better */
	16.	G=Gm, GFit=GmFit
	17.	End if
	18.	if (fitness(G) > fitness(Best)) /* ">" means better */
	19.	reward-alleles (G) ;
	20.	penalize-alleles (Best) ;
	21.	Best=G, BestFit=GFit;
	22.	Overwrite Best's genes by G's Genes with probability POV
	23.	else
	24.	reward-alleles (Best) ;
	25.	penalize-alleles (G) ;
	26.	end if
	27

. End while 28. return Best

TABLE I .

 I EXPERIMENTAL RESULTS OF SGA, RNPSGA, PESGA, AND NESGA OVER 25 INDEPENDENT RUNS ON 30 TEST FUNCTIONS OF 30 VARIABLES WITH 300,000 FES. "MEAN " AND "STD"INDICATE THE AVERAGE AND STANDARD DEVIATION OF THE FUNCTION ERROR VALUES OBTAINED IN 25 RUNS, RESPECTIVELY.

TABLE II .

 II WILCOXON'S RANK SUM TEST RESULTS AT A 0.05 SIGNIFICANCE LEVEL BETWEEN BOTH PESGA AND NESGA AND EACH OF SGA AND RNPSGA.

	Algorithm		peSGA	neSGA
	Compared to	SGA	RNPSGA	SGA	RNPSGA
	Function		Comparison Results
	1
	2
	3
	4
	5
	6
	7
	8
	9
	10
	11
	12
	13
	14
	15
	16
	17
	18
	19
	20
	21
	22
	23
	24
	25
	26
	27
	28
	29
	30

TABLE III .

 III ALGORITHMS' RANKING ACCORDINF TO FREEDMAN TEST.

	THE BEST RANKS ARE GIVEN IN BOLD STYLE
	Function	SGA	RNPSGA	peSGA	neSGA
	1	2.0400	3.0800	3.2400	1.6400
	2	2.6800	3.1600	3.1600	1.0000
	3	2.2000	2.4400	2.7200	2.6400
	4	2.9600	2.8800	2.5600	1.6000
	5	2.5200	2.8400	3.5600	1.0800
	6	2.2800	2.4800	2.3200	2.9200
	7	3.2000	2.9600	2.8400	1.0000
	8	1.5600	2.5200	2.8400	3.0800
	9	2.2000	2.2000	2.6400	2.9600
	10	1.0400	2.8400	2.9200	3.2000
	11	2.3200	2.4000	2.1200	3.1600
	12	2.8400	3.0400	3.1200	1.0000
	13	1.1600	2.8000	3.0000	3.0400
	14	2.4000	2.1600	2.6400	2.8000
	15	2.5600	2.8800	2.4000	2.1600
	16	2.8800	2.4400	2.7200	1.9600
	17	2.4800	2.0400	3.0800	2.4000
	18	2.7200	2.7600	2.5600	1.9600
	19	2.7200	2.8400	2.4800	1.9600
	20	2.4000	3.1600	2.2800	2.1600
	21	2.3200	2.1200	2.7600	2.8000
	22	1.9200	2.6000	2.4800	3.0000
	23	2.0000	3.6400	3.3600	1.0000
	24	2.6000	3.0800	2.3200	2.0000
	25	1.7600	2.9600	2.8000	2.4800
	26	2.0800	2.3200	2.8000	2.8000
	27	2.2800	2.6400	2.5200	2.5600
	28	1.3600	2.5600	2.8800	3.2000
	29	2.1600	3.0800	2.2800	2.4800
	30	2.3600	2.7200	3.2800	1.6400
	Rank	2.2667	2.7213	2.7560	2.2560

TABLE IV

 IV

			.	HYDROID ROBOT PARAMETERS
	Links	Length		Mass	Inertia moment	center of gravity
		[m]		[Kg]	[kg.m 2 ]	[m]
	Foot	Lp=0.207	0.678 0.001	sx=0.135
		lp=0.072				sy=0.321
		hp=0.064
	Tibia	l1=0.392	2.188 0.028	0.1685
	Shin	l2=0.392	5.025 0.068	0.1685
	Trunk l3=0.543	29.27 0.815	0.01921

STD 8.8391e+006 5.9489e+006 3.0021e+007 3.5475e+007 6.1453e+003 5.5397e+003 5.1538e+001 1.5675e+001 2.0093e+001 5.3231e-002 4.7183e+000 1.3987e+000 1.2749e+000 3.3043e-001 5.4338e+000 2.1343e+000 1.6187e+001 4.7513e+000 3.4264e+001 4.0662e+001 6.6327e+002 2.5548e+002 5.6671e-001 1.8193e-001 2.2917e-001 6.8267e-002 3.3928e-001 6.3867e-002 4.2601e+000 2.3018e+000 3.1788e+00 3.0139e-001 1.0812e+006 1.1697e+006 1.6340e+005 6.1673e+005 2.9898e+000 1.2102e+000 5.6601e+003 5.9977e+003 2.7184e+005 3.7984e+005 8.4679e+001 6.4120e+001 3.3494e+002 6.8668e+000 1.5606e+002 2.8575e+001 1.8789e+002 2.1699e+001 1.0140e+002 2.9798e+000 4.0278e+002 3.7082e+001 4.9379e+002 6.6058e+001 4.7815e+004 1.2608e+005 2.3169e+003 1.1451e+003 1.0998e+005 7.2818e+004 1.8399e+003 2.1664e+003 5.3290e+003 5.4175e+003 3.5803e+000 2.1576e+000 1.9201e+001 3.9164e+000 3.8869e+000 1.5154e+000 1.6336e-001 8.0134e-002 1.0825e+001 4.7535e+000 2.0098e+001 9.1924e+000 3.1268e+002 1.1674e+002 6.8359e+002 2.2289e+002 1.0878e-001 6.0947e-002 5.6372e-001 1.7779e-001 4.4715e-001 2.1789e-001 1.5945e+000 1.0458e+000 2.9497e+000 5.1414e-001 3.9386e+004 4.3241e+004 6.2094e+003 6.2300e+003 1.6216e+000 9.4023e-001 3.7832e+003 4.6321e+003 1.7304e+004 1.0660e+004 1.7066e+002 5.0939e+001 3.2945e+002 7.6936e-005 1.4621e+002 2.4652e+001 2.02426e+02 4.06186e+00 1.1627e+002 3.6552e+001 4.0532e+002 9.2968e+001 6.8643e+002 1.4203e+002 2.1025e+002 9.4819e+000 1.1709e+003 2.7562e+002 1.6953e+006 1.4220e+006 2.6658e+003 2.8912e+003 5.6906e+003 4.6509e+003 1.5174e+001 1.8683e+001 1.9999e+001 1.6783e-003 4.5043e+000 1.4108e+000 1.6520e-001 1.0558e-001 1.0546e+001 4.3505e+000 1.9248e+001 9.1792e+000 2.8065e+002 1.0261e+002 6.6865e+002 2.6422e+002 1.1815e-001 8.6073e-002 4.5187e-001 1.6304e-001 5.6759e-001 2.8138e-001 3.1908e+000 1.9813e+000 2.8185e+000 4.6450e-001 6.1347e+005 5.1031e+005 6.2715e+003 7.0412e+003 1.6717e+000 1.0542e+000 3.1857e+003 5.8837e+003 5.6421e+005 7.3073e+005 1.8277e+002 1.0888e+002 3.2945e+002 6.4674e-005 1.4306e+002 2.7629e+001 2.0266e+002 1.2439e+000 1.0833e+002 2.7038e+001 4.0654e+002 3.9069e+001 6.7349e+002 1.2378e+002 7.8504e+005 1.0499e+006 1.3910e+003 3.7766e+002 2.0293e+006 1.4482e+006 4.1337e+003 3.4341e+003 7.4229e+003 4.8457e+003 1.8634e+001 2.1433e+001 2.0000e+001 6.8063e-004 5.1598e+000 1.5977e+000 1.3037e-001 6.8026e-002 1.2178e+001 4.7615e+000 2.3839e+001 8.8700e+000 3.7919e+002 1.4329e+002 8.5289e+002 2.7343e+002 1.5126e-001 1.1456e-001 5.6223e-001 1.9712e-001 4.5673e-001 2.6867e-001 2.7588e+000 1.1566e+000 2.8998e+000 4.3076e-001 1.0057e+006 6.0329e+005 7.3055e+003 8.1281e+003 1.9974e+000 8.1092e-001 2.5409e+003 2.8518e+003 6.0261e+005 7.6800e+005 1.9536e+002 1.0623e+002 3.2945e+002 2.7065e-005 1.3960e+002 1.4609e+001 2.0305e+002 1.1586e+000 1.2428e+002 4.2572e+001 4.1369e+002 4.6694e+001 6.7465e+002 2.0022e+002 4.2333e+005 8.4798e+005 1.6133e+003 2.9799e+002

________________________________________________

H. Results and discussions

We present in Fig. 7 the obtained optimal trajectory using neSGA. A proportional derivative control law is used to control articulations' trajectories. Despite the high complexity of this task, the humanoid robots are hybrid highly nonlinear under actuated machine with nine DOF, a satisfying results are obtained. neSGA assures stability and minimal energy consumption.

Hence, the proposed algorithms including neSGA, peSGA can be used in realworld applications. In the future works we can focus on the parameters adaptation of neSGA and peSGA.


In this paper, we propose a new variant of fish school search algorithm, called rate-learning based fish school search algorithm (RL-FSSA), that uses a new refinement strategy to guide the population of fishes towards the best solutions. The fish motion is based on a collective sensorimotor behavior. Their learning ability determines the fish with the best position (i.e., leader). In each cycle, the best fish ever found is refined by a rate learning based process. This refinement allows a local search and decreases the effect of the hazard through time, and the whole algorithm goes from diversification towards intensification. In order to evaluate its performance, the proposed algorithm is tested under a set of benchmark functions. The numerical tests include unimodal and multimodal functions. The results show the high performance of RL-FSSA compared to the standard version FSSA. Furthermore, the computational time of the fish school seach algorithm is reduced. For the validation, the algorithm is used for the trajectory planning and control of the mobile robot in an environment containing obstacles.

I.

INTRODUCTION

Evolutionary algorithms are stochastic search and optimization algorithms [START_REF] Huang | An Improved Artificial Fish Swarm Algorithm based on Hybrid Behavior Selection[END_REF]. Their origin can be traced back to the middle of the twentieth century [START_REF] Filho | A Novel Search Algorithm based on Fish School Behavior[END_REF][START_REF] Kennedy | Particle swarm optimization[END_REF][START_REF] Yu1 | An Improved GAFSA with Adaptive Step Chaotic Search[END_REF]. They are inspired from the theory of evolution [START_REF] Ae. Eiben | Introduction to evolutionary computing[END_REF]. In the last two decades, researchers have proposed other algorithms inspired from the social and biological behaviors of species. In the literature we talk about swarm based intelligence [START_REF] Zhenga | A Winner Determination Algorithm for Combinatorial Auctions Based on Hybrid Artificial Fish Swarm Algorithm[END_REF][START_REF] Xing | Fish Inspired Algorithms[END_REF][START_REF] Zhong-Hui | Hybrid Algorithm Based on Artifiial Fish Swarm Algorithm and Tabu Search in Distribution Network Reconfiuration[END_REF]. The swarm bioinspired algorithms tend to imitate a given social behavior of some insects, weed, and birds. This includes particle swarm optimization (PSO), ant colony (AC), bat colony (BC), invasive weed optimization (IWO) and fish school behavior.

The fish school behavior has inspired the works of Carmelo J. A. Bastos Filho who had proposed the fish school search algorithm (FSSA) [START_REF] Bastos Filho | A novel search algorithm based on fish school behavior[END_REF][START_REF] Bastos Filho | Fish school search[END_REF]. It is known that this global search and optimization method has a good ability of global convergence, associated with its simplicity and easy implementation [START_REF] Bastos Filho | Fish school search[END_REF]. As a self-organizing system, it hasn't sensitivity to the initial guess of the solution.

To sum up the FSSA strategy, we summarize it through 5 steps.

T. Aguercif is with the Laboratoire de Technologie Industrielle et de l'Information (LTII), Faculté de Technologie, Université de Bejaia, 06000 Bejaia, Algérie (aguercifthafsouth@gmail.com ).

L. Tighzert is with Laboratoire de Technologie Industrielle et de l'Information (LTII), Faculté de Technologie, Université de Bejaia, 06000 Bejaia, Algérie (e-mail: ltighzert@gmail.com ).

1. Initialization of the population, the weight of the fish, the individual movement, the collective instinctive movement, and the collective volitive movement.

2. Individual movement: Each fish updates randomly his position. They add to each position a random number multiplied by a predetermined step called the step individual. It allows each fish to reach a neighborhood position. This is formulated as

ind i i S rn x N     ) 1 , 1 ( (1)
Where i N is neighborhood position. i x is a current position.
rn : randn(-1,1) is a normal rand vector of dimension D, D is the dimension of the problem at hand. The movement only occurs if the new position Ni has a better fitness than the current position xi. Fitness difference (Δf) and displacement (Δx) are evaluated according to:

) ( ) ( i i i x f N f f    (2) i i i x N x    (3)
3. Feeding operation: The fishes are attracted towards the food. The weights of population are updated as

f f t w t w i i i      max ) ( ) 1 ( (4)
4. Collective instinctive movement: The school of the fish is directed towards the same direction denoted m(t), with

         N i i i N i i f f x t m 1 1 ) ( (5)
All fishes update their positions as follows

) ( ) ( ) 1 ( t m t x t x i i    (6)
B. MENDIL is with the Laboratoire de Technologie Industrielle et de l'Information (LTII), Faculté de Technologie, Université de Bejaia, 06000 Bejaia, Algérie (e-mail: bmendil@yahoo.fr).

C. Fonlupt is with Laboratoire d'Informatique Signal et Image de la Côte d'Opale (LISIC) Université du Littoral -Côte d'Opale BP 719 62228 CALAIS Cedex -France (cyril.fonlupt@univ-littoral.fr )

Rate-Learning Based Fish School Search Algorithm for Global Optimization

Thafsouth Aguercif, Lyes Tighzert, Boubekeur Mendil, and Cyril Fonlupt 5. Collective volitive movement: This movement is to determine the final position of all fish school. Sv : The step-vol parameter defines a contraction or a dilation of the swarm towards or outwards the barycenter of the school

     i i N i i w t w x t b ) ( ) ( 1 (7) )) ( ), ( (
)) ( ) ( ( ) 1 , 0 ( ) ( ) 1 ( t b t x Dis t b t x rn Sv t x t x      (8)
Researchers have developed several variants for FSSA. Some have proposed hybridization of the FSSA and the PSO [START_REF] Zhou | Hybrid Artificial Fish School Algorithm Based on Mutation Operator for Solving Nonlinear Equations[END_REF] with application for solving non-linear equations. Others have proposed an enhancement based on genetic operator [START_REF] Bastos-Filho | An Enhanced Fish School Search Algorithm[END_REF]. The binary representation has been proposed by A. G. João et al in [START_REF] Sargo | Binary Fish School Search Applied to Feature Selection: Application to ICU Readmissions[END_REF] and the chaos artificial fish swarm [START_REF] Huang | A new ICA Method based on Chaos Artificial Fish Swarm Algorithm[END_REF], to cite just a few.

The performance of FSSA are discussed and analyzed in several works. It is well know that it is a very competitive one compared to other bio-inspired, swarm and population-based search algorithms. The FSSA has been applied in several domain including combinatorial optimization [START_REF] Cai | Artificial Fish School Algorithm Applied in a Combinatorial Optimization Problem[END_REF] and image processing for the multi-threshold image segmentation [START_REF] Jiang | Image Segmentation with Improved Artificial Fish Swarm Algorithm[END_REF].

In this paper, we propose a refinement process to guide the fish school towards better positions (solutions). The proposed approach adds a local search algorithm to the FSSA to enhance its performance. It can be viewed as a mimetic fish school search algorithm. To well present our work, the remaining of this paper is organized as fallows. The second section discusses and analysis the proposed algorithm. In order to evaluate the proposed algorithm, numerical experiments are performed and presented in the third section. The fourth section gives a real world validation in the field of mobile robots, where the proposed algorithm is used for optimal path planning while avoiding obstacles. The section five concludes this paper.

II. RATE LEARNING BASED FISH SCHOOL ALGORITHM

In this section, we introduce a new variant of FSSA. The proposed algorithm is called rate learning-based fish school search algorithm (RL-FSSA). The main idea behind RL-FSSA is to use a refinement strategy to guide the population towards best regions in the search domain. As the refinement imitates the local search algorithm, the RL-FSSA can be classified under the umbrella ground of mimetic algorithms.

The sensorimotor capability of fishes allows them to interact with their environment and between them. Each fish can move individually (autonomy) and depending on the entire swarm (collective movement). In this study, we associate for the fishes the ability of learning. This learning stage determines the moving of the best fish. It is designed as refinement procedure based on a rate-learning capacity.

The rate learning based refinement allows to RL-FSSA more diversification at the beginning and intensification over time. Like FSSA, the proposed algorithm use the same steps as well as individual movements, feeding phase, collective instinctive movement, and collective volitive movement. However, the RL-FSSA uses an extra refinement phase. In each cycle, the algorithm seeks for the position assigned to the best fish ever found and refines it by a rate learning process. Here, the proposed "rate learning" operation is inspired from incremental learning [START_REF] Jin | A population-based incremental learning method for constrained portfolio optimisation. InSymbolic and Numeric Algorithms for Scientific Computing (SYNASC)[END_REF] defined by a rate-memorizing size. The main difference between our approach and incremental learning model resides in the fact that the rate learning strategy of our approach is a self-adaptive parameter that decreases through iterations. During the refinement cycles, the best fish is oriented to a new position which is computed using the memorized current position and a randomized number. [START_REF] Bastos Filho | A novel search algorithm based on fish school behavior[END_REF] Where:

) 1 , 1 ( ) 1 ( 1        rn RL x RL x t best t best
1  t best
x is the new position assigned to the best fish, t best

x is the actual position of the best fish, RL in [0,1] is the rate learning value , and rn is a normal randomized number between [-1,1].

During the refinement cycle, the rate learning value, RL, is dumped (scaled) by a dumping parameter, chosen between 0.95 and 0.999, denoted by RLdump in the pseudo of RL-FSSA shown in Fig. 1 dump RL RL RL   [START_REF] Bastos Filho | Fish school search[END_REF] The Fig. 1 shows in more details the pseudo code of the proposed algorithm. In the initialization step, we initialize all the input parameters, the size of school fish, the maximum number of iterations, and the initial rate learning RL. The fish weighs are initialized to be one for the entire population (line 9 of the pseudo code). The main loop contains the individual movement, feeding step, the swimming collective movement and volitive movement. The latter can be attracted/repulsed towards/outwards the barycenter depending on the rate of success in the school's last movement. In order to enhance the performance of FSSA, a refinement phase is added. In each cycle, the refinement process is executed ten times. In the first iterations, the RL-FSSA refines the best fish towards new position without taking much care on the current position. As the rate learning is decreased (dumped), the effect of the hazard will decrease through iterations, and the algorithm goes from diversification towards intensification. The parameters of this algorithm, like step individual size and step volitive size are also self-adaptive, as it is shown by lines 13 and 14 of the pseudo-code of Fig. 1. In the next section, we discuss the performance of RL-FSSA.

III. BENCHMARKING OF THE RL-FSSA

In order to evaluate our algorithm, we take from the literature a set of numerical benchmark functions. This set includes unimodal and multimodal functions. The benchmark mathematical formulas, their interval and optimum positions are shown in TABLE I. The two first benchmarks are unimodal and the others are multimodal functions. The dimension D of those functions is fixed to be 30.

We compare the proposed algorithm to the standard FSSA using the benchmark cited above. Each test repeated 15 times. The parameter setting of the numerical experiments are described as fallow:

-nPop: number of fish in the school : 30 -The number of fitness evaluation is set to 300 000

        N i i i N i i f f x t m 1 1 * ) ( 25. ) ( ) ( ) 1 ( t m t x t x i i    26.
end for 27. /* The collective-volitive movement */ 28.

for i=1:N do

29.     i i N i i w t w x t b ) ( * ) ( 1 30. )) ( ), ( (
)) ( ) ( ( ) 1 , 0 ( ) ( ) 1 ( t b t x Dis t b t x rn Sv t x t x      31.
end for 32.

Best= seek for the best fish position ever found 33. /* refinement the best fish position */ 34.

while

BENCHMARK TEST FUNCTIONS AND FUNCTION PARAMETERS

Name

Equation Search space Optimum

Sphere function

x ∈[-100, 100]

[0 0 …0] D D=30 Rosenbrock function         1 1 2 2 2 1 5 ) 1 ( ) ( 100 ) ( D i i i i x x x x f x ∈[-30, 30] [11 …1] D D=30
Rastrigin's function

      1 1 2 6 ] 10 ) 2 cos( 10 [ ) ( D i i i x x x f  x ∈[-5.12,5.12] [0 0 …0] D D=30
Griewank's function

1 ) cos( 4000 1 ) ( 1 1 2 7        D i i D i i i x x x f x ∈[-600, 600] [0 0 …0] D D=30
1. Ackley function

      ) 1 2 . 0 exp ( 20 ) ( 1 2 9 D i i x D x f e x D i   2 )) 2 cos( 1 exp (  x ∈[-32,32] [0 0 …0] D D=30
-Initial individual step = 1 -Final individual step =0.1 -Initial weigh is set to be 1. Those common parameters are for both RL-FSSA and FSSA and they are taken from the original FSSA paper [START_REF] Bastos Filho | Fish school search[END_REF]. The additional RL-FSSA refinement procedure parameters are set as follows:

-Initial rate learning = 0.5 -Rate learning dump = 0.99 -The number of refinement per generation is 10. The results of the numerical experiments are given in TABLE II. It gives the averaged value of the fitness found, for each function, for 25 independent runs. It gives also the processing time and the standard deviation. The best results are given in bold style. From the averaged values, we can study the quality of convergence. The standard deviation gives the repartition of the results around the averaged values. It is the most indicator of stability. The time processing gives the rapidity of the compared algorithm. From Table II, one can see that the proposed algorithm surpasses the original FSSA algorithm in terms of convergence quality, stability and rapidity. In the next section, we give an application of the proposed algorithm in mobile robot domain.

   D i i x x f 1 2 1
) (

IV. OPTIMAL PATH LEARNING AND OBSTACLE AVOIDANCE

A. Problem and Objective Definition:

Let us consider a mobile robot in an environment containing multiple obstacles. The objective is how to find an optimal path to reach the goal position while avoiding obstacles. This is handled in this work, using the RL-FSSA. The aim is to find a sequence of intermediary points between the initial and the goal positions. The trajectory is computed based on spline functions. The algorithm, iteration after iteration, will tune just the positions of the intermediary points. Each obstacle is augmented in order to consider the robot as a point (See Fig2). Two different environment are considered here.

B. Results and discussions

Fig. 2 gives the obtained results for each environment. The obtained optimal trajectory for each case is given on the left of the figure while the evolution of the optimization process is shown on the right of the figure. We note that cost function is the length of the trajectory. Although the environment contains multiple obstacles, the RL-FSSA has provided an optimal solution to avoid them. In this paper, an improved version of the FSSA is introduced. The algorithm is called the rate learning based fish school search algorithm (RL-FSSA).This new variant makes use of a refinement strategy to enhance the performance of FSSA. The refinement strategy is based on a self-adaptive rate-learning parameter (RL). The fishes will learn to move to new positions while taking increasingly the influence of the current position. Therefore, the memory guidance increases through iterations. Hence, the algorithm goes from diversification towards intensification.

The proposed algorithm is tested under a set of benchmark functions. Our study is focalized on the convergence quality, the accuracy, the stability and the speediness of the proposed algorithm. The results show the high performance of RL-FSSA compared to the FSSA original version. For the validation, the proposed algorithm is used for the optimal path planning of mobile robot with obstacle avoidance. Although the environment contains multiple obstacles, the RL-FSSA provided optimal solutions to avoid them.

indS

  : is an individual step.

____________________________________________________________________________________________________________________________________Algorithm1:

  Rate Learning based Fish School Search Algorithm ____________________________________________________________________________________________________________________________________ 1. /* initialization */ 2. it=0 /* generation index */ 3. MaxIt = maximum number of generation



Fig. 2 .

 2 Fig.2. Results of the expirements, In the left : Optimal trajectyory obtained by RL-FSSA, in the right the coreponding cost of the trjectory (optimization process)

  Pseudo code of Rate Learning Fish School Search Algorithm.

	Fig.1 .
				j ˂ =10
	35.	NewBest		1 (		RL	)	*	Best		*	randn
	36.		if	fitness	(NewBest	)	>	fitness	(Best	)	do
	37.			Best 	NewBest
	38.	end if
	39.	RL		RL		RL
									dump
	40. Return the best solution;
	_______________________________________________________________________________________________________

TABLE I .

 I

TABLE II .

 II RESULTS OF THE EXPERIMENTS: MEANS, TIME PROCESSING, AND STANDARD DEVIATION OF RL-FSSA AND FSSA

	Functions	F1	F2	F3	F4	F5
	MEAN	0.0013	14.4554	9.7261	7.2946	5.8421
	RLFSSA
	TIME	5.5913	6.4131	6.3567	6.6422	6.7465
	STD.DEV	5.9907E-4	14.4542	13.5654	12.4800	11.5343
	MEAN	0.2940	47.8327	61.1636	45.8782	36.8707
	FSSA
	TIME	5.7001	6.5556	6.5207	6.8154	6.9250
	STD.DEV	0.0310	65.7377	57.1026	56.0933	53.3076


S Dufour

G Vinsard

J.-P Brancher

Lemta

Lateral forces on a magnet falling in a conducting pipe

Keywords: Eddy currents, Magnets, Velocity, Finite Element Methods MSC2010 -35Q61, 35A10, 35A20, 78M10 Paper type -Research paper

Purpose -The purpose of this paper is to investigate the lateral forces during the fall of a magnet in a conducting pipe, when the direction of magnetization of the magnet is fixed. If the direction of magnetization is not parallel to the axis of the pipe, lateral forces occur and a decentration of the magnet happens. Design/methodology/approach -The problem is studied numerically, with a Th 3D FE formulation well-suited for the problem. Computational results are compared with experimental results. Findings -The physical model is given and the main force coefficients analyzed. The lateral forces and the decentration phenomenon are studied as a function of the main parameters (thickness and radius of the pipe). Originality/value -The direction of magnetization is a key parameter to analyze the dynamics of a magnet motion inside a conducting pipe, when the radii of the pipe and the magnet are not so close. This analysis with a fixed direction of magnetization allows to quantify the lateral forces and the decentration, and is a first step to understand the complete motion which includes the rotation which can be linked to the decentration.

Introduction

The fall of a magnet in a conducting pipe is a classical academic illustration of forces generated by eddy currents (Levin et al. 2006[START_REF] Partovi | Electrodynamics of a magnet moving through a conducting pipe[END_REF][START_REF] Roy | Demonstration of Lenz's law: Analysis of a magnet falling through a conducting tube[END_REF]). The corresponding engineering application is the building of magnetic dampers [START_REF] Bae | Modeling and experiments on eddy current damping caused by a permanent magnet in a conductive tube[END_REF]. In these cases, the radii of the magnet and of the pipe are close.

Analytical studies of the steady-state solution have been made over the years. All of them consider that the magnet falls along the axis of symmetry with a direction of magnetization of the magnet parallel to its axis, since analytical solutions can be found due to the axisymmetry. Then the different parameters are analysed, such as the conductivity [START_REF] Iniguez | Study of the conductivity of a metallic tube by analysing the damped fall of a magnet[END_REF], Levin et al. 2006), the pipe thickness [START_REF] Donoso | Magnet fall inside a conductive pipe: motion and the role of the pipe wall thickness[END_REF], and the excentricity with the assumption of an off-axis longitudinal motion [START_REF] Donoso | Damped fall of magnets inside a conducting pipe[END_REF].

However when the radius of the pipe is slightly larger than the one of the magnet, it has been noticed by experiment, that the position of the magnet (and the direction of magnetization) oscillates around the axis of symmetry. After having checked many reasons, an explanation is that a dynamical equilibrium state occurs either when the starting direction of magnetization is not parallel to the axis of symmetry, or when the starting position is not on the axis of symmetry.

When the radius of the pipe is even larger, then the magnet rotates around itself at an almost constant speed and it oscillates around a mean position which is not the axis of symmetry. The fall of a magnet near a plane represents in fact the limit case of a pipe of infinite radius [START_REF] Davis | Drag force on a magnet moving near a thin conductor[END_REF].

A full dynamical model of the motion of the magnet (valid in all these cases) requires to consider five degrees of freedom (dof): six for a solid minus the one of the rotation of the magnet around its axis of symmetry, that is the three dof of the position of the center and the two dof of the direction of magnetization. It requires a huge amount of work to build the expression of the system. The generalized force, which corresponds to the orientation of each of these dof, has to be computed, and the crossed effects have to be taken into account.

A step in such a work is to fix the direction of magnetization (which is the case of magnetic dampers), and to investigate only the forces along the axis of symmetry (brake forces) and the ones perpendicular to this axis (lateral forces): this step is the aim of the paper. Moreover, the validity of the computation will be controlled by experiment.

As the motion of the magnet mainly occurs in the plane described by the axis of symmetry and the starting position of the center of the magnet, the direction of magnetization is only considered in this plane. Lateral forces occur when the direction of magnetization is neither parallel nor perpendicular to the axis (even if the magnet is centered on the axis of the pipe).

This can be understood qualitatively in terms of equivalent wires where the magnet is considered as the superposition of two coils, one with its normal directed along the axis and the other with orthogonal normal. In the motion along the axis the normal directed coil induces circular loop currents with a normal also directed along the axis. The value of the mutual inductance between these last currents (considered as coils) and the coil directed orthogonally depends on the position of the magnet with respect to the axis; and in particular its derivative, which varies as the lateral force exerted on the magnet when the currents are fixed, is not zero on the axis.

Computing the forces may be done by different ways, analytically or numerically. On the one hand, when the direction of magnetization is perpendicular to the axis of the pipe, full analytical solutions are quite elaborate to find because of the complex paths of the current lines (Fig. 3, current lines obtained numerically).

On the other hand, a finite element classical magnetodynamic model may be used. However, it is a bit difficult to handle the computations at any scale of magnet and pipe, to take into account the lateral motion, and the different directions of magnetization, as they are linked to the limitations of the FE package that is used (in particular the mesh generator).

Thus, a mixed approach, with suitable 3D finite element formulations, valid for low induced magnetic fields, has been searched, in order to consider only the pipe. These formulations have been used with Taylor expansions, to find the dependencies of the brake and lateral forces as functions of the speeds. It allows to find the coefficients of the Taylor expansion with successive computations, like with analytical methods.

The different parameters of the problem are for the pipe, its inner radius, thickness and conductivity, and for the magnet its radius and magnetization intensity. The dependence of the forces with respect to the conductivity and magnetization intensity is linear. The radius of the magnet will be fixed, and the pipe's inner radius and thickness will be the effective studied parameters, to consider the lateral and brake forces of a magnet falling in a pipe as a function of the decentration and direction of magnetization.

R 1 R 2 x † k x k z k y (Y • k z )k z R ϕ d Y H e Figure 1.

Geometry and notations

Position of the Problem

The axis of the pipe (interior and exterior radii R 1 and R 2 ) is along the z-direction (Figure 1). The position of the center of the spherical magnet (radius R) is in a fixed frame Y (t) and its direction of magnetization d.

The motion of the magnet generates magnetic flux variations and then electromotive forces that drive eddy currents in the pipe. For completeness it would be necessary to take into account both the longitudinal and lateral motions as source of eddy cur-rents; and then it would be impossible to find an expression of the force that only depends on the state variables of the motion.

Fortunately, the essential part of the motion is the longitudinal motion with a speed almost constant; then the lateral motion or the fluctuations around this constant speed are low and the corresponding eddy currents are negligible.

The source field of a spherical magnet [START_REF] Jackson | Classical Electrodynamics[END_REF] is in this fixed frame:

h s = M R 3 |y -Y (t)| 3 d • (y -Y (t)) |y -Y (t)| 2 (y -Y (t)) - d 3 (1)
where y and Y are respectively the considered position and the position of the center of the magnet, M is the magnetization intensity and d its direction in space. The position of the center of the magnet is and its speed are:

Y = x † k x + z k z and v = ẋ † k x + ż k z (2)
As the speed in the vertical direction is almost constant, it is interesting to consider a moving frame by displacing vertically the center of the fixed frame at speed of

Ẏ (t) • k z = -żk z (3)
In this moving frame the current point x is linked to y by:

x = y -(Y • k z ) k z (4)
and then the position of the magnet is X = x † k x . The magnetic source field of the magnet is (outside the magnet):

h s = M R 3 |x -X(t)| 3 d • (x -X(t)) |x -X(t)| 2 (x -X(t)) - d 3 (5)
d is assumed to be in the (x, z) plane as the motion itself:

d = sin ϕ k x + cos ϕ k z (6)
The source electric field (EMF) corresponding to the magnet motion is

e s = µ 0 v × h s = µ 0 ( ẋ † k x + ż k z ) × h s ( 7
)
where h is the magnetic field generated by the eddy currents density j, that is:

∇ × h = j ; ∇ • h = 0 (8)
the current density verifies:

∇ • j = 0 (9)
then it is linked to the EMF by the Ohm law:

j = σ(e s + e i + ∇φ) ( 10
)
where σ is the conductivity, φ the electrical scalar potential (needed to insure the free-divergence of the induced current density) and e i the EMF corresponding to the induced current motion:

e i = µ 0 ż k z × h (11)
There is no term in µ 0 ẋ † k x × h s in (11) since the induced current only moves in the z-direction. The rotational of (10) in the pipe leads to:

∇× j σ = µ 0 ∇× (v × h s ) + µ 0 ∇× ( żk z × h) (12)
And finally, the force exerted on the magnet is the opposite of the (Laplace) force exerted on the pipe:

f = -µ 0 D j × h s dx 3 (13)
where the induced magnetic field is not considered because the induced currents have no global effects on themselves. The steady-state assumption allows to express this force as a function of the state variables x † , ẋ † , z, ż of the magnet (in addition to the geometrical parameters R 1 and R 2 of the pipe and of its electrical conductivity).

And the next section explains the used method to obtain this result through numerical computations.

Numerical Model

We use the current function T in the conducting region D such as:

j = ∇× T (14)
(to satisfy the free-divergence of j in D), submitted to the conditions

∇ • T = 0 in D T × n = 0 on ∂D (15)
which insure j • n = 0 on ∂D and the uniqueness of T . The magnetic vector potential (MVP) is computed by using the Biot-Savart law:

a(x) = µ 0 4π D ∇× T (y) |x -y| dy 3 (16)
The induced magnetic field h solution of ( 8) is then obtained with:

h = 1/µ 0 ∇ × a (17)

Electromagnetic fields

In a certain sense, the use of both current function T and magnetic field h is a mixed formulation [START_REF] Dular | Magnetostatic and magnetodynamic mixed formulations compared with conventional formulations[END_REF]. It is in fact possible to compute the values of these variables in one step (strong coupling), with classical T -Ω or a -φ methods. The main advantage to consider separately the computation of T and h [START_REF] Dufour | Iterative computation of induced currents in a three-dimensional conducting domain[END_REF] is that only the conducting domain can be considered, whereas other methods need to mesh the domains outside the conducting domain, in particular the magnet (which raises numerical difficulties to get the appropriate forces when the magnet has a lateral motion).

Let us consider the solution U , which depends on a given vector field W , of the generic variational problem: ∀ T ′ verifying ( 15)

D (∇× U -σµ 0 W ) • ∇× T ′ dx 3 = 0 (18)
under the constraint:

∇ • U = 0 in D; (19)
This dependence of U with respect to W is summarized by the introduction of a shorthand notation (functional operator) G that associates the solution U of (18-19) to the considered vector field W :

U = G(W ) (20)
Then the current function T defined by ( 14) with ( 15) and ( 12) is

T = G(v × h s + ż k z × h) (21)
Nodal finite elements with a penalty method (Ribaki et al., 1990) to enforce the freedivergence ( 19) are used to compute T . The MVP a is a function of the current function T with the help of the Biot-Savart law (16). Then a is computed on the nodes of the tetrahedral mesh of D. The Gauss quadrature rule may be used but has the drawback that the integral contains singularities for the tetrahedron belonging to the node where a is computed. A treatment with locally-corrected quadrature rules for these tetrahedrons may be done [START_REF] Strain | Locally corrected multidimensional quadrature rules for singular functions[END_REF] but we have preferred to separate and treat analytically (with formulae issued from a symbolic computation) the singular parts of the integrals to keep simultaneously a low order of Gauss quadrature and enough accuracy. The next step is to compute the induced magnetic field: it is obtained as piecewise on the element vector field by the direct differentiation (17) of the linear approximation of the MVP in the element obtained by the Biot-Savart law (16). These operations are summarized once again by the introduction of the shorthand notation:

h(x) = F (T (x)) = 1 4π ∇ [x] × D ∇ [y] × T (y) |x -y| dy 3 (22)
The processes where T is computed for a given h with (21) (electric step) and h for a given T with ( 22) (magnetic step) may not converge. Such lack of convergence occurs if the magnetic field h play a major role in the EMF. Fortunately in the case where the speed of fall corresponds to a low magnetic Reynolds number σµ 0 R 1 | ż| (less than unity for speeds below 2 m/s), the convergence happens and moreover is quick.

From a practical point of view, the iterative finding of (T , h) which is done by ( 21)-( 22) is eased by a Taylor expansion of the vectors fields as a function of the speeds ( ẋ † , ż):

T = ∞ n=0 żn (T nx ẋ † + T nz ż) h = ∞ n=0 żn (h nx ẋ † + h nz ż) (23)
Then the zero-order fields verify:

T 0x = G(k x × h s ) and h 0x = F (T 0x ) T 0z = G(k z × h s ) and h 0z = F (T 0z ) (24)
They correspond to the case where the induced EMF is neglected. The next order terms verify: for j ∈ {x, z} and n ≥ 0

T (n+1)j = G(k z × h nj ) and h (n+1)j = F (T (n+1)j ) (25)

Force components

The iterative computation of T allows to find directly the force as a function of the speeds, since the force (obtained with ( 13)) is linear with respect to T and then is (formally) a series as:

f = ∞ n=0 żn (f nx ẋ † + f nz ż) (26)
where the coefficients depend on ϕ and x. The components of the Laplace force (26) are such as: for (i, j) ∈ {x, z} 2 and n ≥ 0

f ni • k j = D 1 σ ∇× T 0j • ∇× T ni dx 3 (27)
For n ≥ 1, these components can also be written as:

f ni • k j = µ 0 D (h (n-1)i × (∇× T 0j )) • k j dx 3 (28)
The advantage of (28) with respect to ( 27) is that only the (n -1)-order field components are needed for the computation of the n-order force components.

The Taylor expansion of the force as a function of the speeds is well suited to analyse the dependencies of each coefficient on ϕ and x † . Since the lateral ẋ † and longitudinal ż speeds are low, ( 26) is limited to the second order only as a function of ż and ẋ † : higher orders have been investigated and are negligible. The first order components are computed with ( 27) and the second order with (28) for n = 1: only the computation of the zero-order fields (T 0x , h 0x , T 0z , h 0z ) is necessary.

As the source field ( 5) is a sinusoidal function of ϕ these components depend sinusoidally on 2ϕ. The dependence on x † is more complicated, however as our objective is to determine the tendencies of the motion near the centered position x † = 0, we bound the dependence to the first order on x † . We obtain then

F x ≃ -ẋ † α x (1 -β x cos(2ϕ)) + γ sin(2ϕ) ż -κ x † ż2 (29)
F z ≃ -ż α z (1 + β z cos(2ϕ)) + γ sin(2ϕ) ẋ † -κ x † ẋ † ż (30)
where F x = f • k x and F z = f • k z are the lateral and axial forces and where the constants (α x , β x , α z , β z , γ, κ) can be computed by with (T 0i , h 0j ).

The terms -α z (1 + β z cos(2ϕ)) and -α x (1 -β x cos(2ϕ)) are exactly f 0z • k z and f 0x • k x , computed with (27); the α z/x , β z/x coefficients are identified by using two computations with two directions of magnetization: one in the fall (ϕ = 0) and the other in the lateral (ϕ = π/2) directions. We obtain 0 < β z/x < 1 (in accordance to the fact that the fall speed is lower for ϕ = 0 than ϕ = π/2 for the z coefficient). Since ∂ x h s | ϕ=0 = ∂ z h s | ϕ=π/2 , the fields T nx | ϕ=0 and T nz | ϕ=π/2 are the same. It leads to the relation

α z (1 -β z ) = α x (1 -β x ).
The term γ sin(2ϕ) is f 0z • k x or equivalently f 0x • k z due to the symmetry of ( 27); it is the source term of the decentration. It happens when the set direction of magnetization is neither parallel nor perpendicular to the axis of the tube, even if the magnet is centered in the pipe. To explain the presence of this term, we can consider that the main component of T 0x is its radial component T 0xr , then γ ≃ µ 0 D ∂ z T 0xr | π 4 h sr | π 4 dx 3 (h sr being the radial component of h s ) can be computed by:

γ ≃ µ 0 D (∂ z T 0xr | π 2 h sr | 0 + ∂ z T 0xr | 0 h sr | π 2 )dx 3 (31)
Figure 2 shows the different fields of ( 31 And last, the term -κ x † ż2 is f 1z •k x ; On the contrary to the previous coefficients which can be found with T 0(x/z) and ( 27), the computation of κ needs h 0(x/z) and ( 28). The term -κ x † ż2 is analogous to a restoring force that tends to center the magnet on the axis. The accuracy of the Th formulation has been verified by a comparison with a classical T -Ω formulation: the agreement is good. The power loss has been computed in both cases (with the same mesh for D) for several thickness (1 mm ≤ e ≤ 8 mm): the difference is less than 0.5%.

T 0xr | ϕ= π 2 T 0xr | ϕ=0 h sr | ϕ=0 h sr | ϕ= π 2 - + + - + - + - - - - + + ∂ z T 0xr | 0 h sr | π 2 ∂ z T 0xr | π 2 h sr | 0 ∂ z T 0xr | 0 h sr | π 2 + ∂ z T 0xr | π 2 h sr | 0

Centered Magnet

The magnet is supposed to be centered in this section, and the investigated parameters are the thickness and the inner radius of the pipe (e = 0.5 -8 mm, R 1 = 8 -16 mm). Figure 3 shows the current lines of T 0zr , the main component (along the radial direction) of T 0z for two directions of magnetization on one half of the tube. For a small thickness (e = 0.5 mm), T 0zr is almost constant as a function of r. If the direction of magnetization is in the fall direction (ϕ = 0), the field lines are axisymmetric. For ϕ = π/2, the cut plane is parallel to the direction of magnetization. T 0zr is an odd function of z and even function of y. For a large thickness (e = 8mm), T 0zr decreases as a function of r, and it is almost negligible for r > 3 mm. The thickness of the pipe is an important factor: the problem is three-dimensional for large thickness, whereas a model reduction in the r direction may be acceptable for lower thickness. In this case a shell current approximation may be used, where the induced current is assumed to be concentrated in a shell and the thickness is taken into account in an average treatment of (12). Figure 4 (left) shows the brake force coefficients of (30) α z (1 + β z cos(2ϕ)), for ϕ = 0 and ϕ = π/2, versus the thickness of the pipe, for a constant inner radius of the pipe. For small values, the coefficients are almost linear with respect to the thickness: the r-dependence of the currents is low. In this case the shell current approximation (small dots) may be used: it valid for a thickness smaller than 2 mm.

f 0z • k z (kg .
Above a thickness of 6 mm the brake force coefficients are almost constant: the induced currents vanish at the outer radius for these depths.

The brake force coefficients of (30) are investigated as a function of the direction of magnetization for several pipe radii R 1 (Fig. 4 (right)). The variation of the brake force between maximum and minimum as a function of the direction of magnetization is about one third.

The computation has also been compared to experiment. The measurement of the brake force for a given direction of magnetization has been investigated. The problem is to insure a fixed direction of magnetization. Instead of a fall, a system of rope pulley and counterweight has been built to lift the magnet, and control the speed. The frictional force of the pulley has been measured with the help of a non-magnetized sphere. It is then possible to get the brake force as a function of the direction of magnetization. These experimental values are displayed with the larger dots. The model predicts quite well the variation of the brake force as a function of the direction of magnetization for different inner radii of the pipe (8 -10 -16 mm).

f 0j • k x (kg .
Figure 5 (left) shows the lateral force coefficients of (29) α x (1 -β x cos(2ϕ)), for ϕ = 0 and ϕ = π/2. The same current shell approximation can be made for values of thickness lower than 2 mm. If the lateral speed is of the same order of magnitude as the vertical speed, then the main term of the lateral force is in this case f 0x • k x : it is of the same order of magnitude as the brake force (Fig. 5 (right)). The decentration force (term in γ) is about one tenth of the brake force for ϕ = π/4, which corresponds to its maximum value.

Decentration of the magnet

For a constant decentration, the lateral forces are the decentration term γ sin (2ϕ) ż balanced by the restoring force -κ ż2 x † .

This first-order approximation of the restoring force on x † was made in section 3.2 to obtain a simple expression of the restoring force near the axis: it comes from the Taylor expansion of h s (5) to the first order as a function of x † . The computations can also be made by using an effective decentration x † for h s (i.e. by entering the effective lateral position X in (5)).

The comparison between the first-order approximation and the less approximated model is shown (Fig. 5) for some values of internal radius of the pipe. The first-order approximation is valid for low decentrations (in particular for the pipe with R 1 = 8 mm). For higher decentrations, the restoring force becomes a cubic function, and a third-order approximation could be more accurate.

The linear approximation gives as computed lateral equilibrium position x † = 8 mm for R 1 = 12 mm (beyond the validity of the first-order approximation of the restoring force) and 6.3 mm for R 1 = 16 mm for ϕ = π/4. The experimental values are x † = 5.5 mm for R 1 = 12 mm (the magnet hits the pipe) and 7 mm for R 1 = 16 mm.

Conclusion

The force exerted on a magnet which falls in a conducting pipe, with a fixed direction of magnetization, depends on this direction of magnetization and has two components. The first component is directed along the symmetry axis of the pipe and tends to brake the fall. The lowest speed is found when the direction of magnetization is parallel to the symmetry axis and the highest when it is perpendicular: the relative difference between them is about one third.

The second component, labelled lateral, is perpendicular to the axis of symmetry, and tend to decentrate the magnet when the direction of magnetization is neither parallel nor perpendicular to the axis of symmetry, even if the term of decentration in the lateral force is less than 10% of the brake force. The decentration increases with the increase of the radius of the pipe with respect to the one of the magnet.

The chosen model of computation has allowed the quantification of these force components and their parametric analysis as a function of the speeds. It has been verified by experiment as far as possible.

The validation of this model when the direction of magnetization is fixed will allow to study the full motion of the magnet, i.e. its fall with decentration and rotation, with the same methodology. The current density and magnetic fields will be found iteratively, as a function of the speeds (rotating speeds included), and the analysis of the forces and torques with respect to the decentration and direction of magnetization will allow to find the ordinary differential set of equations that the motion of the magnet verifies.

  Figure2shows the different fields of (31) with their directions (in the inner radius of the pipe |x • k r | = R 1 ). Due to the symmetries of the fields (h sr | ϕ=0 and T 0xr | ϕ=π/2 are symmetric with respect to the plane z = 0, whereash sr | ϕ=π/2 and T 0xr | ϕ=0 are anti-symmetric), f 0x | ϕ=0 •k z , which is approximately equal to D ∂ z T 0xr | 0 h sr | 0 dx 3 , vanishes. For the same reasons f 0x | ϕ=π/2 • k z is also null. It is not the case for the term ∂ z T 0xr | π2 h sr | 0 : it has two vortices, centered on z = 0. The term ∂ z T 0xr | π 2 h sr | 0 has four vortices in the opposite direction, thus the integral of their sum γ is not zero.And last, the term -κ x † ż2 is f 1z •k x ; On the contrary to the previous coefficients which can be found with T 0(x/z) and (27), the computation of κ needs h 0(x/z) and (28). The term -κ x † ż2 is analogous to a restoring force that tends to center the magnet on the axis.

  Figure 2. Fields used for the computation of the decentration term

Figure 3 .

 3 Figure 3. Lines fields of T 0z • kr (main component of the current function)

Figure 4 .

 4 Figure 4. Brake force coefficients f 0z • kz as a function of the pipe thickness (left) and the direction of magnetization (right) without decentration

  Figure 5. Lateral force coefficients f 0j • kx as a function of the thickness (left) and the direction of magnetization (right) without decentration

  Figure 6. Restoring force coefficients f 1z • kx as a function of x † (e = 1 mm)


Gérard Vinsard

Esteban Saatdjian

Stéphane Dufour

No-Slip Motion of a Spherical Magnet on the Top of a Conductive Plate

Keywords: Conducting materials, eddy currents, finite element methods, magnets

I. INTRODUCTION

T HE coupling between the partial differential equations (PDEs) that describe the electromagnetic field and the ordinary differential equations (ODEs) that give the motion of a moving conducting object can be a particularly difficult problem especially when the control parameters of the ODEs use variables that issue from the field problems.

Most coupled mechanical-electromagnetic problems concern devices which have only one mechanical degree of freedom [START_REF] Weidermann | Lorentz Force and Joule Heat Induced in an Electrically Conducting Plate Moving With Time-Dependent Velocity Under the Influence of a Homogeneous Magnetic Field[END_REF]- [START_REF] Zec | Fast Technique for Lorentz Force Calculations in Non-destructive Testing Applications[END_REF]- [START_REF] Paudel | General 2-D Steady-State Force and Power Equations for a Traveling Time-Varying Magnetic Source Above a Conductive Plate[END_REF], notable exceptions are magnetohydrodynamical flows. Most cases describe either the field-circuit coupling [START_REF] Melgoza | Strong Coupling of Electromagnetic Transients and Finite Element Magnetic Field Solvers[END_REF] or the mechanical coupling [START_REF] Carpentier | Modeling of Magneto-Mechanical Coupling Using Magnetic Volume Integral and Mechanical Finite-Element Methods[END_REF].

If a moving spherical magnet is on top of a conductive plate, eddy currents are induced and tend to slow down the spherical magnet until it stops. This is a transient problem. The position and the magnetization direction are the five mechanical degrees of freedom. The control parameters of the ODEs are the Lorentz force and the torque, which are calculated by solving the transient EM field problem. Similarly, the electromagnetic field depends on the position of the moving magnet. The aim here is to find the dynamics of the magnet, and to validate the accuracy of results by experiments.

Although this is, at first glance, an academic problem, its solution strategy will be used in future work on granular flows. These flows may jam due to Coulomb's friction forces. If spherical magnets, of the same size as the grains, are introduced and submitted to forces like the ones studied here, the flow may be unblocked.

II. GOVERNING EQUATIONS

The spherical magnet (radius R) is assumed to roll without slipping on a horizontal plane. The position X = X(t) k x + Y (t) k y + R k z of its center of gravity is obviously time-Manuscript received June 24, 2015. Corresponding author: G. Vinsard (email: gerard.vinsard@univ-lorraine.fr). We benefited from the savoir faire of Raoul Lenoir in preparing this paper. dependent, as is the (planar) velocity

˙ X = Ẋ k x + Ẏ k y . The magnetization is directed along d(t) = -sin θ sin ψ k x + sin θ cos ψ k y + cos θ k z (1)
where ψ, θ are the precession and nutation Euler angles (Fig. 1). The intrinsic rotation ϕ around the axis of direction d has to be considered to describe the motion of the magnet but has no influence on the magnetic field.

Due to the no-slip condition, the instantaneous rotation vector ω is a function of the planar velocity components with only one extra kinematically free component ω z :

ω(t) = -Ẏ (t)/R k x + Ẋ(t)/R k y + ω z (t) k z
(2) induced by the motion of the magnet. These eddy-currents are determined by using a timedependent -Ω formulation. The magnetic field due to the magnet is given by the dipolar expression (outside the magnet):

h s = ∇ x (Ω s ) ; Ω s ( x) = - M d • -→ xX 4π | -→ xX| 3 ; -→ xX = x -X (4)
where M (A m 2 ) is the magnetic moment of the magnet. The time derivative of h s is ˙ h s = ∇ x ( Ωs ) where

Ωs = M 4π | -→ xX| 3 d • ˙ X -3 ˙ X • -→ xX | -→ xX| 2 -→ xX -ω × -→ xX (5)
Notice that this term, which is the time-derivative of a kind of magnetomotive force, is a function of all the variables of the dynamical system (3) except ϕ.

In the T -Ω formulation, the system of PDEs to be solved both in the conductive plate D and outside the plate E 3 -D is (D the conducting domain with a conductivity σ and a permeability µ 0 , ∂D its boundary and n the outward normal field):

∇× 1 σ ∇× T + µ 0 ∂ t T + ∇Ω + µ 0 ˙ h s = 0 in D (6) ∇ • T = 0 in D (7) T × n = 0 on ∂D ; T = 0 in E 3 -D (8) ∇ • µ 0 ∇Ω + T = 0 in E 3 (9)
The eddy current density and magnetic field are j = ∇× T and h = T + ∇Ω and the total Lorentz force and torque exerted on the plate are respectively (d x 3 = dx dy dz)

-f = D ∇× T × µ 0 h s d x 3 -Γ = D -→ xX × ∇× T × µ 0 h s d x 3 (10)
Due to Newton's third law, they are exactly equal and opposite to the force and torque f and Γ exerted on the magnet (labelled D s ), whose expressions are

f = ∂Ds -µ 0 M( d • ∇Ω) n d x 2 Γ = Ds µ 0 M( d × ∇Ω) d x 3 (11)
The system (6-9) has to be solved in time. A time step τ is given; T , Ω at discrete times n τ are denoted T n , Ω n and an estimation

˙ h * n s (described below) of ˙ h s at this time n τ is used. The full time-implicit weak form version of the T -Ω system (6-7-8) is: ∀ T ′ such that T ′ × n = 0 on ∂D, ∇ • T ′ = 0 on D D ∇× T n+1 σ • ∇× T ′ + µ 0 τ T n+1 + ∇Ω n+1 • T ′ d x 3 = D µ 0 τ T n + ∇Ω n • T ′ d x 3 - D µ 0 ˙ h * n+1 s • T ′ d x 3 (12)
Practically the constraint ∇ • T = 0 in D is taken into account by an extra penalty term. Since the only contribution of ∇ • T in E 3 is the product between the jump of T • n on ∂D and the surface Dirac distribution, the Biot and Savart form of (9) reduces to:

Ω n+1 ( x) = µ 0 4π ∂D T n+1 ( y) • n | x -y| d y 2 (13)
Even if Ω can be computed in E 3 , the computation is only required in the domain D to find T with (12). It allows us to use a fixed mesh (and then to disregard the motion of the magnet in the mesh) as well as a fixed frame of reference. For a given ˙ h s at the time (n + 1) τ , the computation of ( T n+1 , Ω n+1 ) can be performed from a given set of ( T n , Ω n ) by (12-13) with the finite element method. To avoid handling a non-sparse matrix, (12-13) is computed iteratively: T n+1 is computed with (12) for a given Ω n+1 , and Ω n+1 with (13) for a given T n+1 . Practically, if the time step is not too large, convergence is obtained with a limited number of iterations, which depends on the used convergence criterium. The computation of (13) is done once and the corresponding matrix is stored. The singularities contained by (13) are handled with closed-form expressions.

If ˙ h * n+1 s really corresponds to ˙ h n+1 s
, the value of ˙ h s at the time (n+1) τ , the components of the dynamical system (3) have to be computed implicitly. It is possible to solve such a system implicitly, but it requires an iterative process to handle the non-linearities. An estimation of ( X n+1 , d n+1 , ω n+1 ) is first required to compute ˙ h n+1 s , then T n+1 is obtained. The forces and torques are deduced and an implicit scheme of (3) is used to find a new estimate of ( X n+1 , d n+1 , ω n+1 ). Such a method of solution would be time-consuming due to the numerous degrees of freedom (including the T -Ω ones), and the fact that the problem is not time-periodic, even if diagonally implicit Runge-Kutta (DIRK) methods were to be used [START_REF] Nicolet | Implicit Runge-Kutta Methods for Transient Magnetic Field Computation[END_REF]- [START_REF] Pries | Steady-State Algorithms for Nonlinear Time-Periodic Magnetic Diffusion Problems Using Diagonally Implicit Runge-Kutta Methods[END_REF].

To overcome this drawback, a kind of semi-implicit method is used: the system of ODEs ( 3) is solved explicitly, whereas (12-13) is solved partially implicitly. ( X n+1 , d n+1 , ω n+1 ) are found with a time step of (3) from ( X n , d n , ω n ), with an explicit high order Runge-Kutta method (Dormand-Prince 5), where the components of force and torque computed with (10) are frozen at time n τ , i.e. computed with ( T n , h n s ). Then, as ( X n+1 , d n+1 , ω n+1 ) are considered as given, ˙ h * n+1 s can be deduced as ˙ h s ( X n+1 , d n+1 , ω n+1 ). (12-13) are then solved implicitly with the source term ˙ h * n+1 s , which allows to find ( T n+1 , Ω n+1 ).

Let us introduce the magnetic co-energy W , Joule losses p j , given supplied electrical power p s of induced currents

W = E3 µ 0 T + ∇Ω 2 d x 3 ; p j = D ( ∇× T ) 2 σ d x 3 p s = - D µ 0 ˙ h s • T d x 3 (14)
and kinetic energy E k , supplied mechanical power p m to the magnet

E k = m 2 ˙ X 2 + 2 5 R 2 w 2 ; p m = f • Ẋ + Γ • ω (15)
The weak forms of (6-9) and the Hamiltonian form of the dynamical system (3) allow us to write two half-power balance equations:

dW dt + p j = p s ; dE k dt = p m (16)
The coupling which results from the dependence of ˙ h s on mechanical variables (5) and of ( f , Γ) on electrical variables (10) corresponds to:

p s + p m = 0 =⇒ d dt W + E k + p j = 0 (17)
At the discrete level, the two-half power balances (16) are globally satisfied. But, due to the approximation on ˙ h s , this is less so for the total balance equation (17), the accuracy of this equality has to be monitored.

III. RESULTS

A spherical Nd-Fe-B magnet (R = 6.35 mm, M = 1 A m 2 ) is thrown on a copper plate (thickness 5 mm, σ = 50M S/m), by using a helical compression spring. It helps to ensure a reliable initial speed ˙ X(t = 0) = ẋ0 k x ( ẋ0 = 3.7 m/s). Even if it is a low effet negligible during it motion, due to the earth's magnetic field the equilibrium direction of magnetization of the motionless spherical magnet is θ 0 = 25 • in Nancy (France), which fixes one initial condition. The other one can be chosen as ψ 0 = 0, for which the shot direction is perpendicular to the north-south axis. The trajectory of the magnet has been filmed with a high-speed camera (500 frames/s, 6 px/mm). Figure 2 shows snapshots of the magnet from t = 0 to t = 100ms, every 20ms. The terminal position, t = 160ms, is represented by the cross.

The motion of the magnet has been computed with the help of (3-12-13) and compared with the experiment, to determine the validity of the model (Fig. 3). The model (plain line) is in a good agreement with the experimental data (dashed line): As can be seen on Fig. 3, in addition to the main motion along k x , there is also a motion of much smaller amplitude along k y (about 0.8cm). We have checked that this does not occur with a non-magnetic sphere. This lateral motion is due to the fact that θ 0 = 90 • . If this were not the case, the starting direction would be along k y , and Y would remain at zero. The experimental velocity Ẋ has been estimated with the change in position between snapshots (Fig. 4 dashed line). We have checked that the motion of the spherical magnet cannot be obtained with a single viscous-friction force along the x-direction.

The computed velocity Ẋ (Fig. 4 plain line) is of the same order of magnitude. At the end of the motion, there is a slightly higher difference between computation and experiment, two reasons are invoked. The first is the earth's magnetic field b e , whose torque Γ e = Ds M( d × b e ) d x 3 has to be taken into account. The intensity of the magnetic field due to the eddy currents is compared with the intensity of the earth's magnetic field (0.05mT ) (Fig. 5). Although this effect is not completely negligible, it is too small to explain the difference. The magnet's rotation due to the earth's magnetic field occurs in experiments when t = 150ms: both fields are of the same order of magnitude. The second reason is the frictional contact: the rolling friction (or adhesion phenomenon) plays a role when the velocity is low (An analogous phenomenon has been observed with a non-magnetic sphere ; we have not investiged further due to the complexity of these tribology aspects).

The chosen time step is τ = 1.2ms. In the worst case, when the velocity is a maximum (3.7m/s), it corresponds to 8 steps per rotation of the sphere. This is sufficient to ensure accuracy with a 5th-order Runge-Kutta method and does not lead to a too long time of calculation (2 hours with a 2.8GHz processor). To verify this statement we have computed the difference

E k (t) -E k (0) - t 0 p m (τ ) dτ (18)
and have obtained an error less than 0.1% at the end of the motion (E k (0) is the initial energy stored in the spring). If -p s is substituted to p m in (18), a bigger error is obtained, but it is less than 0.8% (Fig. 6): it validates the assumption of using Finally the starting direction of magnetization is a parameter which affects the trajectory. In the case considered here, the shot direction was chosen perpendicular to the north-south axis (ψ 0 = 0 • ). This introduces a symmetry default which generates a slight deviation along k y but d becomes quickly aligned to k y and then remains parallel to the plate.

If now a symmetric configuration is chosen, where the shot direction is aligned to the north-south axis (and then the initial direction of magnetization is in the ( k x , k z ) plane), a more complex kinematics is obtained (Fig. 7). When the initial velocity is high (3.7m/s) there is a complex rotating motion during the third of the time that leads d to be aligned with k y . But when this initial velocity is low (1m/s), d remains in the ( k x , k z ) plane: the sphere rolls without deviating.

The stopping distance (X when the magnetic field due to the induced current becomes equal to the terrestrial magnetic field) is approximately 24cm for the case (ψ 0 = 0 • , 3.7m/s); this distance is 20% smaller in the case (ψ 0 = -90 • , 3.7m/s). The same comparison made with an initial velocity 1m/s gives a ratio of 50%. The same ratios have been obtained by experiment.

IV. CONCLUSION

The major aim here was to solve a problem where an electromagnetic field is coupled to the motion of a solid particle. The difficulties lie essentially in the choice of the approximations on the coupling variables: force and torque on one side, position and speed on the other side. We decided to freeze these variables during a time-step for each problem, even if the accuracy of the numerical schemes differ. Comparison with experiment shows that this approach is adequate. In order to reduce the computation time, an embedded Runge-Kutta method with an energy criterion (difference between kinetic energy and accumulated electric power) will be tested.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Snapshots of the magnet at a set of discrete times (∆t = 20ms).

Fig. 4 .

 4 Fig. 4. Velocity Ẋ vs time (dashed line: experimental, plain line: computed).

Fig. 5 .

 5 Fig. 5. Magnetic induction due to the eddy currents in the center of the magnet µ 0 | ∇Ω| vs time (plain line) and earth magnetic induction (dashed line).

Fig. 6 .

 6 Fig. 6. Kinetic energy E k (dotted line), work of the electrical force (dashed line) and sum of both terms (plain line) versus time.

Fig. 7 .

 7 Fig. 7. Components of the magnetization direction vs time ( d• ky in plain line, d • kx dashed line, d • kz dotted line). Initial conditions: ψ 0 = -90 • , θ 0 = 25 • ; ẋ0 = 3.7m/s (bottom) and ẋ0 = 1.0m/s (top).