Jean-Philippe Boucher

Christophe Clanet

David Quéré

Frédéric Chevy

The popsicle-stick cobra wave

The cobra wave is a popular physical phenomenon arising from the explosion of a metastable grillage made of popsicle sticks. The sticks are expelled from the mesh by releasing the elastic energy stored during the weaving of the structure. Here we analyse both experimentally and theoretically the propagation of the wave-front depending on the properties of the sticks and the pattern of the mesh. We show that its velocity and its shape are directly related to the recoil imparted to the structure by the expelled sticks. Finally we show that the cobra wave can only exist for a narrow range of parameters constrained by gravity and rupture of the sticks.

Physics of metastable states is a classical topic of statistical physics [START_REF] Landau | Statistical Physics[END_REF][START_REF] Langer | Statistical theory of decay of metastable states[END_REF]. A well known route to relax towards equilibrium is via a nonlinear front which propagates with a constant speed, such as in viral spread [START_REF] Fisher | The wave of advance of advantageous genes[END_REF][START_REF] Kolmogorov | Study of the diffusion equation with growth of the quantity of matter and its applications to a biological problem[END_REF], bio-chemical reactions [5] or combustion [START_REF] Ya | The mathematical theory of combustion and explosions[END_REF][START_REF] Clavin | Combustion waves and fronts in flows[END_REF]. In mechanics, the domino race provides an example of such a process for a non-connected network [START_REF] Stronge | The domino effect: a wave of destabilizing collisions in a periodic array[END_REF][START_REF] Stronge | The domino effect: successive destabilization by cooperative neighbours[END_REF]. For entangled structures, the question of the optimization of the strength of grillages has been addressed [START_REF] Rozvany | Grillages of maximum strength and maximum stiffness[END_REF][START_REF] Tarnai | Duality between plane trusses and grillages[END_REF] especially because of its role in construction [START_REF] Stylianidis | Robustness assessment of frame structures using simplified beam and grillage models[END_REF] but their stability remains an open question. The same type of question also arises in biological systems, such as in the microtubule catastrophe [START_REF] Cheeseman | Molecular architecture of the kinetochore-microtubule interface[END_REF][START_REF] Desai | Microtubule polymerization dynamics[END_REF]. Microtubules are assemblies of GDP tubulin arranged in a tubular shape ending with a cap of GTP tubulin. The loss of this cap triggers a rapid depolymerization driven by the release of the stored mechanical strain [START_REF] Mitchison | Dynamic instability of microtubule growth[END_REF][START_REF] Mitchison | Microtubule assembly nucleated by isolated centrosomes[END_REF][START_REF] Bowne-Anderson | Regulation of microtubule, growth and catastrophe: Unifying theory and experiment[END_REF]. Here we study a macroscopic version of such a system, namely the so-called "popsiclestick cobra wave" [18], obtained by releasing a mesh of sticks woven according to Fig. 1.b.

To generate a cobra wave, the whole structure is loaded by the geometrically-constrained bending of the individual sticks and is held together by the red and blue sticks at the end of the mesh (see Fig. by one the freed sticks. Due to the asymmetry of the weaving, two very different dynamics occur depending on which stick was initially removed. When the red stick is taken away first, the sticks are expelled upwards and by reaction they pin down the rest of the mesh to the ground (see [19]). The outcome is dramatically different when the blue stick is removed. In this case, the sticks are expelled downwards and they raise the whole structure as presented in Fig. 2.a (see also the movie provided in [19]). After a few hundreds of milliseconds, the shape of the wave reaches a steady state (Fig. 2.b), and propagates at a few meters per second. Both the shape and the velocity remain the same until the wave-front reaches the end of the grillage (Fig. 2.c). In this letter, we combine experimental and theoretical approaches to characterize the velocity and shape of the cobra wave in the steady state.

Since the lifting force raising the lattice originates from the recoil imparted by the expelled sticks, the global dy-namics of the wave is set by the ejection rate γ and the momentum transferred during the expulsion M v where M is the mass of a stick and v the velocity of a stick right after expulsion. The time γ -1 taken by a stick to exit the mesh is given by L/v, where L is the length of individual sticks (Fig. 1.a). Taking E the Young's modulus, w the width and e the thickness of individual sticks (Fig. 1.a), v can be estimated from the balance between the kinetic energy M v 2 of a stick after ejection and the bending energy Ewe 5 /L 3 stored in each stick blocked by the lattice.

From the previous scaling analysis, we readily deduce the velocity v 0 of the wave-front. Indeed, since the sticks are expelled one by one, we have v 0 = aγ/2 where a is the spatial period of the pattern (Fig. 1.b). Noting that, up to a geometric factor depending on the angle θ, we have a ∝ L, both v and v 0 scale as

v v 0 = b(θ) E ρ e L 2 . ( 1
)
where ρ is the mass-density of a stick and b(θ) is a scaling factor that depends on the geometry of the mesh. With c = E/ρ the speed of sound in the material, we find that v 0 ∝ c(e/L) 2 . In Fig. 3, we confirm experimentally this scaling for six kinds of wooden sticks (the values of the mechanical and geometric parameters of the different stick models are given in [19]) and we observe that indeed v ∝ v 0 . As expected, the speed does not depend on the width of the sticks and increases quadratically with the ratio e/L. In Fig. 3.i, one can see that the speed of the cobra decreases with the angle of the lattice θ. This trend can be easily understood qualitatively by noting that the velocity of the wave is proportional to the spatial periodicity a = L cos θ/3. We now focus on the shape of the wave. The height of the cobra can be understood quantitatively within a generalized version of Euler's elastica theory, where the mesh profile results from a competition between elasticity, gravity and recoil imparted by the expelled sticks [START_REF] Audoly | Elasticity and geometry: from hair curls to the non-linear response of shells[END_REF]. We treat the mesh as a linear continuous medium characterized by a flexion modulus K [19] and we describe the expulsion of the sticks by a force F 0 and a torque C 0 exerted at the free end of the grillage. Describing the shape of the mesh by a profile r(s, t), where s is the curvilinear abscissa (Fig. 4), the local force F (s, t) and torque C(s, t) are given by

F = -K∂ 3 s r C = K∂ s r × ∂ 2 s r. (2)
In steady state, the shape of the cobra is constant and moves at the velocity v 0 . We therefore have r(s, t) = r(s = s-v 0 t) and writing Newton's Law for an infinitesimally small element of the mesh leads to the following dynamical equation where µ = 2M/a is the linear mass density of the cobra, T the longitudinal tension, τ the tangent unit vector and R the ground reaction. We assume that the contact with the ground occurs for s ≤ 0, so that z(s ≤ 0) = 0 and R(s ≥ 0) = 0, and s max is the total mesh-length rising above the ground. Projecting Eq. ( 3) on the tangent and normal directions, these equations can be recast into a closed equation for the curvature Γ = |∂ 2 s r| 1 2

µv 2 0 ∂ 2 s r = µg -K∂ 4 s r + ∂ s (T τ ) + R. (3)
d 3 Γ 2 dα 3 + dΓ 2 dα = µg K 2 sin α Γ + cos α Γ 2 dΓ dα . ( 4
)
where α is the local angle between the mesh and the horizontal axis. This equation is of third order in α and thus requires three boundary conditions to be solved. We obtain these conditions by writing the stress at the free end α max = α(s max ) of the mesh, namely:

C 0 = KΓ (5) F = T + KΓ 2 (6)
F ⊥ = - K 2 dΓ 2 dα ( 7
)
where F = F 0 • τ and F ⊥ = F 0 • n and the right-hand side terms are taken at α = α max . A fourth condition is required by the fact that, contrary to the elastica problem where the length of the beam is fixed, we must here determine self-consistently the mesh length rising above the ground. To close the system, we therefore impose the usual mobile contact-point condition Γ(α = 0) = 0 that assumes that there is no adhesion energy between the mesh and the ground [START_REF] Burridge | Peeling, slipping and crackingsome one-dimensional free-boundary problems in mechanics[END_REF].

On the one hand, the forces can be calculated from the momentum transfer between the lattice and the expelled sticks and we have F ⊥ = µv 0 v sin ψ and F = µv 0 (v 0v cos ψ).

On the other hand, the torque exerted at the free end of the cobra can be neglected. Indeed, assuming that all the elastic energy is converted into rotational energy of the sticks, we have the upper bound C max = γ √ 2IE el where I is the moment of inertia of a stick. In Eq. ( 5) and ( 6), the force must be compared to C 2 / K. We have

C 2 KF γ 2 I(Ke 2 /L 3 )) KγM v 0 e L 2 1,
where K K is the flexion modulus of a single stick [19] and we have used the fact that v v 0 γL and I M L 2 . We thus see that for thin sticks, the torque does not affect much the shape of the cobra.

