Jean-Guy Caputo
email: caputo@insa-rouen.fr

Denys Dutykh
email: denys.dutykh@univ-savoie.fr

Nonlinear waves in networks: model reduction for sine-Gordon

Keywords: sine-Gordon equation, kink, breather, networks, Josephson junction MSC: [2010] 35R02 (primary), 34B45 (secondary) PACS: [2010] 05.45.Yv (primary), 74.81.Fa (secondary) sine-Gordon equation, kink, breather, networks, Josephson junction

To study how nonlinear waves propagate across Y-and T-type junctions, we consider the 2D sine-Gordon equation as a model and examine the crossing of kinks and breathers. Comparing energies for different geometries reveals that, for small widths, the angle of the fork plays no role. Motivated by this, we introduce a 1D effective model whose solutions agree well with the 2D simulations for kink and breather solutions. These exhibit two different behaviors: a kink crosses if it has sufficient energy; conversely a breather crosses when v > 1ω, where v and ω are respectively its velocity and frequency. This methodology can be generalized to more complex nonlinear wave models.

Introduction

The propagation of nonlinear waves in networks is a very common problem. Examples are the nerve impulse traveling in arrays of neurons [14], the motion of the pulse wave in the arterial circulatory system [START_REF] Nichols | McDonald's Blood Flow in Arteries[END_REF] or the propagation of waves in the electrical power grid [START_REF] Backhaus | Getting a grip on the electrical grid[END_REF]. In general the problem is difficult to tackle because both the equation of motion and the geometry are complex. A first direction is to consider a simpler geometry like a Y-junction, see Fig. 1. Another simplification is to study what happens for a linear wave equation. In this context, a number of researchers have examined so-called quantum graphs where the Schrödinger equation is solved on a network. See [START_REF] Gnutzmann | Quantum graphs: Applications to quantum chaos and universal spectral statistics[END_REF] for a review. For these linear systems, the scattering formalism can be employed and this gives the reflection and transmission coefficients for a harmonic wave. This is detailed specifically for a Y-junction and for the Klein-Gordon linear wave equation in [3].

In many cases however, the nonlinearity cannot be neglected. For fluid systems, note the works by Bona and Cascaval [START_REF] Bona | Nonlinear dispersive waves on trees[END_REF] and Mugnolo and Rault [START_REF] Mugnolo | Construction of exact travelling waves for the Benjamin-Bona-Mahony equation on networks[END_REF] who used the Benjamin-Bona-Mahoney (BBM) shallow water equation to describe a fluid network. The authors used the fact that the BBM equation is unidirectional, hence, most of the energy is propagated downstream. For the Boussinesq equation in a junction, Nachbin and Da Silva Simoes [12] used a conformal map technique. However, all these studies do not provide a simple understanding of the behavior of the waves, in particular one cannot see easily how energy travels across the network.

To address these issues, a first step is to consider a simpler model. For instance, before tackling the propagation of shallow water waves in a river basin, for which there are two variables, the water elevation and the potential flow, it is useful to consider a simpler nonlinear hyperbolic equation. The sine-Gordon equation is precisely a simple nonlinear hyperbolic equation that admits localized solutions. It is also a Hamiltonian system in any dimension, integrable in 1D so that one can compare the numerical solutions with their exact counterpart as they propagate in a 1D channel. Finally, the sine-Gordon equation is an excellent model of an extended Josephson junction between two superconductors [2]. Here, to simplify the issue and to keep the Hamiltonian framework, we exclude external actions on the network, like a current or a magnetic field, that are commonly used in Josephson junction arrays.

Consider the 2D sine-Gordon equation defined in a Y-junction such as shown in Fig. 1. A first work on the problem is by Gulevich and Kusmartsev [8] who examined numerically how kinks propagate in such a system, in the context of Josephson junctions. They showed that the kink needs a sufficient velocity to cross the branch. Here we follow up on this and define a 2D symmetric junction parameterized by the angle θ between the branches and by the widths w 1,2 of the branches. This setup can describe a Y-junction up to a Tjunction (θ = 180 • ). We solve the 2D problem using the FreeFem++ finite element library [START_REF] Hecht | New development in Freefem++[END_REF]. It is important that the energy is conserved by the code; for this we found a suitable time discretisation. The study of the propagation of a kink in the junction yielded the

w 0 w 1 w 1 x y O ¡ 2 L Figure 1. Sketch of the computational domain Ω.
immediate result, that for small widths, there is no dependence of the velocity on the angle of the fork for the full 2D simulation. We therefore introduce a 1D effective partial differential equation to capture the essential features of the 2D propagation. This model incorporates the junction, using the ideas of graph Laplacian [5]; its solutions agree well with the 2D solutions. For the kink propagation in a junction we confirm the existence of a critical velocity given approximately by the simple energy conservation argument. Below this velocity, the kink gets reflected by the fork. Above it, it passes through the junction and splits into two kinks that propagate in the two different branches. For breathers there are two parameters, ω the frequency and v the velocity. For a given velocity v, junction crossing is only possible above the frequency ω ≈ 1 -v indicating a nonlinear resonance. After the breather passes through the junction it gives rise to new breathers in the branches. We characterize these using their energy density and estimate their velocity and frequency. We always observe an up-shift of the frequency and a slight down-shift of the velocity.