Eq. ( 4) can be solved numerically in the general case using the shooting method and the height of the cobra can be obtained from

H = αmax 0 sin α Γ(α) dα. (8)
The analysis of equations [START_REF] Kolmogorov | Study of the diffusion equation with growth of the quantity of matter and its applications to a biological problem[END_REF](5)[START_REF] Ya | The mathematical theory of combustion and explosions[END_REF][START_REF] Clavin | Combustion waves and fronts in flows[END_REF][START_REF] Stronge | The domino effect: a wave of destabilizing collisions in a periodic array[END_REF] shows that H follows the general scaling When gravity can be neglected, Eq. ( 4) can be solved analytically and yields Γ = Γ 0 cos(ψα) -cos(ψ) with Γ 0 = 2µvv 0 / K. We then obtain α max = 2ψ and the dimensionless height can be expressed in terms of the elliptic integrals E and K [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF] with

H = K µv 0 v h ψ Λ = g K µv 3 0 v 3 . ( 9
) • • • • • • • • • • • • 0.0 0.
h ψ (0) = 2 sin(ψ) [2E(sin(ψ/2)) -K(sin(ψ/2))] . (10)
For large values of Λ, gravity becomes dominant and the cobra does not rise as high. In this regime, α max → 0 and we can therefore neglect the lower order derivatives in each sides of Eq. ( 4) leading to the simplified expression

d 3 Γ 2 dα 3 = 2µg KΓ 2 dΓ dα . (11)
This equation can be solved analytically leading to an asymptotic behaviour h ψ 2 sin 4 ψ/3Λ 3 . The asymptotic behaviours obtained in both the weak and strong gravity regimes can be understood by a straightforward argument. We note first that Eq. ( 7) and ( 8) lead to the following scalings

F ⊥ K Γ 2 0 α max H α 2 max Γ 0 . ( 12
)
We can then distinguish two regimes. For small g, the height is saturated and

α max 1, hence Γ 0 F ⊥ / K and H Γ -1 0 K/F ⊥ .
The scaling for H yields the condition h ψ 1 for weak gravity. Using Eq. ( 1), we can express v 0 and F ⊥ with K, e and L. We then obtain a simple scaling for H L 2 /e, which does not FIG. 6: Set of parameters (e, L) for which a cobrastick wave can be observed (red region). This region is delimited by the two conditions given in Eq. ( 14) and ( 16) : the limit set by gravity Lmax = 36Ee 4 /ρg 1/5 and the breaking limit

Lmin = 1.5 × 3 E/σ * e.
Filled blue dots : sets of parameters for which the cobrastick wave is observed. Open blue dots : sets of parameters for which the cobrastick wave could not be observed (the sticks are too small and therefore break). (i) Schematics of the shape of a stick in the lattice with the most probable breaking region.

depend any more on the elasticity of the mesh. This purely geometric scaling stems from the fact that, when gravity is negligible, stick elasticity provides both the thrusting and restoring forces responsible for the shape of the mesh. For heavy sticks, the lattice is almost horizontal and the value of α max is set by the balance between F ⊥ and the weight. The length of the cobra being s max α max /Γ 0 , we have thus the additional condition

F ⊥ µgα max Γ 0 . (13)
Combining Eq. ( 12) and ( 13) yields the condition α max Λ -2 and h ψ 1/Λ 3 . The transition between the two regimes occurs for Λ 1.

We now compare the previous model to our measurements. We measured the velocity v and the angle ψ at which the sticks are expelled. We observe that for almost all stick models, ψ varies between 50 o and 70 o . In Fig. 5, we compare our measurements to the predicted value h ψ=60 o without any adjustable parameter. Except in the weak-gravity regime, we observe a relatively good agreement between experiment and theory. We attribute the saturation of the height of the cobra-wave for small Λ to the strong curvature of the mesh (in this regime the radius of curvature is only a few times larger than stick length), leading to a breakdown of the underlying assumptions of the theoretical model. For instance the validity of the continuum approximation for the description of the mesh, or the linear approximation for the bending energy. Friction can also play a larger role, and the strong deformation can weaken the structure, preventing it from reaching its predicted height.

Finally, we discuss the condition of existence of the cobra-wave. The first requirement is that the curvature energy stored in a single stick (E el = 18Ewe 5 /L 3 ) should overcome the gravitational energy (E g = ρgweL 2 /2). This leads to an upper bound for the length L of the sticks:

L < L max = 36Ee 4 ρg 1/5 . ( 14
)
However, the length L of the sticks cannot be too small because if so it becomes impossible to build the lattice: the sticks either break or slide over each other destroying the lattice. The breaking condition is derived from a simple scaling law for the bending stress in a beam which sets an upper limit for the curvature of a stick in the lattice (Fig. 6.i) :

C ∼ e d 2 < σ * Ee , (15)
where the length d = (Lw)/3 is defined in Fig. 6.i and σ * is the bending stress at rupture of the material. We then get a lower bound for the length L of the sticks

L > L min b ∼ 3 E σ * e + w. (16)
For wooden sticks, these two conditions set the boundaries of the cobra-wave region of existence. The phase diagram (e, L) is plotted in Fig. 6 with the region of existence of the cobrastick wave in red, assuming the width w to be negligible compared to the length L of the sticks. In particular, this region is semi-infinite, the maximal length being L * ∼ 10 m which corresponds to a maximal thickness e * ∼ 25 cm. We identify a third condition on L which is related to the friction of the sticks and can be obtained from Coulomb's law of friction

L > L min f ∼ 3e f + 4w, ( 17
)
where f is the coefficient of static friction. It is noticeable that the relative position of the breaking threshold and the friction threshold depends on the material of the sticks. For instance, for sticks made of PVC, unlike wooden sticks, the friction threshold is higher than the breaking threshold and is therefore the one that sets the lower bound for the length of the sticks.

In conclusion, we have shown that the shape of the popsicle-stick cobra wave was the result of a competition between on the one hand elastic and gravitational restoring forces and on the other hand the thrust provided by the expulsion of the sticks. Depending on the relative importance of gravity, we identified two asymptotic regimes. In particular, for negligible gravity, the cobra rises at a height which is solely set by the weaving pattern and the dimensions of single sticks. Finally, we showed that the Cobra wave can only exist in a narrow region of the parameter space bounded by gravity and rupture of the sticks. In future work, we will study more carefully the local expulsion dynamics of individual sticks to clarify the role of the weaving pattern.

Supplementary materials

I. CHARACTERISTICS OF THE STICKS

We measured the Young modulus E of our sticks through measurements of the load on a force sensor for given deflections of the sticks and we found E 15 GPa. We estimated the bending stress at rupture of the sticks σ * 120 MPa experimentally through a three-point flexural test.

Sticks L (mm) w (mm) e (mm) M (g) for which the cobrastick wave was not observed and appear in Fig. 6.

II. ELASTIC MODULUS OF THE MESH

The elastic properties of the mesh are described by a bending energy

E el = K 2 Γ(s) 2 ds, ( 18
)
where Γ is the local curvature of the mesh and K the effective bending modulus. K is determined experimentally by measuring the load on a force sensor for given deflections of the mesh. In Fig. 7, we plot the measured value K normalized by the bending modulus of a single stick K = EI ⊥ (where I ⊥ = we 3 /12 is the second moment of area) as a function of the weaving angle θ. For θ = π/4, we find in particular K 0.7K. Assuming that the sticks are simply wrapped with an angle θ on a cylinder of radius Γ -1 , one obtains K = 3K cos(θ) 3 . This scaling is close to the fitted behaviour K/K = k 0 cos(θ) 3 , with k 0 = 2.2. The discrepancy with the expected value k 0 = 3 may come from the torsion of the sticks in the mesh or from the fact that even in the absence of macroscopic bending of the mesh, individual sticks are already bent by the weaving of the mesh.

FIG. 1 :

 1 FIG. 1: (a) Picture of a wooden stick with its characteristic parameters : length L, width w, thickness e, mass M , density ρ and Young's modulus E. (b) Schematics of the lattice with definition of the angle θ of the lattice and the period a of the pattern. The blue and red sticks are the sticks that end the lattice. The construction of the lattice starts with the blue stick (number 1), then the sticks are added one after the other (according to the numbering for the first four sticks), with an alternation of orange and green sticks up to the final red stick.

FIG. 2 :

 2 FIG. 2: Time-lapse photographies of the cobra wave obtained with sticks of type-1 (see [19]) for θ = 45 o , from videos taken at 1000 fps with a Photron-Fastcam high-speed camera. The different colors represent the wave at different instants. (a) At the beginning (∆t = 100 ms between two consecutive images); (b) During the stationary phase (∆t = 70 ms), with v0 2.2 m/s the velocity of the wave front ; (c) At the end (∆t = 70 ms). Scale bars are 10 cm long.

FIG. 3 :OFIG. 4 :

 34 FIG. 3: Speed of the wave-front v0 (filled dots) and speed of the expelled sticks v (open dots) as a function of a characteristic speed c(e/L) 2 for θ = 45 o and six different kinds of sticks (see [19]). The black line corresponds to the fit v0 = 3.95 c(e/L) 2 and the dashed line to the fit v = 5.46 c(e/L) 2 . (i) Speed of the cobra wave v0 as a function of the angle of the lattice θ for sticks of type-1.

FIG. 5 :

 5 FIG. 5: Dimensionless height h ψ as a function of the dimensionless number Λ = g K/ (µv 3 0 v 3 ) which compares gravity and elasticity, for the six different kinds of sticks ([19]). The solid line corresponds to the prediction of Eq. (4-8) for ψ = 60 o . The shaded band corresponds to the observed 10 o variations of the ejection angle. The red dashed line represents the large Λ expansion h ψ 2 sin 4 ψ/3Λ 3 for ψ = 60 o .

  k in g th re s h o ld L im it s e t b y g ra v it y

FIG. 7 :

 7 FIG.7: Normalized bending modulus of the mesh K/K as a function of the weaving angle θ obtained experimentally by measuring the force on a load cell for imposed deflections of the mesh. Solid line: fit of the data using the theoretical function K/K = k0 cos(θ) 3 with k0 = 2.2, close to the expected value k0 = 3 derived from a purely geometric argument.

TABLE I :

 I Geometrical and mechanical properties of the wooden sticks used in this study. Sticks 7 and 8 are sticks

Acknowledgments

This work was inspired by one the problems of the 2016 edition of the International Physicists' Tournament (IPT). The authors acknowledge fruitful discussions with Basile Audoly, Daniel Suchet and the ENS and École Polytechnique IPT teams.


Bijan Mohammadi ---- Montpellier