The article is organized as follows. In Section 2 we derive the 1D effective model from the 2D sine-Gordon equation defined in the fork. In Section 3 we recall the energies for the kink and the breather and show how they can be used to estimate a critical velocity. Section 4 introduces energy conserving discretisations for the finite element 2D problem and the 1D effective equation. Their solutions are compared in section V for both kink and breather initial conditions. Conclusions are presented in Section 5.

The 1D effective sine-Gordon model

We consider the 2D sine-Gordon equation

φ tt -∆φ + sin φ = 0, (2.1)
on a bounded domain Ω ⊂ R 2 with Neuman boundary conditions:

∇φ • n = 0,
where n is an exterior normal. The t subscript indicates the time derivative and ∆ is the usual Laplacian in spatial coordinates. This equation conserves the energy

E = ˆΩ 1 2 φ 2 t + 1 2 |∇φ| 2 + (1 -cos φ) dx dy. (2.2)
This can be checked easily by multiplying (2.1) by φ t , integrating over the domain and using the Stokes formula for the spatial operator.

Since the boundary conditions of the 2D problem are homogeneous Neuman it is natural to assume that the solution is uniform in the transverse direction. In other words we keep only the first transverse Fourier mode. Then equation (2.1) reduces in each branch to a 1D sine-Gordon equation,

φ i tt -φ i xx + sin φ i = 0, i = 1, 2, 3 , (2.3)
where the label i corresponds to the three branches as shown in Fig. 2. These equations are coupled at the apex by two conditions; one is the continuity of φ i

φ 1 (x = l) = φ 2 (x = 0) = φ 3 (x = 0) , (2.4)
and the other is the flux conservation or the Kirchoff law

-w 1 φ 1 x + w 2 φ 2 x + w 3 φ 3 x = 0, (2.5)
where φ i x is the normal velocity in branch i. Let us now briefly justify this flux relation. For that, consider the fork domain F obtained by taking the normals at the different branches as close as possible to the fork as shown in Fig. 3. Integrating the two-dimensional equation (2.1) on F yields

ˆF (φ tt + sin φ) dxdy - ˆ∂F ∇φ • n ds = 0, (2.6)
where n is the normal to the edge of the domain ∂F . The second term is equal to the left hand side of (2.5). The first one is of order w 2 where w is the typical width of the branches.

In the limit of small w, such that w i → 0 with w 2 /w 1 and w 2 /w 1 constant, the first term vanishes while the second one remains, yielding (2.5).

The numerical scheme used to solve this 1D effective model is described below (see Section 4); it is a finite difference approximation. The junction corresponds to the four nodes highlighted on Fig. 2; these are labeled as 1, 2, 3 for the three branches and are connected to the central node 4. The outer nodes are the last nodes updated by the PDE solver; let us name the value of the solution there φ 1 , φ 2 , φ 3 for each branch. The value at the central node φ 4 can be computed from the interface conditions (2.4) and (2.5). Using a forward finite difference approximation for φ i

x we get from (2.5)

-w 1 (φ 4 -φ 1 ) + w 2 (φ 2 -φ 4 ) + w 3 (φ 3 -φ 4 ) = 0,
where we have assumed the same space step on the three branches and used the notation φ i ≡ φ(x i ). We have also omitted the j index corresponding to the different branches. We then obtain

φ 4 = w 1 φ 1 + w 2 φ 2 + w 3 φ 3 w 1 + w 2 + w 3 .
(2.7)

Theoretical considerations

In 1D the sine-Gordon equation is integrable, see for example [14]. It has two families of localized exact solutions, the kink

φ(x, t) = 4 arctan [exp(γ(x -vt))] , (3.1)
φ(x, t) = 4 arctan √ 1 -ω 2 cos(ωγ(t -vx)) ω cosh( √ 1 -ω 2 γ(x -vt)) , (3.2)
where the Lorentz factor γ is given by

γ = 1 √ 1 -v 2 . (3.3)
Let us first consider the kink, its energy is

E k = 8γ. (3.4)
The energy of the breather depends also on the frequency, it is given by

E b = 16γ √ 1 -ω 2 . (3.5)
In two dimensions the equation is not integrable. In addition there is the complication of the boundaries. Therefore the only relations that can be used are conservation laws, and in particular the conservation of energy. When the kink is in branch 1, its energy is 8w 1 γ because it is homogeneous in the transverse direction. Similarly in branch 2, it has energy 8w 2 γ 2 . The conservation of energy reads

w 1 8 1 -v 2 1 = 2w 2 8 1 -v 2 2 , (3.6)
where we assumed w 2 = w 3 . This expression gives a critical velocity v 1 for which v 2 = 0:

v k = 1 - w 1 2w 2 2 . (3.7)
This formula was derived in [8] and compared successfully to the 2D numerical results for a fixed angle and widths w 3 = w 1 , 0 < w 2 /w 1 < 1. In the next section we confirm this estimate by numerical simulations and show its limitations.

A similar argument for the breather yields the following result for the parameters {v 1 , ω 1 } in the bottom branch and the parameters {v 2 , ω 2 } in the top branches

v 2 1 -1 ω 2 1 -1 = w 1 2w 2 2 v 2 2 -1 ω 2 2 -1 . (3.8)
This gives a critical velocity v 1 for which v 2 = 0

v k = 1 - ω 2 1 -1 ω 2 2 -1 w 1 2w 2 2 .
(3.9)

The practical application of the previous formula is difficult because ω 2 remains unknown. Note however that for small amplitudes, i.e. in the linear limit ω 1 = ω 2 so that we recover (3.7) for the critical velocity.

Numerical methods

We first consider the propagation of a kink in Y-and T-junctions. As expected the kink gets reflected if it does not have enough energy (velocity). Also the motion depends very weakly on the angle. To illustrate this fact, we show in Fig. 4 a kink propagating in a T-junction and crossing it. We take the same kink and run it into a Y-junction. This is shown in Fig. 5. One can see that the time intervals for propagation are about the same. This is confirmed by examining the evolution of the energy in the branches 1 (bottom) and 2 (left), see Fig. 6. Note also that a very small amount of energy, typically 5% of the total energy, is left in branch 1 once the kink has crossed over into branches 2 and 3. The solution of the 1D effective model is plotted with points in Fig. 6; it agrees very well with the 2D solution and this confirms the informal asymptotic reduction from 2D to 1D of the previous section.

We now compare systematically the 2D solution with the one of the 1D effective equation. To validate this approximation, we conduct a parametric study choosing w 1 , w 2 and w 3 such that w 1 = 1, w 2 = w 3 = w 1 + α where α = -0.3, -0.1, 0.1, 0.3. The results for the critical velocity as well as the estimate (3.7) are reported in the Table 1. The 2D and 1D models are very close even for α > 0. On the other hand the energy estimate is a lower estimate for α > 0, The 2D and 1D effective results reveal that the kink crosses the junction but that there are oscillations. The front seems to oscillate and then reshape as it enters more into branches 2 and 3. We do not see this effect when α ≤ 0. Despite this, the values are all within a 10% interval of error.

To check the 2D-1D reduction even further, we conducted simulations with a large width w 1 = w 2 = w 3 = 10 for large and small velocities. We found that for a small velocity v 1 = 0.75 the kink gets reflected and significant transverse oscillations occur. These transverse oscillations are the ones that propagate along the equal phase contour lines; they were studied by Gulevitch et al. [START_REF] Gulevich | Shape Waves in 2D Josephson Junctions: Exact Solutions and Time Dilation[END_REF] who showed that their dispersion is ω = k. In this situation, the angle becomes important and of course the 2D and the 1D models disagree. On the other hand, for a large velocity v 1 = 0.96, the kink crosses and the transverse oscillations remain small. Then the 2D and the 1D models are very close.

For the breather, things are more complicated because of the additional parameter, the frequency. The energy criterion is not sufficient, the breather needs to have the adapted frequency in order to cross. To illustrate this, we consider a breather of initial velocity v 1 = 0.4 and different frequencies. Fig. 7 shows the energy E 1 in branch 1 as a function of time for frequencies ω 1 = 0.5, 0.7, 0.725 and 0.75. The breather does not cross until the frequency reaches 0.75; for this frequency the initial energy is smaller than for ω 1 = 0.5. The initial energy will not determine crossing, the initial frequency will. Another interesting effect is that for ω 1 = 0.725 the energy E 1 almost goes to 0 in the time interval [200; 300] -as if the breather crossed the junction-and then it returns almost to its initial value, indicating reflection. We checked the position of the breather and found that it stays with an interval of size 10 close to the junction. This long oscillation close to the junction could indicate a nonlinear bound state associated to the fork. Fig. 8 shows the crossing vs. non crossing in the parameter space (v, ω). The crossing (resp. reflection) of the breather is indicated by the + (resp. ×) sign. The calculations were done both with the 1D effective model and the full 2D equations and the results always agreed. For large enough velocities, the breather crosses independently of its frequency. On the other hand, for frequencies close to one, the breather crosses even for small velocities. This situation is close to the linear case for which we expect always some energy transfer to the other branch [3].

There is always a small reflection from the fork. For example, we show the time evolution of the energy of a breather in Fig. 9. Notice how the energy in branch 1 does not drop to zero as for the kink. There is a remainder.The small scale oscillations present in Fig. 9 are the ones seen in Fig. 7; at this time we do not know their origin. To characterize the breathers in the other branches is difficult because the wave oscillates. We found that plotting the energy density

dE = 1 2 φ 2 t + 1 2 φ 2 x + 1 -cos φ, (4.1)
gives a good indication of the position of the breather. Let us analyze in more details the specific configuration where a breather of speed v = 0.8 and frequency ω = 0.3 crosses the junction. Fig. 10 shows the energy density for three different times in the branch 1 (top panel) and in the branch 2 (bottom panel) after the breather has passed the junction. Then the energies in branch 1 and branch 2 are respectively E 1 = 2.16 and E 2 = 13.23.

The velocities estimated by a least square fit on the center of mass of the breather density are respectively v 1 = -0.75 and v 2 = 0.6. They are lower than the initial velocity to accommodate for the crossing of the breather. The frequencies of the breathers in branches 1 and 2 can be estimated, they are respectively ω 1 = 0.996 and ω 2 = 0.75. All these parameters are very different from the initial breather parameter making the scattering of a breather much more complex than the one of a kink. Using the parameters above we can plot the fitted breathers and compare them with the numerical solution. Fig. 11 shows in the top panel, branch 1 before the breather crosses. There the analytical solution matches perfectly the numerical one. The middle and bottom panels show respectively branch 1 after the crossing and branch 2. Here the agreement is not as good but remains quite acceptable. The reflected breather in branch 1 (middle panel) has a small amplitude of about 0.2. Its frequency as seen from Table 2 is 0.99 so that it is very close to the dispersion curve ω = √ 1 + k 2 . This can explain the dispersion observed.

To conclude this study we examine systematically the influence of the breather frequency on its crossing. We took v 1 = 0.8 and chose ω 1 = 0.3, 0.5, 0.7 and 0.9. The results are reported in Table 2.

Conclusion

We analyzed numerically and theoretically how a 2D sine-Gordon kink or breather crosses a Y-or T-junction. Comparing the energies in the different branches for both cases revealed that the angle of the junction plays almost no role in the dynamics for thin trees. This suggested to introduce a 1D effective model where at the junction, we satisfy continuity of the solution and flux conservation. The solutions of this effective model accurately reproduce the 2D solutions.

The parameters for the kink to cross obey the simple relation obtained from the conservation of energy. There is a critical velocity below which no crossing is possible.

Breather crossing is more complex because there are two parameters: the velocity v and frequency ω. For equal widths of the branches, we observe crossing when v > 1 -ω. Table 2. Velocities and frequencies for the crossing of a breather of initial velocity v 1 = 0.8 and different frequencies ω 1 = 0.3, 0.5 (top rows) and ω 1 = 0.7, 0.9 (bottom rows). The columns indicate the branches, 1, "1-return" and 2. The label "1-return" corresponds to branch 1 after the collision. (2) Then the breather gives rise to two other breathers in the two upper branches that we characterize using the energy density and the value of the energy. These new breathers are slower than the initial condition and are also up-shifted in frequency. We always observe a small reflexion at the crossing into the first branch. This study can be extended by considering more branches. Another interesting extension would be to add a source at the junction, enabling to control the crossing. It would be useful to understand how to transpose this to another application, like the reflexion of shallow water waves.

Figure 2 .

 2 Figure 2. Sketch of the tree geometry for the 1D effective model.

Figure 3 .

 3 Figure 3. The fork region F .

Figure 4 .Table 1 .

 41 Figure 4. Motion in a T-junction. Snapshots of a kink starting in branch 1 with a velocity v 1 = 0.75. The times are t = 1350 (a) and t = 3000 (b) .

Figure 5 .

 5 Figure 5. Motion of a kink in a 90 • Y-junction. Snapshots of a kink starting in branch 1 with a velocity v 1 = 0.75. The times are t = 900 (a), t = 1500 (b) and t = 3500 (c). The widths are w 1 = 1, w 2 = w 3 = 0.7.

Figure 6 .

 6 Figure 6. Time evolution of the energy for the kink motion in branches 1 and 2 for the T-junction in full line (red online) and for the Y-junction in dashed line. The energy for the 1D effective model is plotted with points. The parameters are the same as in Figs. 4 and 5.

Figure 7 .

 7 Figure 7. Energy in branch 1 as a function of time for the 1D model and a breather of initial velocity v 1 = 0.4 and frequencies ω 1 = 0.5, 0.7, 0.725 and 0.75. Notice the crossing for ω 1 = 0.75. The widths are w 1 = w 2 = w 3 = 1.

Figure 8 .

 8 Figure 8. Parameter space (v, ω) for the crossing of breathers for the 2D and the 1D models. The crosses (red online) indicate crossing and the squares reflection. The widths are w 1 = w 2 = w 3 = 1.

Figure 9 .

 9 Figure 9. Time evolution of the energy in branches 1 and 2 for a breather for the 2d partial differential equation in full line (red online) and the 1d effective model in dashed line. The parameters are w 1 = w 2 = w 3 = 1, v 1 = 0.8, ω 1 = 0.3, x 0 = 10.