HEURISTIC GLOBAL OPTIMIZATION AND IDENTIFICATION OF GENERAL PARETO FRONTS

Keywords: Heuristic optimization, Multidisciplinary, multicriteria, Optimization, data fusion, Pareto front

We recall some of the multi-criteria and multi-disciplinary optimization formulations with emphasis on algorithmic and implementation issues. We use a partial controllability result which states that the vicinity of any point in the optimization space can be reached in finite time from any other point in the space following a trajectory solution of a particular second order dynamical system. We use this continuous point of view to propose a class of heuristic optimization algorithms also addressing global search issues. For multi-criteria problems this means that any point on a Pareto front can be reached even if the front is non-convex as far as the optimization algorithm has inertial features. The paper also recalls the conditions under which a Pareto front can be accessed by descent methods applied to multicriteria problems using the weighted sum approach. The question of the density of the distribution of the points over the front is also addressed. The approach is illustrated on academic examples as well as on a multi-criteria distributed inverse problem for a heterogeneous data fusion problem.

Introduction

The paper is dedicated to two problems of interest in optimization. First, we address the question of how descent methods can be improved to perform global search ? Once this question discussed, one sees why global search methods are necessary to reach points on general non-convex Pareto fronts where descent methods fail. Addressing these questions lead to the development of a class of heuristic global optimization methods as discretization of special 2 nd order continuous dynamical systems.

Descent algorithms can be interpreted as discrete forms of first order systems [START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF]: a descent method solves a Cauchy problem for such a system. Now, suppose we know the infimum of the functional and the question of global optimization is then where it is reached ? Global optimization becomes then a boundary value problem with an initial condition in the admissible space and a final condition on the functional. Therefore, solving a global optimization problem with such methods leads to an over-determined problem. In other words, if the final condition is given, the initial condition should be prescribed and should be a variable of the problem as well. Intuitively, this suggests that a descent method is willing to succeed if it starts somewhere in the optimal attraction basin (i.e. where the infimum is).

The paper discusses how a controllability result can be used in the development of heuristic optimization algorithms using recursive shootings in order to solve a boundary value problem. We recover the philosophy behind the Bock's multiple shooting algorithm [START_REF] Bock | A multiple shooting algorithm for direct solution of optimal control problems[END_REF] and the Lion's parareal method [START_REF] Maday | A parareal in time discretization of pde's[END_REF].

The paper then considers the application of these ideas to multi-criteria problems. We will see that weighted sums used to transform a multi-criteria problem into a mono-criteria one can still be used to recover general Pareto fronts and not only convex ones. This is interesting because this method is widely used and appreciated by engineers.

Finally, we discuss the question of the distribution and localization of the points on the Pareto front. Indeed, we will see that despite general fronts can be reached, the points are often oddly distributed with concentrations in some regions of the front while we might want some kind of uniform distribution of the points on the front. To address this issue, we propose to couple the different trials introducing a particular repulsion force in the corresponding second order dynamical systems. This point is however not central in the paper and the concepts of the paper can be used together with other proposals such as the normal boundary intersection or the normal constraint methods [START_REF] Das | Normal-boundary intersection: A new method for generating pareto optimal points in multicriteria optimization problems[END_REF][START_REF] Ismail-Yahaya | The normalized normal constraint method for generating the pareto frontier[END_REF]. The proposed method has however the advantage to remain compatible with the weighted sum formulation and easy to explain in engineering environments.

Multi-criteria optimization

The literature on multi-criteria optimization is vast [START_REF] Das | A closer look at drawbacks of minimizing weighted sums of objectives for pareto set generation in multicriteria optimization problems[END_REF][START_REF] Das | Normal-boundary intersection: A new method for generating pareto optimal points in multicriteria optimization problems[END_REF][START_REF] Deb | Multi-objective genetic algorithms: Problem difficulties and construction of test problems[END_REF][START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF][START_REF] Pareto | Manuale di Economia Politica[END_REF][START_REF] Stadler | A survey of multicriteria optimization, or the vector maximum problem[END_REF]. Let us briefly recall some notions and take this opportunity to introduce the notations which will be used throughout the paper. We consider n real functions {J i (x), i = 1, .., n}. One way to look at multidisciplinary optimization is minimizing a vector [START_REF] Jahn | Vector Optimization: Theory, Applications and Extensions[END_REF]: [START_REF] Ismail-Yahaya | The normalized normal constraint method for generating the pareto frontier[END_REF] min

x∈O ad ⊂IR N J, J = (J 1 (x), .., J n (x)) t ,
which brings the question of order relation in a vector space which is in turn related to the norm one chooses for the optimization: shall we consider an L 1 , L 2 or ... L ∞ norm ? or shall we consider several or all of them ? There is no single answer to this question. Instead, one should look for trade-off between conflicting objectives and Pareto equilibrium [START_REF] Pareto | Manuale di Economia Politica[END_REF].

2.1. Multicriteria optimization algorithms. Let us recall a few strategies to solve problem (1) which has a wide literature dedicated to [START_REF] Ehrgott | Multicriteria Optimization[END_REF][START_REF] Jahn | Vector Optimization: Theory, Applications and Extensions[END_REF][START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF]. This paper is not on introducing a new algorithm for this purpose. The algorithms can be cast into two classes considering either the functionals separately or through a unique representation, after weighted summation for instance.

A natural representative of the first class is the min-max approach where one looks for minimizing sequences for the following problem starting from an initial value x 0 ∈ O ad :

(2) min

x∈O ad max j=1,...n J j (x) J j (x 0 ) .
To make sure the sequence x p is a descent sequence for each of the vector functional components, one should require a new proposal x p+1 not to degrade any of the current achievements. This means that x p+1 should be solution of the following constrained minimization problem:

min x∈O ad max j=1,...n J j (x) J j (x 0 ) such that C i = J i (x) -J i (x p ) ≤ 0, i = 1, ..., n, i = j.
This can be achieved solving the KKT conditions for the Lagragian L j (x, Λ) = J j + Λ t C with C ∈ IR n-1 the vector of constraints and Λ ∈ IR n-1

+

is the vector of positive Lagrange multipliers. This Lagrangian might change at each iteration. Looking to achieve the KKT conditions means building x p+1 aiming at vanishing ∇ x L j = ∇ x J j + Λ t ∇ x C with the complementary condition Λ t C = 0 making a distinction between an inactive constraint where C i < 0 and λ i = 0 and an active one where C i = 0 and λ i > 0. In this former situation we have λ i = -∂J j /∂C i meaning that the Lagrange multiplier gives the sensitivity of J j with respect to a perturbation of this particular active constraint. Once the optimum is reached, this information permits a better construction of the Pareto front.