Figure 10 .

 10 Figure 10. Snapshots of the energy density of a breather at different times in branches 1 (top) and 2 (bottom). The times are indicated on the plots. The parameters are w 1 = w 2 = w 3 = 1, v 1 = 0.8, ω 1 = 0.3. The initial position of the breather is x 0 = 10.

Figure 11 .

 11 Figure 11. Snapshots of the breather analytical solution (dashed line)together with the numerical solution (continuous line), in branch 1 before the collision (top panels a,b), in branch 1 after the collision (middle panels c,d) and in branch 2 (bottom panels e,f ). The corresponding times are t = 20.2, 40.4, 80.8, 90.9 for panels (a, b, c, d) and t = 80.8, 90.9 for panels (e,f ).

(bottom panels e,f ). The corresponding times are t = 20.2, 40.4, 80.8, 90.9 for panels (a, b, c, d) and t = 80.8, 90.9 for panels (e,f ).December 2012. The authors thank D. Mitsotakis and G. Sadaka for useful discussions on the finite element numerical method. The authors acknowledge the Centre de Ressources Informatiques de Haute Normandie where most of the calculations were done.

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgments D. Dutykh acknowledges the support from ERC under the research project ERC-2011-AdG 290562-MULTIWAVE and thanks INSA de Rouen for its hospitality during his visit in


Damien Borgetto
email: oatao@listes-diff.inp-toulouse.fr

An Energy Efficient Approach to Virtual Machines Management in Cloud Computing

come   L'archive ouverte pluridisciplinaire

OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible. This is an author-deposited version published in : http://oatao.univ-toulouse.fr/ Eprints ID : 15198

The contribution was presented at CloudNet 2014: http://cloudnet2014.ieee-cloudnet.org/


Christophe Mollaret
email: mollaret@irit.fr

Fréderic Lerasle
email: lerasle@laas.fr

Isabelle Ferrané
email: ferrane@irit.fr

Julien Pinquier
email: pinquier@irit.fr

A PARTICLE SWARM OPTIMIZATION INSPIRED TRACKER APPLIED TO VISUAL TRACKING

Keywords: visual tracking, particle swarm, particle filter, video analysis, RGB-D sensors

Visual tracking is dynamic optimization where time and object state simultaneously influence the problem. In this paper, we intend to show that we built a tracker from an evolutionary optimization approach, the PSO (Particle Swarm optimization) algorithm. We demonstrated that an extension of the original algorithm where system dynamics is explicitly taken into consideration, it can perform an efficient tracking. This tracker is also shown to outperform SIR (Sampling Importance Resampling) algorithm with random walk and constant velocity model, as well as a previously PSO inspired tracker, SPSO (Sequential Particle Swarm Optimization). Experiments were performed both on simulated data and real visual RGB-D information. Our PSO inspired tracker can be a very effective and robust alternative for visual tracking.

INTRODUCTION

Visual tracking is one of the most studied topics in the research field of video analysis [START_REF] Smeulders | Visual tracking: An experimental survey[END_REF]. It has many uses such as video surveillance, vision-based control and human-computer interface. These techniques are more and more used in human-Machine Interaction (HMI) [START_REF] Chen | Real-time visual tracking with adaptative particle filter for humanmachine interaction[END_REF], which needs to fuse a lot of visual information which is given by many sensors and/or algorithms, which will also be considered as sensors. For our application, we need to characterize the user's will to communicate with a machine also known as "awareness" in Human-Machine Interaction community [START_REF] Drury | Awareness in human-robot interactions[END_REF].

In visual tracking, two kinds of tracking algorithms are mainly used. They are Kalman filters and Particle filters. Kalman filters are widespread in industry and research for multiple use such as 3D-modeling, path following and object tracking. Chen in [START_REF] Chen | Kalman filter for robot vision: A survey[END_REF] has made a survey on different techniques and algorithm implementations. However the original Kalman algorithm is limited to linear problems and be-This work is partly supported by a grant from the ANR (Agence Nationale de la Recherche) with reference ANR-12-CORD-003.

comes more complex if we need to apply variants of this algorithm (unscented, extended) on multimodal and non-linear cases. Meanwhile, particle filtering is a very simple algorithm adapted to multimodal and non-linear target tracking. Nevertheless, it requires an increase in the number of particles to perform efficient estimation, raising its complexity exponentially with the dimension of the estimated state vector. This effect, known as dimensionality curse, is very important in particle filtering techniques. This becomes a problem in visual tracking task since we need to approach real time tracking by decreasing the computing time, i.e. decreasing the number of particles.

For many years, Particle Swarm Optimization (PSO) has been gaining attention [START_REF] Poli | Analysis of the publications on the applications of particle swarm optimisation[END_REF], [START_REF] Kulkarni | Particle swarm optimization in wireless-sensor networks: A brief survey[END_REF], since this algorithm is able to solve non-linear, multimodal and high-dimensional optimization problems. Contrary to other particle-based algorithms, particles interact with each other with a "social" and "cognitive" component in the update equation. This leads to finer particle behaviour during the process of optimization. Furthermore, the computation cost remains low since only one computation of fitness function per iteration and per particle is needed, which is, usually, the bottle-neck part in visual tracking problems.

In this paper, a new approach on PSO based tracking is proposed, merging particle filter versatility with a lower computational cost, adapting to real time constraints. The idea is also to keep the "PSO behaviour" of particles which seems more interesting than the usual "particle filter behaviour". Because of its social and cognitive components, particles of PSO algorithm interact with each other, contrary to particle filter. We think that this behaviour could lead to a more efficient estimation of the target as there is no particle degeneration in PSO. Therefore, a resampling step becomes needless and we can boost the algorithm execution. Finally, like the particle filter with constant velocity model, the PSO has an implicit estimation of state vector speed, which seems to model the system dynamics better than a standard particle filter. But contrary to the former particle filter, we do not need to double the state vector dimension to perform the estimation, since only the speed in PSO has a random component.

In a first section we will present some related works, then, we will focus on the formalism of the Particle Swarm Optimization algorithm and our new PSO-based tracker. Next, we will implement and compare four algorithms on simulated data to evaluate the robustness and efficiency of our new tracker. Finally we will perform tracking on real recorded RGB-D data with ground truth in our applicative context, human awareness characterisation in HMI. We used two modalities outputs which worked for us as "sensors". The first one was provided by OpenNI/NITE library [START_REF] Primesense | Openni/nite, open-source sdk for 3d sensors[END_REF] for the shoulder orientation and the detection of head position. The second sensor was developed by Fanelli in [START_REF] Fanelli | Real time head pose estimation with random regression forests[END_REF] for research purposes and enables the full head pose detection on depth informations.

RELATED WORKS

Particle Swarm Optimization algorithms has been used in a wide range of application for many years. Poli in [START_REF] Poli | Analysis of the publications on the applications of particle swarm optimisation[END_REF] proposed an overview of them, and it appears to have been implemented for parameter tuning in particle filter, fuzzy systems, video analysis, and image analysis. This algorithm presents good performances with high dimension problems, and seems to suffer less from dimensionality curse. Moreover, it is computationally lower than Sequential Monte-Carlo (alias SMC) methods such as particle filtering applied in optimization problems. However, this algorithm belongs to the meta-heuristics branch such as genetic algorithms, and its convergence cannot be proven.

For past years some investigations have aimed at adapting PSO in a tracking framework. Zhang first implemented it in [START_REF] Zhang | Sequential particle swarm optimization for visual tracking[END_REF], where he added PSO iterations in a particle filtering framework. However, the computation cost is proportional to iteration numbers and the algorithm finally become much more costly than particle filtering techniques, despite the fact that it performs a more precise tracking. Recently, Li in [START_REF] Li | A video tracking method based on niche particle swarm algorithm-particle filter[END_REF] replaced the resampling step in particle filter by one PSO iteration. This is as computationally costly as particle filters. However, they lose the formalism of the resampling step which has to conserve the estimated MMSE while repositioning particles.

Particle filter is also one of the only particle based tracker, and is widely used in the vision literature for visual tracking problems. There are many variants of this algorithm and one of them is used by Sedai in [START_REF] Sedai | Robust perception of an interaction partner using depth information[END_REF] for human pose tracking which is getting closer to our context. The Annealed Gaussian Process Guided Particle Filter (AGP-PF) that they used is better than basic Particle Filter, but much heavier computationally speaking and complicated to implement.

Thus, a PSO-based tracker is created, more in original PSO algorithm respect, mainly in terms of particles behaviour. It is also wanted lighter in computational resources than a particle filter while being more efficient. Finally, the algorithm aims to remain simple, while some improved Par-ticle Filter implementations become very heavy in code and processing cost.

FORMALISM

Four algorithms are compared for this work. Two SIR trackers with two different dynamic models, and two PSO based trakers namely [START_REF] Zhang | Sequential particle swarm optimization for visual tracking[END_REF] and ours.

SIR-RW: The first one, the more usual, is a SIR-RW where "RW" stands for Random Walk dynamic model.

SIR-CV:

The second one is a SIR-CV, where "CV" abbreviates Constant Velocity dynamic model. Contrary to the previous one, the state vector dimension is doubled in order to add its first derivative as a component of the state vector. Thus, the state vector velocity is also estimated. When s (i) t is the state vector in SIR-RW algorithm, we have to model the SIR-CV state vector by x

(i) t = [s (i) t ṡ(i) t ],
(i) t = s (i) t-1 + ω * ṡ(i) t-1 + w t where w t is gaussian noise as well, s (i)
t is the state vector position part and ω is a multiplicative constant modelling inertia.

SPSO:

The third algorithm is the SPSO (Sequential Particle Swarm Optimization) as formalised by Zhang in [START_REF] Zhang | Sequential particle swarm optimization for visual tracking[END_REF], and its implementation will not be recalled in this paper.

TPSO: Finally, the fourth one is our modified PSO for tracking (alias TPSO). This algorithm is described in details in this section after a short reminder on PSO algorithm for a better understanding of our contribution.

In our approach, optimization and tracking are partially the same problem: in optimization, the algorithm has to converge to a fixed global maximum in a certain number of iterations. However, in tracking, also called dynamic optimization, the algorithm has to converge to a moving global maximum (the target), and iterations are represented by frames (one frame = one iteration). Regarding this, SPSO is not a proper tracker since it iterates more than once for each frame. Considering this, we took the PSO and modified the parts were the target is mentionned in the equations to take into account the target motion. In algorithm 1, the PSO algorithm principe is reminded. This algorithm being adapted for regular optimization and not dynamic optimization, algorithm 2 shows our adaptations to tracking.

In algorithms 1 and 2, U (a, b) is a uniform distribution between a and b, N is the number of particles, b lo and b up are lower and upper boundaries of search space, f (.) is a fitness function, or likelihood for SIR algorithms. r p and r g are factors that respectively ponderate social and cognitive terms randomly in the next line. ω models particles inertia. φ p and φ g are constant weights of cognitive and social factors respectively. The criterion can be a maximum number of iterations or a threshold on the MAP estimate fitness f (g).

Algorithm 1: Particle Swarm Optimization (PSO)

Result: MAP estimate : x = g for i=0 to N (initialization) do

1 x i ∼ U (b lo , b up ) 2 p i ← x i , g ← argmax(f (p i )) 3 v i ∼ U (-|b up -b lo |, |b up -b lo |) 4
while criterion is not met (iterations for optimization) 5 do for i=0 to N do 6 r p , r g ∼ U (0, 1)

7 v i ← ωv i + ψ p r p (p i -x i ) + ψ g r g (g -x i ) 8 x i ← x i + v i 9 if f (x i ) > f (p i ) then 10 p i ← x i 11 if(f (p i ) > f (g)) g ← p i 12
In our algorithm 2, lines 8 and 10 of PSO are modified to add a "prediction" component represented by a dynamic model, and there is now only one iteration per frame as explained before. The d(.) function represents the filter dynamic model. For our experiments, a random walk was considered(d(p i ) = p i + w where w is a gaussian noise). We could have set the dynamics model differently like the constant velocity model used in the second particle filter presented in this paper but chose not to. We added this function assuming that the optimal positions relative p i , and global g change, as the target position changes each frame. Then, f (p i ) has to be re-computed as it will have moved since the last measurement. These two minor changes allowed us to create a robust dynamic optimizer (or tracker), using PSO formalism, keeping the algorithm extremely simple.

EVALUATIONS ON SYNTHETIC DATA

Our new tracker was then processed on synthetic signals while one free parameter was changed at a time, to evaluate and compare its accuracy versus the other algorithms. The synthetic signal is a multidimensionnal signal constructed from concatenation of sinus functions, fast value changes and straight lines, in order to make a challenging target to follow.

Our observation model is a multivariate gaussian with a diagonal covariance matrix. The measurements are samples of synthetic signals added to a gaussian random noise with the same covariance matrix as the observation model, to make this optimal. The similar parameters of the four trackers are set equals: number of particles N , process noise σ in dynamic model for the four algorithms; ω for SPSO, TPSO and SIR with constant velocity model; ψ p and ψ g for SPSO and TPSO.

The following effects of parameters are evaluated on tracking process: particle number, noise power added to the observations and finally, dimension of the state vector. More-

Algorithm 2: PSO based Tracker

Result: MAP or MMSE estimate for each frame for i=0 to N (initialisation)

do 1 x (0) i ∼ U (b lo , b up ) 2 p (0) i ← x (0) i , g (0) ← argmax(f (p (0) i )) 3 v (0) i ∼ U (-|b up -b lo |, |b up -b lo |) 4
for t=1 to frame M (iterations on video sequence) do 5 for i=0 to N do 6 r p , r g ∼ U (0, 1)

7 v (t) i ← ωv (t-1) i + ψ p r p (d(p (t-1) i ) -x i ) + 8 ψ g r g (d(g (t-1) ) -x i ) x (t) i ← x (t-1) i + v (t) i 9 if f (x (t) i ) > f (p (t-1) i ) then 10 p (t) i ← x (t-1) i 11 if(f (p (t) i ) > f (g (t-1) )) g (t) ← p (t) i 12 MAP estimator : x(t) = g (t) or 13 MMSE estimator : x(t) = N i=0 f (pi) (t) N i=0 f (pi) (t) x (t) i 14
over, we computed the N ef f estimate (Number of efficient particles), in SIR and TPSO. In SIR, it is given by the equation