2.2. Geometric interpretation. This construction can be understood geometrically. To be admissible, a search direction d needs to be orthogonal to the following convex cone:

(3)

S p = {x | x = p i=1 β i ∇ x C i , β i > 0 | C i = 0} ⊆ IR p ⊂ IR n-1 ,
where we have p active constraints. At the solution, ∇ x J is orthogonal to this cone. Before working on the cone, let us start defining a local orthonormal basis {q i=1,...,p } for Sp defined by (3) but with β i ∈ IR. Sp is a vector space and an orthonormal basis can be defined for Sp orthonormalizing the columns of ∇ x C by the Gram-Schmidt procedure. Now, we consider the vectors q i = ±q i and the sign chosen such that < q i=1,...,p , ∇ x C k >≥ 0, if C k = 0 for k = 1, ..., p (i.e. pointing inside the cone for an active constraint). {q i=1,...,p } are the generators of the cone S p deduced from a basis of Sp . If the generators cannot be defined, the problem is found having no solution as at least two of the constraints are incompatible with their gradients parallel and pointing in opposite directions. These generators permit to define an admissible search direction d at iteration n following:

(4) d = ∇ x J j - p i=1 χ i < q i , ∇ x J j > q i ,
with χ i = 0 if < q i , ∇ x J j >≥ 0 and χ i = 1 if < q i , ∇ x J j >< 0. At the optimum, the Lagrange multipliers for active constraints correspond to 0

< λ i = -< q i , ∇ x J j >.
This construction can also be obtained applying the Uzawa algorithm for the solution of: min

x∈O ad max Λ∈IR n-1 + L j (x, Λ),
at each iteration of problem [START_REF] Tappeta | Ability of objective functions to generate points on non-convex pareto frontiers[END_REF]. The solution of (2) belongs to the Pareto front as there is no way to improve one of the criteria without degrading the others. This approach requires the gradients of all the functionals and, for large dimensional optimization problems, this means the solution of n adjoint problems [START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF] at each iteration which can be very costly, especially with sophisticated state equations solvers. With black-box solvers and if the size of the problem is small finite differences are often used instead of adjoint because of their natural parallelism and simplicity together with a weighted sum formulation to only handle one functional.

2.3.

Weighted sum formulation. The classical method for multi-objective optimization remains the weighted sum method because of its simplicity [START_REF] Das | A closer look at drawbacks of minimizing weighted sums of objectives for pareto set generation in multicriteria optimization problems[END_REF][START_REF] Das | Normal-boundary intersection: A new method for generating pareto optimal points in multicriteria optimization problems[END_REF][START_REF] Stadler | A survey of multicriteria optimization, or the vector maximum problem[END_REF]:

(5) J(x) = Σα i J i (x), Σα i = 1, α i ∈ [0, 1].
A sampling in α i enables to reach various points in the Pareto front using a descent method for J(x). Indeed any point x ∈ O ad minimizing J is a Pareto equilibrium.

Actually any point x ∈ O ad dominating x would yield an improved weighted sum J(x) < J(x) in contradiction with the minimizing property of x.

The converse is true in a convex setting. That is, any Pareto equilibrium x is a minimum point of some weighted sum J given by ( 5), provided that the functionals J i and the set O ad are convex. Indeed set

C = {ξ ∈ IR n , ∃x ∈ O ad , ξ i ≥ J i (x), ∀i ∈ {1, . . . , n} }, K = {ξ ∈ IR n , ξ < J i (x), ∀i ∈ {1, . . . , n} }.
These are convex nonvoid subsets of IR n , and

K is open. Moreover C ∩ K = ∅, since any point ξ ∈ C ∩ K would verify J i (ξ) < J i (x)
in contradiction with x being a Pareto equilibrium. Hence, there exists a hyperplane separating sets C and K (see figure 1); that is, there exist a vector α = (α 1 , . . . , α n ) ∈ IR n \ {0} and a scalar b such that

ξ ∈ K ⇒ α ξ t ≤ b, (6) ξ ∈ C ⇒ α ξ t ≥ b. (7)
Let ξ be some fixed point in K. For all λ ∈ IR n with nonnegative components, ξ -λ lies in K too, and with (6) we have αξ t -αλ t ≤ b. Hence necessarily a i ≥ 0 for all i ∈ {1, . . . , n}, and there is no harm in supposing n i=1 α i = 1.

Figure 1. Hyperplane separating sets C and K.

Let K denote the closure of K in IR n . Since (J 1 (x), . . . , J n (x)) ∈ K ∩ C, we have n i=1 α i J i (x) ≤ b ≤ n i=1 α i J i (x), hence b = n i=1 α i J i (x).
Finally with ( 7) we obtain in particular

x ∈ O ad ⇒ n i=1 α i J i (x) ≥ n i=1 α i J i (x),
which means that x is a minimum point of J.

A shorter proof may be given under reinforced assumptions. If x is a Pareto equilibrium, it is a solution to the constrained optimization problem:

inf {J 1 (x), x subject to J i (x) ≤ J i (x), i ∈ {2, . . . , n}}. Indeed J i (x) ≤ J i (x), for i ∈ {1, . . . , n}, entails J i (x) = J i (x) for i ∈ {1, . . . , n}.
Under the extra assumption (see [START_REF] Michaille | Variational Analysis in Sobolev and BV Spaces: Applications to PDEs and Optimization[END_REF] e.g.)

i∈{2,...,n} {x ∈ O ad , J i (x) < J i (x)} = ∅,
there exist nonnegative multipliers α i for i ∈ {2, . . . , n} such that x minimizes the functional J 1 + i∈{2,...,n} α i (J i -J i (x)), and the functional J 1 + i∈{2,...,n} α i J i as well. Now divide the latter by 1 + i∈{2,...,n} α i to recover a weighted sum J such as given by [START_REF] Bock | A multiple shooting algorithm for direct solution of optimal control problems[END_REF].

It is reported that this approach does not find Pareto optimal solutions in nonconvex regions of the Pareto front [START_REF] Tappeta | Ability of objective functions to generate points on non-convex pareto frontiers[END_REF]. This is linked to the fact that when the tangent to the Pareto front is parallel to one of the axes the corresponding functional has a plateau which cannot be crossed by a descent method not having inertial features. Now consider the situation of two functionals and a Pareto front P (J 1 , J 2 ) = 0. In cases, where the front is non-convex but the tangent never parallel to one of the axes the implicit function theorem can be applied and one can find J 1 = j 1 (J 2 ), J 2 = j 2 (J 1 ). In these cases, as shown in [START_REF] Tappeta | Ability of objective functions to generate points on non-convex pareto frontiers[END_REF], one can still use a descent method for increasing power m in the functional:

J(x) = αJ m 1 (x) + (1 -α)J m 2 (x), α ∈ [0, 1].
The aim is to find the lowest m making the Pareto front P (J m 1 , J m 2 ) convex. But, this approach has the disadvantage of degrading the conditioning of the optimization problem.

In the sequel we show that the failure of the weighted sum approach is due to a lack in global minimization feature in descent methods and propose an alternative solution which is not submitted to the mentioned limitation on the tangent to the Pareto front.

Controllability and global optimization

Consider the minimization of a functional J(x) ∈ IR, x ∈ O ad . We suppose the problem admissible (i.e. there exist at least one solution x m to the problem: J(x m ) = J m ). Most minimization algorithms can be seen as discrete expressions for the following Cauchy problem [START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF]:

(8) M (x(t))x t = -d(x(t)), x(t = 0) = x 0 .
M is aimed to be positive definite and M -1 d is built to be an admissible direction. Global solution of (8) means, for instance, finding x(T ) for finite T such that

(9) M (x(t))x t = -d(x(t)), x(0) = x 0 , J(x(T )) = 0.
This is an over-determined boundary value problem and it tells us why one should not solve global optimization problems with methods which are discrete form of ( 8). This limitation can be somehow removed if one can provide an initial condition in the attraction basin of the global optimum. To remove this overdetermination one shall consider, for instance, the following second order system:

(10) ηx tt + M (x(t))x t = -d(x(t)), x(t = 0) = x 0 , J(x(T )) = J m .
This can be solved using solution techniques for boundary value problems with free surface. An analogy can be given with the problem of finding the interface between water and ice which is only implicitly known through the iso-value of zero temperature [START_REF] Mohammadi | semi-deterministic global optimization method and application to the control of burgers equation[END_REF].

3.1. Practical impacts. In practice the final condition is not aimed at being exactly realized. Rather, one would like, for a given precision δ in the functional, to build at least one trajectory (x(t), 0 ≤ t ≤ T δ ) passing for finite T δ in the ball B δ (x m ). This can be summarized as:

(11) ∀ δ > 0, ∃ (v, T δ ) ∈ O ad × [0, +∞[ such that J(x v (T δ )) -J m ≤ δ.
If J m is unknown T δ defines the maximum calculation complexity wanted. This is how optimization algorithms are actually used: one chooses the maximum calculation/memory effort one can afford and one looks for the best solution for this complexity. In these cases, setting J m = -∞, one retains the best solution obtained over [0, T δ ]. In other words, one solves:

(12) ∀ (δ, T δ ) ∈ IR + × [0, +∞[, ∃ (v, τ ) ∈ O ad × [0, T δ ], such that J(x v (τ )) -J m ≤ δ.
The following theorem formalizes this controllability question with M = Id and d = ∇J.

Theorem 1: Let J : IR N → IR be a C 2 -function such that min IR N J exists and is reached at x m ∈ IR N . Then for every (x 0 , δ) ∈ IR N × IR + , there exists (σ, t ) ∈ IR N × IR + such that the solution of the following system: [START_REF] Michaille | Variational Analysis in Sobolev and BV Spaces: Applications to PDEs and Optimization[END_REF] ηx

tt (t) + x t (t) = -∇J(x(t)), t ≥ 0, x(0) = x 0 , x t (0) = σ,
with η ∈ IR, passes at time t = t into the ball B δ (x m ).

Proof :

We assume x 0 = x m (x 0 = x m is a trivial case). Let ε > 0, we consider first the initial value problem: [START_REF] Ivorra | Optimization strategies in credit portfolio management[END_REF] ηy tt (t) + εy t (t) = -ε 2 ∇J(y(t)), t ≥ 0, y(0) = x 0 , y t (0) = (x m -x 0 ), with ∈ IR + \{0}.

Let us show that y(.) passes at some time into the ball B δ (x m ):

• If ε = 0, we obtain the following system :

(15) ηy tt,0 (t) = 0, t ≥ 0, y 0 (0) = x 0 , y t,0 (0) = (x m -x 0 ). System (15) describes a straight line of origin x 0 and passing at some time θ by the point x m , i.e. y 0 (θ ) = x m . • If ε = 0, system (14) can be rewritten as (w 2 ) t = f (η, w, ε), with w = (w 1 , w 2 ) t = (y(t), ηy t (t)) t and f = -ε/ηw 2 -ε∇J(w 1 (t)) satisfying the Cauchy-Lipschitz conditions (J is a C 2 function). Then lim ε→0 |y ε (θ )y 0 (θ )| = 0. Thus, for every δ ∈ IR + \{0}, there exists ε δ such that for every ε < ε δ :

(16) |y ε (θ ) -x m | < δ.
Finally, let us consider the change of variables t = ε δ t and x(t

) = y ε δ ( t ε δ
). Then system (14) becomes:

(17) ηx t t (t ) + x t (t ) = -∇J(x(t )), t ≥ 0, x(0) = x 0 , ẋ(0) = ε δ (x m -x 0 ).
Let ϑ = ε δ θ . Under this assumption, x(ϑ) = y ε δ (θ ). Thus, due to ( 16): |x(ϑ)x m | < δ. We have found σ = ε δ (x m -x 0 ) ∈ IR N and t b = ϑ ∈ IR + such that the solution of system (13) passes at time t b into the ball B δ (x m ). 2

3.2. Algorithmic considerations. This theorem is insightful and permits to derive a global minimization algorithm which requires δ > 0, the infimum J m and 0 ≤ T δ < ∞ [START_REF] Mohammadi | A multi-layer line search method to improve the initialization of optimization algorithms[END_REF]. If J m is unknown, one sets J m = -∞. Again, T δ shall be seen as a maximal admissible complexity for the calculation. We then minimize the functional h δ,T δ ,Jm : O ad → IR + :

(18) h δ,T δ ,Jm (v) = min A (J(x v (τ )) -J m ),
where

A = {x v (τ ) ∈ O ad , τ ∈ [0, T δ ]}.
Hence, global minimization becomes a nested minimization problem where one looks to improve the initial impulse v for a 'core' minimization sub-problem to generate the trajectory (x v (τ ), 0 ≤ τ ≤ T δ ) by:

(19) ηx tt + M (x(t))x t = -d(x(t)), x(t = 0) = x 0 , x t (t = 0) = v, t ∈ [0, T δ ].
In cases one would like to use a core minimization algorithm which is a discrete form of a first order dynamical system (i.e. η = 0), [START_REF] Medic | Nsike unstructured solver for laminar and turbulent incompressible flows simulation[END_REF] becomes:

(20) M (x(t))x t = -d(x(t)), x(t = 0) = v, t ∈ [0, T δ ].
In other words, the external minimization problem aims at improving the initial condition for the core optimization sub-problem. Both the initial condition and initial impulse can therefore be considered as optimization variables by the external optimization problem. To solve the external minimization problem, our approach is based on constructing a multi-level reduced-order parametric model by 'data assimilation' or 'learning' [START_REF] Mohammadi | A multi-layer line search method to improve the initialization of optimization algorithms[END_REF][START_REF] Mohammadi | semi-deterministic global optimization method and application to the control of burgers equation[END_REF][START_REF] Ivorra | Optimization strategies in credit portfolio management[END_REF][START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF]. This algorithm has shown great efficiency for cases where the functional is coercive and comes from perturbation of a convex functional [START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF]. However, we keep its discussion to a minimum as it is not central to our topic here.

Application to the identification of Pareto fronts

We saw how to link any couple of points in the admissible space in finite time with some given accuracy by a trajectory solution of a second order dynamical system. If the problem is admissible a point in a Pareto front corresponds to a point in the admissible control space. Theorem 1 states that the vicinity of this point can be reached in finite time from any other point in the control space following a trajectory solution of the second order system (13). The algorithm above can therefore be used to recover any point in a Pareto front with:

(21) h δ,T, Jm (v) = min A J(x v (τ )) -J m ,
regardless of whether the front is convex or not. The unknown components of J m are set to -∞ as in mono-objective minimization. One looks for the best solution at given complexity. In all cases, this means that one shall use a minimization method which has global search features to cross local minima in the case of non-convex Pareto fronts. Let us illustrate this construction on a simple example.

4.1.

A first illustration. We consider the functionals

J 1 = x 1 + x 2 and J 2 = (1/x 1 + x 2 ) + β(exp(-100(x 1 -0.3) 2 ) + exp(-100(x 1 -0.6) 2 )) for x ∈ [0.1, 1] 2 .
Functional J 1 is bilinear and J 2 has local minima for non zero β. We consider the cases of β = 0 and β = 3. One expects a local minimization algorithm to be captured by local minima for the latter situation showing that the corresponding Pareto front generation requires a global minimization algorithm.

Figure 2 shows examples of convex and non-convex and non-connected Pareto fronts for these two functionals when the parameter space has been uniformly sampled. The figure also shows the front identified by a steepest descent method with optimal step size and when enriching the method with second order dynamics features as presented in the paper. This improves the global search capabilities of the method and permits to reach the non-convex region in the front. We have solved 20 minimization problems based on uniform sampling of the weighting parameter α ∈ [0, 1] in the weighted sum with the steepest descent method which corresponds to a discrete expression of [START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF]. We have then solved the same minimization problems with a discrete expression of [START_REF] Medic | Nsike unstructured solver for laminar and turbulent incompressible flows simulation[END_REF] with M = Id, d = α∇J 1 + (1 -α)∇J 2 and η = 1.

4.2.

Improving the sampling of the Pareto front coupling dynamical systems . We saw that the vicinity of any admissible point in the Pareto front can be reached in finite time from any point in the admissible domain. As we are interested in a discretization of the Pareto front, this means we proceed with several minimization in parallel as described in previous sections. However, these minimization by our dynamical systems can interact in order to increase their global search features. This coupling can also be used to achieve a more uniform distribution of the points over the Pareto front to avoid the kind of situation one sees in figure 2 even for a convex Pareto front where points accumulate in some regions regardless of whether the dynamics is first or second order. A natural way to introduce the uniformity of the distribution is to ask for a point to be pushed away the positions already occupied by others. This can be done introducing a repulsive force in our dynamics. This couples the different paths using cross-information and involving a global gradient. [START_REF] Pareto | Manuale di Economia Politica[END_REF] to be compared to lower images in figure 2.

More precisely, we consider q individuals (x j , j = 1, ..., q) with their motions prescribed by q systems (19): [START_REF] Pareto | Manuale di Economia Politica[END_REF] ηx j tt + M x j t = d j + G j , x j (t = 0) = x 0 , x t (t = 0) = v j . d j is based on the gradient of the functional J j given by (5) for a particular choice of the n-uplet α j ∈ [0, 1] n such that n i=1 α j i = 1. G j is a global gradient representing the interaction between the individuals:

G j = - q k=1,k =j J j -J k ||x j -x k || 2 (x j -x k ), for j = 1, ..., q.
This repulsive force aims at forcing the search to take place in regions not already visited. The right-hand side in [START_REF] Pareto | Manuale di Economia Politica[END_REF] is therefore a balance between local and global forces.

This construction aims at achieving similar behavior than methods such as the normal boundary intersection [START_REF] Das | Normal-boundary intersection: A new method for generating pareto optimal points in multicriteria optimization problems[END_REF] or the normal constraint [START_REF] Ismail-Yahaya | The normalized normal constraint method for generating the pareto frontier[END_REF]. It has the great advantage to be easy to implement and to be naturally parallel. Figure 3 shows the impact of this correction on the Pareto fronts given in figure 2. In the sequel, we will apply this combination to a more sophisticated multi-criteria and multi-physics data fusion problem.

Application to multi-criteria and multi-physics data fusion

As previously, our multi-criteria problem involves a vector of functionals:

J = (J 1 (u 1 (x)), ..., J n (u n (x))) t ∈ IR n ,
for a same set of optimization parameters x ∈ O ad ⊂ IR N . For simplicity we assume all the parameters are public and do not split the optimization parameters between public or global and private or local or disciplinary sets. One natural way to achieve this splitting in our context is to consider territory splitting strategies and to affect or not a given parameter x i=1,...,N to a given discipline j represented by J j=1,...,n following its sensitivity dJ j /dx i . The difficulty is that these sensitivities are a propri unknown. This splitting is an important step in the formulation of a multidisciplinary optimization problem but it is not central in our discussion [START_REF] Desideri | Cooperation and competition in multidisciplinary optimization[END_REF][START_REF] Ph | The topological asymptotic for pde systems: The elasticity[END_REF][START_REF] Martins | Multidisciplinary design optimization: A survey of architectures[END_REF]. 5.1. Multicriteria and multiphysics data fusion. Our interest is on data assimilation/fusion and m represents the different data set to be assimilated. In this case, J i reads for instance as:

(23) J i (u i (x), u * i ) = u i (x) -u * i 2 Σ -1 i + T i (x), with T i (x) ≥ 0.
T i (x) represents regularity control requirements or additional constraint. The norm . Σ -1 i uses the inverse of the covariance matrix Σ i of the different target observations

u * i .
The problem is multi-physics as getting the functional J i necessitates the solution of a different state equations F i (u i (x)) = 0 for the description of the different state variable u i (x) ∈ IR M . One expects M to be very large. For instance, in the case of a distributed optimization problem involving the solution of PDE models, M typically represents the size of the mesh. 5.2. An example of vector and scalar fields assimilation. Consider the general problem of mixing a passive scalar (here the temperature) in a fluid. The question is to find the flow field which permits to realize a given distribution for a transported quantity with constraint on what the admissible field is. Inversely, one wants to answer the question of source identification from a few flow field and advected quantity observations. The state equations we consider are the incompressible Navier-Stokes equations together and an advection and diffusion equation for the temperature. Details on these equations and the flow solver can be found in [START_REF] Medic | Nsike unstructured solver for laminar and turbulent incompressible flows simulation[END_REF]. The flow field can be controlled by various means as, for instance, injection devices or variations of a magnetic field. The former is taken into account through a modification of the boundary conditions for the velocity and the later by the introduction of a source term in the flow equations. Here we consider a vector x of 12 control parameters with lower and upper bounds on each component. The bounds express the admissibility domain for flow directions and maximum possible velocity amplitude. Figure 4 shows an example of flow velocity intensity and temperature distributions for a centrifugal device in which control of such quantities might of interest for suitable functioning.

The two cost functions we consider are non dimensionalized and are of the form:

J 1 (x) = 1 |Ω|T 2 des Ω (T (x) -T des ) 2 , J 2 (x) = 1 |Ω| u des 2 Ω u(x) -u des 2 ,
where upper bars indicate mean values over the domain. To simplify the presentation we have assumed identity covariance matrices in [START_REF] Sa | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF]. These aim at realizing target flow velocity and temperature distributions in the domain or some part of it.

A sampling of the control space to generate the domain of possible (J 1 (x), J 2 (x)) ∈ IR 2 is too costly even if the use of sparse grids [START_REF] Bungartz | Sparse grids[END_REF][START_REF] Sa | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF] or other intelligent sampling techniques can bring some relief. Indeed, sampling an admissible domain in dimension n would require 2 n evaluations of the functional to have just the bounds for the variations of the cost function due to each of the parameters. Here, with only four points over each of the variables range we need over than 16 millions functional evaluations. In this problem we obviously have no idea of how the Pareto front is and we can only guess the front from the simulations. Therefore, what we would like to mainly see through this example is the exploration capability of optimization methods based on discrete forms of either first [START_REF] Mohammadi | Applied Shape Optimization for Fluids (2nd Edition)[END_REF] or second order (19) dynamical systems together with the repulsion forces given in section 4.2.

Figures 5 and6 show visited (J 1 (x), J 2 (x)) points by first and second order dynamics for a same set of optimization problems corresponding to a 20 points uniform sampling of α ∈ [0, 1] and the weighted functional αJ 1 + (1 -α)J 2 as described in the previous academic example. However, this problem is much more computer intensive.

One can see that second order dynamics have higher exploration capabilities. The functional J 1 and J 2 are apparently non-convex as the first order dynamics do not necessarily reach points on the non-convex part of the front as the trajectories are probably captured by local minima, away from the front. We can also see that it is interesting to combine the information from the first and second order constructions. Indeed, some points on the Pareto front (highlighted by the circle in figure 5) have been identified by the first order dynamics while missed by the second order construction. We do not have an explanation for this behavior yet. Combining the two sets of information, we have drawn an estimation of the Pareto front which we think is non-convex and also non-connected.

Concluding remarks

Different formulations of multicriteria optimization problems have been presented together with a class of metaheuristic optimization algorithms intuited from theoretical behavior of continuous dynamical systems. In that purpose, global optimization has been formulated as solution of boundary value problems for second order systems. Two theoretical results have been given for the continuous systems on sub-optimal partial controllability in finite time and on the conditions under which a Pareto front can be accessed by descent methods applied to multicriteria problems using the weighted sum approach. It has been shown that global search features are suitable for the identification of points on general Pareto fronts. The behavior of minimization algorithms based on first and second order dynamical systems have been discussed on an academic example and a distributed data fusion problem in large dimension. For this latter case a sampling of the parameter space is clearly impossible. Both situations feature non-convex and non-connected fronts. It has been shown that minimization algorithms having inertial features permit to reach non-convex regions in Pareto fronts thanks to their higher global search capabilities. The question of the distribution and localization of the points on the Pareto front has also been discussed to see how to monitor the distribution of the points on the front. The different trials have been coupled introducing a particular repulsion force in the corresponding second order dynamical systems. However, to address this issue, the concepts of the paper can also be applied together with the normal boundary intersection [START_REF] Das | Normal-boundary intersection: A new method for generating pareto optimal points in multicriteria optimization problems[END_REF] or the normal constraint [START_REF] Ismail-Yahaya | The normalized normal constraint method for generating the pareto frontier[END_REF] method.

Figure 2 .

 2 Figure 2. Top: Pareto front for (J 1 , J 2 ) given in section 4.1 for β = 0 and β = 3 after uniform sampling of the parameter space. Solutions by optimization algorithms issued from discretization of a first (middle) and second-order (bottom) dynamical system for a uniform 20 point sampling of α ∈ [0, 1].

Figure 3 .

 3 Figure 3. Pareto front by second-order dynamics accounting for repulsive forces (22) to be compared to lower images in figure 2.

Figure 4 .

 4 Figure 4. Typical flow field u distribution and temperature T isocontours for a set of control parameters x. Two injection locations can be seen on the external wall.

Figure 5 .

 5 Figure 5. Visited (J 1 (x), J 2 (x)) points by first order dynamics together with repulsion forces and estimation (continuous line) of the Pareto front combining information from these simulations and those from the second order dynamics shown in figure 6. The first order dynamics has difficulty in reaching the Pareto front as most of the trials seem to be trapped by local minima.

Figure 6 .

 6 Figure 6. Visited (J 1 (x), J 2 (x)) points by second order dynamics together with repulsion forces and estimation (continuous line) of the Pareto front combining information from from first and second orders dynamics.