N ef f = 1/( N i=1 w 2 i ). In TPSO, an equivalent to N ef f is the equation N ef f = ( N i=1 f (x (i) p )) 2 / N i=1 f (x (i) p ) 2 .
We did not implement this estimator in SPSO as there is more than one X (i) p estimation per iteration and, therefore, no true equivalent to usual N ef f . In Fig. 1. are displayed the RMSE error means of 50 runs. TPSOmap and TPSOmmse are respectively MAP and MMSE estimator of the same TPSO run due to the stochastic nature of the process.

Noise: First we can see (Fig. 1 (c)) that our tracker reacts better to noise, with a tracking error divided by 2 for low SNR, where SNR is Signal to Noise Ratio in dB: SN R = 10log(P owerSignal/P owerN oise).

Particles: Then, we can notice (Fig. 1 (a)) that the number of particles can be divided by 5 while keeping the same precision as SIR algorithms. Thus, even tough we have to compute likelihood two times by iteration and particle, we can have a reduced complexity by reducing tremendously the number of particles.

Dimension: Moreover, we can see (Fig. 1 (b)) that our tracker also better reacts when state vector dimension increases, and estimates the target state vector with a reduced variance in high dimensions.

Neff: Finally, we can explain these results with N ef f estimator for TPSO and SIR algorithms (Fig. 1 (d)). For 500 particles, SIR-RW N ef f degenerates to 50, while SIR-CV acts better with 450 particles, and ours achieve the best with 490 particles. This result is very interesting considering that our algorithm has a simpler dynamic model than SIR-CV and does not need any resampling step to avoid degeneration.

APPLICATION TO VISUAL TRACKING

In our context (cf. Fig. 2), we want to track head pose (yaw, pitch, roll), and shoulder orientation in space as they are key information of machine awareness. Full video sequences can be found here 1 . These modalities are given by OpenNI [START_REF] Primesense | Openni/nite, open-source sdk for 3d sensors[END_REF] for shoulder orientation and head position, and Fanelli [START_REF] Fanelli | Real time head pose estimation with random regression forests[END_REF] for head pose. However we had a problem with angles as they vary from -180 to +180 and can jump between these extremes values. We needed to take care of it during sampling phase. This led us to a The dynamic model is a random walk for all filters, excepted for the SIR-CV, where there is a speed component. The process gaussian noise is modelled by a diagonal covariance matrix as well. Parameters defined in previous sections are set to: ω = 0.9, ψ p = 0.8, ψ g = 1, and N = 100 particles. Finally, data acquisition was performed by a calibrated RGB-D sensor. The ground truth was created with a marker-based Motion Capture (Mocap) disposal which is also calibrated and synchronized temporally with the RGB-D sensor.

A set of 100 runs for each tracker were carried out on the same path, and results are plotted on Fig. 3, representing the tracking error means. The plot (a) displays RMSE mean position error in mm, and (b), the RMSE mean angular error in degrees, since a global error on position and angles would have no sense. We can see that our algorithm achieved better performances on both position and angles, while particle filtering techniques had problems to perform an accurate tracking jointly on position and angles. These results came to comfort our simulation results and showed that our tracker was a suitable solution for robust tracking, keeping a low computational cost with only 100 particles and a relative complexity of 2/5 of the SIR, based on how often the observation function is called in a loop.

CONCLUSION

In this paper, we proposed a new tracking framework based on Particle Swarm Optimization formalism and, demonstrated on both, synthetic and real data, that our algorithm outperforms SIR algorithm while remaining as simple in implementation, and faster than SPSO. Our algorithm also seems to be more robust to state vector dimensionality, noise in observations, and an efficient tracking can be performed with fewer particles, as each one of them carries more information. All of these aspects present our PSO-inspired tracker as a promising alternative for visual tracking and others applications.

In future work, we intend to evaluate the influence of each TPSO parameter on tracking performance and will try to apply it to multi-target tracking.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Simulation results

Fig. 3 .

 3 Fig. 3. Tracking results

  where t is time and i the particle index. The update equations become ṡ(i)

t = ṡ(i) t-1 + v t where v t is a gaussian noise and ṡ(i) t the state vector velocity part. s

  7-state vector s = [x y z θ head φ head ψ head θ shoulders ]. The observation model is a multivariate gaussian model with a diagonal covariance matrix diag(Σ) = [diag(Σ position ) 1 https://projects.laas.fr/riddle/riddle-project/papers/

		2000											100
																						SIR-RW
	RMSE error mean in millimeters	600 800 1000 1200 1400 1800 1600					SIR-RW SIR-CV TPSOMap TPSOmmse SPSO				RMSE error mean in degrees	80 0 20 40 60									SIR-CV SPSO TPSOMap TPSOmmse
		400
													-20
		200
		20 0	22	24	26	28	30	32	34	36	38	40	20 -40	22	24	26	28	30	32	34	36	38	40
						Frame number									Frame number
		(a) Position tracking error		(b) Angular tracking error