P Wenger

D Chablat

KINETOSTATIC ANALYSIS AND SOLUTION CLASSIFICATION OF A CLASS OF PLANAR TENSEGRITY MECHANISMS

Keywords: Tensegrity mechanism, kinetostatic model, geometric design, algebraic computation 1

Tensegrity mechanisms are composed of rigid and tensile parts that are in equilibrium. They are interesting alternative designs for some applications, such as modelling musculo-skeleton systems. Tensegrity mechanisms are more difficult to analyze than classical mechanisms as the static equilibrium conditions that must be satisfied generally result in complex equations. A class of planar one-degree-of-freedom tensegrity mechanisms with three linear springs is analyzed in detail for the sake of systematic solution classifications. The kinetostatic equations are derived and solved under several loading and geometric conditions. It is shown that these mechanisms exhibit up to six equilibrium configurations, of which one or two are stable, depending on the geometric and loading conditions. Discriminant varieties and cylindrical algebraic decomposition combined with Groebner base elimination are used to classify solutions as a function of the geometric, loading and actuator input parameters.

INTRODUCTION

A tensegrity structure is an assembly of compressive elements (struts) and tensile elements (cable, springs) held together in equilibrium [START_REF] Fuller | Tensile-integrity structures[END_REF][START_REF] Skelton | Tensegrity Systems[END_REF][START_REF] Motro | Tensegrity systems: the state of the art[END_REF].Their inherent interesting features (low inertia, natural compliance and deployability) make them suitable in several applications. They can also be used as preliminary models in musculo-skeleton systems to analyze animal and human movements [START_REF] Vogel | Cats Paws and Catapults: Mechanical Worlds of Nature and People[END_REF][START_REF] Levin | The tensegrity-truss as a model for spinal mechanics: biotensegrity[END_REF]. A spine can be modelled by stacking a number of suitable tensegrity modules, see for example [START_REF] Sabelhaus | Mechanism design and simulation of the ULTRA spine: a tensegrity robot[END_REF]. Figure 1 illustrates some configurations of a serial assembly of planar tensegrity mechanisms similar to the ones studied in this paper, suitable for a plane model of a spine or of an animal neck. The frame of this work is a preliminary step of a large collaborative project with the Museum National d'Histoire Naturelle, called "AVINECK: the neck of birds, an arm for robots". This project aims at understanding the functional behavior of bird necks and using them as bio-inspired models for the design of innovative robotic arms. In this project, several bird neck models will be considered, starting with a planar model made of tensegrity mechanisms. A first step of the project is to seek for appropriate simple models that lend themselves to symbolic calculations in order to better understand their kinetostatic behavior and to ease further optimizations for design purposes.

A tensegrity mechanism is obtained when one or several elements are actuated. Most research efforts on tensegrity mechanisms have been conducted recently (see for example [START_REF] Sabelhaus | Mechanism design and simulation of the ULTRA spine: a tensegrity robot[END_REF][START_REF] Aldrich | Time-energy optimal control of hyper-actuated mechanical systems with geometric path constraints[END_REF][START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-dof tensegrity mechanism[END_REF][START_REF] Crane | Kinematic analysis of a planar tensegrity mechanism with pres-stressed springs[END_REF][START_REF] Boehler | Definition and computation of tensegrity mechanism workspace[END_REF] and references therein). The input and output equations of a tensegrity mechanism are derived after solving the equilibrium conditions. These conditions are generally obtained by minimizing the potential energy. Symbolic derivation of minimization conditions often generates complex equations and inequations that are difficult to solve. Due to their relative simplicity, planar tensegrity mechanisms are more suitable for algebraic computations. A planar 2-DOF planar tensegrity mechanism was analyzed by Arsenault in [START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-dof tensegrity mechanism[END_REF], where kinetostatics and dynamics models as well as some simplified workspaces were studied. Recently, Boehler [START_REF] Boehler | Definition and computation of tensegrity mechanism workspace[END_REF] proposed a more complete definition of the workspace of that 2-DOF mechanism as well as a higher-order continuation based method to evaluate it.

This paper goes beyond a first study presented in [START_REF] Wenger | Kinetostatic Analysis and Solution Classification of a Planar Tensegrity Mechanism[END_REF] and analyses in details a family of one-degree-of-freedom planar tensegrity mechanisms made of one base telescopic rod, two crossed fixed-length rods and three connecting springs (see fig. 1). This type of mechanism was studied in [START_REF] Arsenault | Développement et analyse de mécanismes de tenségrité[END_REF] in the particular case of symmetric geometric and loading conditions. The equilibrium configurations were solved for a set of geometric parameters and for one actuator input value. Here, these planar tensegrity mechanisms are analyzed in a more systematic way and in more detail, with the goal of understanding in depth the evolution of the number of stable and unstable solutions as a function of the geometric parameters, the loading conditions and the actuated joint inputs.

It turns out that the algebra involved in the stability analysis may prove very complicated while the planar tensegrity mechanisms at hand are rather simple, thus making a complete classification intractable when standard computer algebra tools are used. To overcome this difficulty, more sophisticated tools are used, namely, discriminant varieties and cylindrical algebraic decomposition combined with Groebner bases. These tools allow classification of the number of stable solutions as a function of various parameters.

Next section 2 introduces the architecture of the family of planar tensegrity mechanisms under study. In section 3, the direct kinetostatic problem is first analyzed for different geometric and loading configurations. A deeper analysis is then conducted in section 4 with the help of appropriate computer algebra tools, where a full classification of the number of stable solutions is established. Finally, section 5 is devoted to the study of solution continuity, namely, in which conditions the mechanism at hand can move while always remaining in the same stable solution. Last section adds some discussions and concludes this work.

MECHANISM ARCHITECTURE AND EQUILIBRIUM CONDITION

The mechanism architecture under study is shown in Fig. 2. It is made of two crossed rigid rods A 1 A 3 and A 2 A 4 of lengths L 1 and L 2 , respectively and three linear springs of stiffness k connecting A 1 A 4 , A 2 A 3 and A 3 A 4 . A reference frame is placed at the fixed point A 1 and the x-axis is along A 1 A 2 . Point A 2 can be translated along the x-axis with a prismatic actuator. This mechanism has three degrees of freedom. The first one, referred to as , is controlled by the actuator to raise or lower the mechanism. The other two variables,  1 and  2 , are generated by the compliant rotations of the two rods A 1 A 3 and A 2 A 4 about A 1 and A 2 , respectively.

A 4 A 3 A 2 A 1  l 2 l 3 l 1 L 1 L 2 F 4 F 3  1  2

Fig. 2: Planar tensegrity mechanism studied

Two vertical forces F 3 and F 4 are applied at nodes A 3 and A 4 , respectively. In this study, zero-free-length springs are considered. This is not a purely theoretical hypothesis since equivalent zero-free-lengths springs can be designed as shown for example in [START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-dof tensegrity mechanism[END_REF] and [START_REF] Arsenault | Développement et analyse de mécanismes de tenségrité[END_REF]. In the absence of friction and neglecting the compliance of the rods and of the prismatic joint, the potential energy U of this mechanism can be written as follows:

3 2 3 3 4 4 1 2 i i k U l F y F y      ( 1
)
where y 3 and y 4 are the ordinates of A 3 and A 4 , respectively. Expressing the spring lengths l i as a function of the other geometric parameters using the cosine laws, U takes the form below:

            2 2 2 1 1 1 2 1 2 1 2 2 2 3 1 1 4 2 2 3 4 cos 2( cos 4 ρcos ρsin ρsin 2 k U L L L L L L F L F L                  (2)
For the mechanism to be in equilibrium, the two derivatives of U with respect to  1 and  2 must vanish simultaneously:

      1 2 2 4 21 sin 2 sin cos 0 kL k F          (3)       2 1 1 3 21 sin 2 sin cos 0 kL k F          (4)
The solutions to the direct kinetostatic problem of the mechanism for a given input , are obtained by solving Eqs (3) and ( 4) for  1 and  2 , respectively.

SOLUTIONS TO THE DIRECT KINETOSTATIC PROBLEM

GENERAL CASE

The first step is to write equations ( 3) and ( 4) in polynomial form by using the tan-half-angle substitutions

t 1 =tan( 1 /2) and t 2 =tan( 2 /2) : 2 2 2 2 2 2 2 3 1 2 2 1 2 2 2 1 2 1 3 1 3 2 2 1 2 2 1 3 2 2 4 2 4 0 2 F t t kL t t kL t t k t t F t F t kL t kL t k t F            (5) 2 2 2 2 2 2 2 4 1 2 1 1 2 1 2 1 1 2 4 1 4 2 1 1 1 2 2 4 2 2 4 2 4 0 2 F t t kL t t kL t t k t t F t F t kL t kL t k t F            (6)
Upon elimination of one of the variables (e.g. t 2 ) in the above two equations, a polynomial of degree six in t 1 is obtained after clearing the factor (1+t 1 2). For each root, Eqs ( 5) and ( 6) can be combined to eliminate the terms of degree 2 and t 2 is then solved with a linear equation. Thus, the mechanism may have up to 6 solutions, of which some of them might be associated with instable equilibriums. Stable solutions are those solutions for which the 22 Hessian matrix H is definite positive, namely, its leading principal minors H(1,1) and det(H) are greater than zero. However, the general symbolic expressions of these minors are very large and stable solutions cannot be sorted out systematically at this stage.

We now inspect particular conditions that decrease the degree of the solution polynomial or lead to interesting, symbolically tractable cases.

SYMMETRIC DESIGN AND EQUAL FORCES

Here, a symmetric design L 1 =L 2 under symmetric loading F 3 =F 4 is examined. This situation was already studied by Arsenault [START_REF] Arsenault | Développement et analyse de mécanismes de tenségrité[END_REF], but under the assumption that the mechanism should always remain with its base and upper sides parallel, namely, all solutions should satisfy y 3 =y 4 , or equivalently,  1 = 2 . Accordingly, the direct kinetostatic problem was solved with only one equation, namely, the derivative of U w.r.t. y=y 3 =y 4 . A 4 th -degree polynomial equation was then obtained in y= y 3 =y 4 and solved for one particular set of parameters. In fact, it is not proven a priori that the solutions are always of the form  1 = 2 . For the sake of exhaustiveness, the direct kinetostatic problem is solved here with  1 not necessarily equal to  2 . Accordingly, the two equations ( 5) and ( 6) are used to solve the equilibrium conditions, in which F 3 and L 2 are replaced with F 4 and L 2 , respectively. To get simpler expressions, the second equation is subtracted from the first one. The tan-half substitution is then applied in the resulting equation and the following new system is obtained:

1 2 4 1 2 1 2 ( )( ( ) 2 (1 ) 0 t t F t t k t t       (7) 2 2 2 2 2 2 2 4 1 2 1 1 2 1 2 1 2 4 2 1 1 1 2 2 2 ( 4 ( (1 2 ( 0 ) ) ) ( ) ) F t t L k t t t t k t t F t t L k t t t            (8)
The first factor (t 1  t 2 ) in [START_REF] Aldrich | Time-energy optimal control of hyper-actuated mechanical systems with geometric path constraints[END_REF] confirms that solutions of the type  1 = 2 do exist but the second factor indicates that solutions with distinct angles may also appear.

The solutions obtained from the first factor of ( 7) are calculated by substituting t 2 =t 1 in eq [START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-dof tensegrity mechanism[END_REF], which yields up to 4 solutions.

The other solutions are obtained upon eliminating t 1 from the second factor of (7) and from [START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-dof tensegrity mechanism[END_REF], which leads to a polynomial of degree 2 in t 2 :

2 2 2 3 3 3 2 2 2 4 2 2 2 4 2 4 2 4 4 ( ) 16 ( 1) 4 4 ( ) 0 F k t L k t F t F k t F k L               (9)
Since t 1 can then been solved linearly using the second factor of ( 7), there are up to two solutions of the type  1 ≠ 2 . Moreover, the two solutions are of the form (t 1 , t 2 ) and (t 2 , t 1 ) since the same polynomial as ( 9) could have been obtained in t 1 by eliminating t 2 instead of t 1 .

We now study the stability of the solutions. The leading principal minors H(1,1) and det(H) of the Hessian matrix must be calculated and their sign must be positive for an equilibrium solution to be stable. The symbolic expression of det(H) for the solutions  1 ≠  2 is surprisingly simple and the sign analysis can be handled easily.

Indeed, solving the second factor of (7) for t 2 and replacing the solution in det(H) leads to the following factored expression:

2 2 2 2 2 2 2 2 2 4 2 4 2 4 24 4( 1) (4 ) ( ) ( 2 ) det( ) t k F k t F t k k t F H            (10)
It can be observed that det(H) is always negative. Thus, the two equilibrium solutions of the type 12   are always unstable. On the other hand, the sign study of det(H) and H(1,1) for the solutions of the type  1 = 2 is not straightforward and requires tedious analyses. Stability analysis of those solutions is left to section 4 with the use of more appropriate tools.

NO EXTERNAL LOADING (F 3 =F 4 =0)

When F 3 =F 4 =0, Eqs (3) and ( 4) become:

    1 2 1 2 sin 2ρsin 0 L       (11)     2 2 1 1 sin 2ρsin 0 L       (12)
It is apparent that  i = 0 or π, i=1, 2 are solutions to the above system, which define four singular configurations.

In such configurations, the mechanism is fully flat and cannot resist any vertical forces. There are two more solutions of the form:

1 1 2 2 1 1 2 2 (( arctan( / ), arctan( / )), arctan( / ), arctan( / )) Q L Q L Q L Q L               (13)
where

  2 2 2 2 1 2 1 2 ( ) 4 )( ) 4 ) LL Q LL      
. On the other hand, when the coordinates of A 3 and A 4 are calculated with the above solutions, one obtains y 3 = y 4 and x 3  x 4 = . This means that the mechanism remains always in a parallelogram configuration, even when L 1 ≠L 2 .

It can be shown easily by reporting  1,2 = 0 or π into det(H) and H(1,1), that three solutions of the four flat ones are always unstable.

Solutions classification

We know from previous section that the solutions are obtained upon solving a polynomial of degree six. In some particular cases, we were able to identify stable solutions. The goal of this section is to classify, in a systematic way, the number of stable equilibrium solutions as a function of the geometric and physical parameters of the planar tensegrity mechanisms. The algebraic problem at hand relies on solving a polynomial parametric system of the form:

  11 ( ) 0, , ( ) 0, ( ) 0, , ) 0 , (

ml n p p q E q      K ¡ K vv v vv (14)
Such systems can be solved in several ways. Discriminant varieties [START_REF] Lazard | Solving parametric polynomial systems[END_REF][START_REF] Corvez | Using computer algebra tools to classify serial manipulators[END_REF][START_REF] Manubens | Cusp Points in the Parameter Space of Degenerate 3-RPR Planar Parallel Manipulators[END_REF] and cylindrical algebraic decomposition [START_REF] Collins | Quantifier Elimination for Real Closed Fields by Cylindrical Algebraic Decomposition[END_REF][START_REF] Corvez | Using computer algebra tools to classify serial manipulators[END_REF][START_REF] Manubens | Cusp Points in the Parameter Space of Degenerate 3-RPR Planar Parallel Manipulators[END_REF] are used here. They provide a formal decomposition of the parameter space through an algebraic variety that is known exactly. These tools have already been applied successfully for the classification of the singularities of serial and parallel manipulators [START_REF] Corvez | Using computer algebra tools to classify serial manipulators[END_REF][START_REF] Manubens | Cusp Points in the Parameter Space of Degenerate 3-RPR Planar Parallel Manipulators[END_REF], a class of problems that can be formulated as in [START_REF] Collins | Quantifier Elimination for Real Closed Fields by Cylindrical Algebraic Decomposition[END_REF]. Shortly speaking, discriminant varieties generate a set of separating hyper-surfaces in the parameters space of the parametric system at hand, such that the number of solutions in each resulting connected component or cell is known and constant. The discriminant varieties can be computed with known tools like Groebner bases using the Maple sub-package RootFinding [parametric]. Once the discriminant varieties are

obtained, an open cylindrical algebraic decomposition is computed to provide a full description of all the cells: the number of solutions in each cell is determined by solving the polynomial system for one arbitrary point in each cell. Finally, adjacent cells with the same number of solutions are merged.

The equilibrium solutions depend on three geometric parameters (the rod lengths L 1 and L 2 and the input variable) and two physical parameters (the spring stiffness k and the forces F 3 and F 4 ). Without loss of generality, however, the lengths parameters can be normalized with L 1 , while F 3 and F 4 can be normalized with k.

Finally the system at hand depends on four independent parameters only. In what follows, L 1 and k are fixed to 1 and 100, respectively, without loss of generality.

NO EXTERNAL LOADING

In this section, the unloaded situation F 3 =F 4 =0 is studied. We were able to show in the preceding section that three of the four flat solutions were unstable but no general information could be obtained about the two non-flat solutions. Since F 3 =F 4 =0, the parameter space is a plane defined by L 2 and . Calculating the discriminant varieties and the cylindrical algebraic decomposition for this case shows that there exists one region in the parameter plane where the mechanisms have two stable solutions. Outside this region, only one stable solution exists. Moreover, as already shown in section 3.1, the mechanisms take the shape of a parallelogram even when L 1 ≠L 2 .

The 2-solution region is of constant width equal to 1 when L 2 >1, while it decreases linearly when L 2 <1. In view of the optimal design, it is thus interesting to note that the largest operation range  = 1 is obtained for L 2 >1.

L 2  2 1.5 1 0.5 0 0 0.5 1 1.5 2 A 4 A 4 A 3 A 1 A 2 A 3

FULLY SYMMETRIC CASE

The mechanisms have now a symmetric design (L 1 =L 2 ) and symmetric loading conditions (F 3 =F 4 ). It was shown in the preceding section that the two solutions  1 ≠  2 were always unstable but no stability information could be brought for the four solutions of the type  1 =  2 . The parameter space is now the plane (, F 4 ), since L 1 =L 2 =1. There exists a region with two stable solutions and two regions with one stable solution, filled in red and blue, respectively. Since we know that the solutions corresponding to

GENERAL CASE

The parameter space is now 4-dimensional since it is defined by (, L 2 , F 3 , F 4 ). One parameter is assigned in order to have a parameter space of dimension 3. First, the design parameter L 2 is fixed. Figure 5a (resp. 4b, 4c)

shows the obtained partition of the parameter space (, F 3 , F 4 ) when L 2 =3/2 (resp. 1, 1/2) for 0<£2 and 0£ F 3 £10, 0£ F 4 £10. The boundaries are three surfaces of degree 12 in  2 and in F 3 , F 4 , defined by equations containing not less than 436 terms. Note that the regions have been also plotted on separate figures to show better the separating surfaces. These plots show that the regions with two stable solutions have a regular shape and are of reasonable size under quite small force magnitudes. It is apparent that when L 2 is increased, the size of those regions is also increased.

It is interesting to see how the size of the 2-solution regions varies when much higher forces are applied. Figure 6 shows how the 2-solution region shrinks and even disappears when the force magnitudes reach certain limits.

Since the discriminant varieties give an algebraic description of the separating surfaces, it is possible to determine exactly the singular points of these surfaces and in particular the coordinates of transition point C. Two parameters are now assigned in order to have a parameter space of dimension 2. Accordingly, the discriminant varieties and cylindrical algebraic decomposition are computed for F 3 =F 4 = 10. Figure 7 (left)

shows the obtained partition of the parameter plane (, L 2 ) for 0<£2 and 0£ L 2 £2. It looks similar to the unloaded case but the boundaries here are three curves of degree 12 in  and in L 2 and their equations are too large to be displayed in this paper. There are two stable solutions in the red region and only one in the blue region. Figure 7 (right) shows two stable solutions for L 2 =3/2 and =7/10.

C

known as a catastrophe, a phenomenon often observed in mechanisms with springs [START_REF] Arsenault | Développement et analyse de mécanismes de tenségrité[END_REF][START_REF] Hines | Catastrophe analysis of the planar two-spring mechanism[END_REF]. Figure 8b shows that when L 2 is slightly increased, the first critical point C 1 is transformed into a fold point. ttand of 2 t in Eq. ( 6)). Globally one comes up with a system yielding six solutions, of which one or two are stable like. The gravity has been omitted in this work. Again, taking into account the rod loads does not modify the nature of the problem. It just adds additional terms in the equations that result in slightly modified equilibrium solutions. The three springs were assumed similar but taking different stiffness's

does not change the nature of equations and of the general results as well. In fact, playing with additional force components or distinct stiffness values increases the dimension of the parameter spaces. This is not problematic

C 1 C 2 C 2 C 1
for the calculation of the discriminant varieties but the computation time of the cylindrical algebraic decomposition increases and graphical representations are more challenging. Using springs with non-zero free lengths make the equations much more complex and might render symbolic calculations intractable, as shown in [START_REF] Wenger | Kinetostatic Analysis and Solution Classification of a Planar Tensegrity Mechanism[END_REF] and already observed in former works [START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-dof tensegrity mechanism[END_REF][START_REF] Crane | Kinematic analysis of a planar tensegrity mechanism with pres-stressed springs[END_REF]. Practically, however, equivalent zero-free lengths springs can be designed without much effort. Finally, interferences between the two crossed rods can be easily avoided by assembling them in two different layers, as shown in [START_REF] Arsenault | Développement et analyse de mécanismes de tenségrité[END_REF] for example.

CONCLUSIONS

The goal of this paper was to investigate in depth the direct kinetostatic solutions of a family of planar tensegrity mechanisms composed of a telescopic base, two crossed rods and three springs. It was shown in this paper that the problem can be treated using suitable computer algebra tools. We have shown that a univariate polynomial of degree 6 must be solved in the general case, resulting in one to two stable solutions. In the unloaded case, there are always two stable symmetric solutions for a range of the input prismatic joint which is of constant width and whose limits vary with the rod lengths. Moreover, the mechanisms remain always in a parallelogram configuration even when their two rod lengths are different. The mechanisms can also reach one flat stable solution, which is singular. Such a stable flat solution might be of interest to store the mechanism when it is not used. When the two external forces and the two rod lengths are equal, there are still 6 solutions, including 4 unstable, non-symmetric solutions.

We have used discriminant varieties and cylindrical algebraic decomposition to study the evolution of the number of solutions as function of various parameters. Those tools have already proved efficient to treat several difficult kinematic issues arising in serial and parallel manipulators. Discriminant varieties and cylindrical algebraic decomposition generate a set of separating curves or surfaces in the parameter space, which are fully defined with algebraic equations.

The next step of this work is to consider and design an assembly of proper planar tensegrity modules stacked in series and suitable actuation to get a first planar model of a bird neck. There are many different choices in the design of the elementary tensegrity module by playing with the choice of the actuated link(s) and the number of springs [START_REF] Arsenault | Développement et analyse de mécanismes de tenségrité[END_REF]. It is possible to actuate a second rod, namely, one of the crossed ones. The classification analysis conducted in this paper can be applied directly to those mechanisms. It is also possible to select a mechanism with only two springs along the two sides and all other segments made with fixed-length rods. Such a mechanism has only one degree of freedom in total (a rotation about the intersection point of the two crossing rods) and can be driven with cables [START_REF] Bakker | Design of an environmentally interactive continuum manipulator[END_REF]. Figure 9 illustrates a possible assembly of three such mechanisms with two driving cables thread through the spring attachment points as described in [START_REF] Bakker | Design of an environmentally interactive continuum manipulator[END_REF]. Depending on the number of parameters to be controlled (position only or position and orientation of the terminal link, stiffness, configuration of the whole serial assembly), the arrangement and actuation of the cables have to be properly chosen.

Fig 1 :

 1 Fig 1: A serial assembly of elementary tensegrity mechanisms mimicking an animal neck (illustration from [7])

Figure 3 (

 3 left) shows a representation of the parameter plane obtained for L 2 and  in[0,[START_REF] Vogel | Cats Paws and Catapults: Mechanical Worlds of Nature and People[END_REF]. The 2solution region is shown in red. The three separating lines that bound the regions are defined by 2L 2 1=0, 2L 2 +1=0 and 2+L 2 1=0, respectively. Here it can be verified with simple geometric arguments that these boundaries define the limit conditions for the mechanisms to become fully flat. In the three 1-solution regions, the mechanisms have one stable fully flat solution. In those regions including their boundaries, the mechanisms are thus in singular configurations. They can move along a horizontal direction but they cannot resist any vertical force. Such solutions could be used to store the mechanisms when there are not in use. In the 2-solution regions, there are two stable non-flat solutions, one being the mirrored image of the other as shown in Fig.3(right).

Fig. 3 :

 3 Fig. 3: Unloading case: parameter plane (left) and the two mirrored stable solutions for =1, L 2 =3/2 (right)

  The computed discriminant varieties and cylindrical algebraic decomposition provides the parameter plane shown in fig.4a(left) for 0<£2, 10£ F 4 £0 (pushing forces) and in fig.4bfor 0<£2, 0£ F 4 £10 (pulling forces).

4 A 4 A 3 A 2 A 3 A 1 Fig. 4a : 1 Fig. 4b :

 4432314a14b Fig. 4a: Symmetric case: parameter plane under pushing forces for 0<£2, 10£ F 4 £0 (left) and two stable solutions for =3/4, F 4 = 10 (right)

  Fig. 4c: Stable solution in each of the two 1-solution regions under a pushing force F 4 =10 (left: =2/10, right: =3/2)

Fig. 5a :Fig. 5c :

 5a5c Fig. 5a: Partition of the parameter space when L 2 =3/2 for 0<£2 and 0£ F 3 £10, 0£ F 4 £10

Fig. 6 :

 6 Fig. 6: The 2-solution region for large ranges of force magnitudes. Point C defines the limit force magnitudes for two stable solutions to exist.

Fig. 8a : 6 DISCUSSION

 8a6 Fig. 8a: Solution continuity with L 1 =L 2 =1, F 3 =F 4 =5 (stable solutions in blue)

Fig. 9 :

 9 Fig. 9: An assembly of 2-spring tensegrity mechanisms.

are unstable. Figure4a, right (resp. 4b, right) illustrates two stable solutions for  =

3/4 and F

Acknowledgments This work was funded by the French ANR project "AVINECK: an arm for the bird", ANR-16-CE33-0025-02, 2017-2020.

 [START_REF] Vogel | Cats Paws and Catapults: Mechanical Worlds of Nature and People[END_REF]= 10 (resp.  = 3/4 and F 4 = 10). Contrary to the unloaded case, those two solutions are not symmetric.

The boundaries between the 2-solution and the 1-solution regions are two curves of degree 6 defined by:

In the two 1-solution regions, it can be shown that the mechanisms operate always with y 3 and y 4 negative under pushing forces. Figure 4c shows stable solutions selected in each of the 1-solution regions. Assuming that a mechanism is controlled to start in a configuration with y 3 and y 4 positive, its operation range  for a given pushing F 4 is thus determined by the 2-solution region only. The operation range decreases when the magnitude of external force increases. It can be verified that for F 4 =0,  reaches its maximal value, equal to 1, in accordance with the result in section 3.1. In the presence of pulling forces (F 4 >0), y 3 and y 4 are positive in the two one-solution regions. Figure 4d shows stable solutions chosen in each of the two 1-solution regions. Like in the preceding case, the operation range is determined by the 2-solution region if the mechanism starts with y 3 >0 and y 4 >0. The operation range reaches its maximal width for L 2 =1, which is the fully symmetric case:

it can be verified that this range is exactly the same as the one calculated in the above section for F 3 =F 4 =10.

Starting from L 2 =1, the operation range then decreases slowly when L 2 increases but the decrease is much more significant when L 2 decreases.

Solution continuity

Section 4 above shows that one or two stable solutions may exist for a given value of the actuator parameter . If a mechanism is actuated to follow a prescribed path, it is necessary to verify that the starting solution remains available during the whole trajectory to avoid any undesirable jump to another solution during motion. In this section, the solution continuity is analyzed. The case with no external forces is clear since the 1-solution regions are associated with only singular configuration that cannot be used for practical motions and the mechanisms must be operated by starting and ending any trajectory in the 2-solution region only. In the fully symmetric case, the 2-solution region is surrounded with two 1-solution regions. Figure 8a shows the evolution of the solutions for a mechanism defined by L


X Giraud

M Budinger

X Roboam

H Piquet

M Sartor

J Faucher

Optimal design of the Integrated Modular Power Electronics Cabinet

Keywords: more electrical aircraft, power electronics, optimal design, multi-physical sizing, combinatorial optimization

This paper deals with a new concept of electrical power distribution system called IMPEC for Integrated Modular Power Electronics Cabinet. The paper defines methods aiming at carrying out an optimal design of IMPEC, the main variables being on one hand the number and size of power electronics module. On the other hand, reconfigurations between these modules and electrical loads are also optimized.The formalization of the problem highlights that designers must deal simultaneously with a combinatorial explosion and a multi-physical system sizing. The main objective of the study is to propose a methodological framework for solving this original optimal design problem. A heuristic-based algorithm is developed to solve this combinatorial optimization problem. A particular attention is paid to develop a weight estimation procedure using generic sizing models. Finally a mapping is performed to identify the best solutions and to highlight the technological components having the most significant sensitivity on the complete system weight.

I. Introduction and context I.1 More Electrical Aircraft

On traditional aircraft such as A320 or A330, systems are powered by 3 different energy vectors: pneumatic, hydraulic and electric (Figure 1 and Figure 2). These 3 power vectors are extracted from primary power sources such as the engines or the auxiliary power unit (APU) that are today all supplied by fuel (kerosene) [START_REF] Moir | Aircraft Systems -Mechanical, electrical, and avionics subsystems integration[END_REF].

In terms of power levels, the pneumatic vector is the most demanding one. During some flight phases, it represents around 80% of the power taken off the engines (Figure 1). This energy vector supplies the engine starting system, the wing anti-ice protection system (WIPS) as well as the environment control system (ECS) by extracting pressurized air taken from the engines (bleed air).

Systems requiring high force at low speed such as flight control surface actuation, landing gear actuation and aircraft breaking are traditionally powered by the hydraulic vector.

Eventually, the electrical power system (EPS) provides power to numerous and various systems (also called electrical loads). On-board an Airbus A330, around 700 electrical loads are embedded. Some of them are essential such as: the fuel pumps, the electro-hydraulic pumps, the fans, the processing blades ; others are installed in order to provide comfort to the passengers: the galleys and the flight-entertainment systems. For the last decade, the use of electricity is rapidly increasing in commercial aircraft. This trend illustrated in the paradigm of the More Electrical Aircraft (MEA) consists in replacing hydraulic and/or pneumatic powered systems by electrical ones [START_REF] Rosero | Moving towards a more electrical aircraft[END_REF] [4] [START_REF] Roboam | More electricity in the air: Towards optimized electrical networks embedded in "more electrical aircraft[END_REF]. As a consequence, the realization of the "allelectrical" aircraft aims at completing simultaneously 2 types of configurations: Bleedless [START_REF] Sinnett | 787 No-Bleed Systems: Saving Fuel and Enhancing Operational efficiencies[END_REF] and Hydraulicless [START_REF] Todeschi | Airbus -EMA's for flight controls actuation systems -2012 status and perspectives, perspectives[END_REF] [8] (Figure 3). By embedding only 2 hydraulic circuits instead of 3 traditionally, the Airbus A380 and A350 are seen as an intermediate configuration towards hydraulicless aircraft [START_REF] Van Den Bossche | A380 primary flight control actuation system[END_REF]. On the other side, the Boeing 787 is a complete bleedless aircraft by supplying electrically: the ECS, the WIPS and the engine starting [START_REF] Sinnett | 787 No-Bleed Systems: Saving Fuel and Enhancing Operational efficiencies[END_REF].

Today, the MEA is seen as a major axis of improvement for the aviation industry to achieve increasingly ambitious objectives: decrease of weight, rationalization of costs, reducing environmental impact... In this frame, several research projects are today launched in order to progress on architecture and technology axis [START_REF]European research project[END_REF] [START_REF]European research project[END_REF].

I.2 High demanding loads requiring power electronics in bleedless configuration

In terms of power demand, the bleedless aircraft configuration is the most challenging one. It is often said that the bleedless effect leads to multiply by 4 the electrical power generation [START_REF] Roboam | New trends and challenges of electrical networks embedded in "more electrical aircraft[END_REF]. For instance, the Boeing 787 comprises 4 generators of 250 kVA whereas a bigger aircraft such as the Airbus A380 "only" embeds 4 generators of 150 kVA.

The ECS is the main contributor of the power step change. The ECS electrification requires installing electrical air compressors consisting of permanent magnet synchronous motors driven by power electronics modules (i.e inverter). The compressors can consume up to 100 kW during some operational cases. As the consequence, the bleedless effect will lead to a significant increase of the power amount transmitted by the power electronics modules (PEM).

I.3 Classical use of power electronics in aeronautics

On today aircraft, an electrical machine is supplied by a dedicated PEM often located close to the actuator. Depending on the voltage levels provided by the EPS, two structures appear (Figure 4) :  When the load is allocated to an AC busbar (Figure 4-A), a rectifying stage is required to create a local HVDC bus (+/-270 VDC).

 When HVDC busbars are provided by the EPS (Figure 4-B), the rectifier stage is removed and the PEM is directly connected to a HVDC busbar.

These two structures where each load has its dedicated PEM have 3 major drawbacks:

 Reliability aspect. The PEM loss leads inevitably to the load loss. This aspects can have an impact on the availability/reliability of the aircraft function carried out by the load.

 Cost aspect. Each PEM is sized by its load. On the Figure 4 example, two different PEM sizes are required: 40 kW and 75 kW. This leads to an increase of the amount of equipment references also called "part numbers". This list of equipment shall be as low as possible in order to increase standardization and reduce costs. Of course, it would be possible to have only one PEM reference of 75 kW but this would lead to an oversizing and therefore a weight increase.

 Utilization rate. A given set of loads or aircraft functions are only supplied during a reduced period of the mission. For example, these 3 systems are not active during the same operational case: the engine starting, the landing gear actuation and the engine thrust reversers…Electrifying each system with one PEM each could lead to a low PEM utilization rate.

Load n°1

Load n°2 In the general context presented here before, the present paper is organized into 6 parts:

 the section n°2 describes the IMPEC concept being the system under study in this paper. The IMPEC organic structure is presented as well as the advantages of this new power distribution concept.

 the section n°3 is dedicated to the presentation of the optimization problem associated to the IMPEC optimal design. The two main sources of problem complexity are identified and illustrated with examples.

 the section n°4 aims at providing an overview of the sizing procedure deriving the overall system weight and features of IMPEC implementations.

 the section n°5 presents the optimization algorithm developed to solve the optimal design problem.

 the section n°6 is a post-optimal analysis aiming to identify the most robust solution forms and the component technology sensitivities on the overall system weight.

II. The IMPEC concept

In aeronautics, a new concept enabling to supply several loads requiring an inverter is under study. IMPEC concept stands for Integrated Modular Power Electronics Cabinet [START_REF] Prisse | New power centre and power electronics sharing in aircraft[END_REF]. It is made of an electrical power center hosting a set of PEMs that have the same features (modularity/standardization) and that can supply several loads included in the power center perimeter (mutualisation).

II.1 State of the art on modularity and mutualisation of electronics in aircraft

Modularity and mutualisation of electronics embedded in Airbus aircraft are only implemented for processing electronics. On conventional aircraft (e.g A320), each system has one or several processing blades covering its own calculation need. With the A380 launched in the mid-2000, the IMA concept (Integrated Modular Avionics) has been implemented in order to minimize the number of equipment references (part numbers) [START_REF] Ramsey | Integrated Modular Avionics : Less is More[END_REF] [START_REF] Prisaznuk | Integrated modular avionics[END_REF] . The main benefits expected from this concept are cost and reliability. The IMA consists of a set of standardized processing blades and a communication network ensuring the data transmission. An important mutualisation is achieved since one processing blade can host applications coming from different systems.

Mutualisation of power electronics is a relatively new concept in aeronautics and more generally in the industrial fields linked to electrical engineering. Two types of mutualisation can be highlighted:

1) A first one consists in using a same PEM to supply different loads at different instants according to the operational needs. This is today applied for the Bleedless aircraft Boeing 787. One PEM (called MCU for Motor Control Unit in Figure 6) ensures the engine starting [START_REF] Roboam | New trends and challenges of electrical networks embedded in "more electrical aircraft[END_REF]. Once the engines are running, the same PEM supplies the ECS electrical compressor for the rest of the mission. Thus, the same PEM supplies two different systems: "Engine starting" and "ECS". 2) A second type of mutualisation consists in associating several PEMs together in order to supply a same load. This concept, introduced in [START_REF] Prisse | New power centre and power electronics sharing in aircraft[END_REF], enables to obtain different levels of electrical production capacities. For instance, with 4 PEMs of 30 kW each, it is possible to form four different power levels: 30, 60, 90 and 120 kW. An additional flexibility is therefore obtained in order to stick as much as possible to the operational needs of the electrical consumers.

II.2 IMPEC features

II.2.1 Architecture and operation modes

IMPEC implements the 2 types of mutualisation concepts presented here before. The IMPEC structure is made of 2 main parts: a set of modular PEMs (all being identical) providing power through HVDC bus and a contactor matrix ensuring the connexions between the PEM and the loads. In comparison with the IMA, the processing blades are the PEMs and the contactor matrix plays the role of the communication network. The main functioning principles are the following ones:  One PEM can only supply one load at the same time;  One PEM can supply different loads during the flight mission: it can supply one load during the "take-off" phase and another load during the "cruise" phase;  One load can be supplied by one or several PEMs at the same time (PEM paralleling);

The Figure 7 shows a principle example of an IMPEC implementation made of 4 PEMs of 30 kW each. An example of the contactor matrix configuration is shown for a case during which all the PEMs are available (right side of Figure 7).

II.2.2 IMPEC components and architecture

The Figure 8 shows the organic architecture of an IMPEC implementation made of only 2 PEMs. According to the power paths formed by the contactors, the 2 PEMs supply 1 load each. In more details, the physical structure of the IMPEC is the following one:  PEMs made of 3 devices: an inverter part based on IGBT modules (component n°1 of Figure 8), a cold plate (2) and a DC link capacitor (3).

 Single phase inductances (4) at the output of each PEM leg. These inductances are not part of the standardized PEM. They are mandatory for the PEM mutualisation (i.e PEM paralleling). Thus, a PEM which is never put in parallel with another PEM will not require having output inductances.

 A contactor matrix [START_REF] Roboam | More electricity in the air: Towards optimized electrical networks embedded in "more electrical aircraft[END_REF] ensuring necessary connexions between modules and loads. A contactor ensures the connexion between one module and one load. As shown on the Figure 8, the closure of the "contactor 1-2" connects the PEM 1 to the load 2.

 A power center chassis ( 6) integrating mechanically all the components listed here before.

 A heat exchanger ( 7) evacuating the heat extracted by the cold plates.

As it will be explained in the section III, the definition of the commutation components (the inductances and the contactor matrix) is strongly impacted by the IMPEC reconfiguration represented by the different connections : module/load. Depending on the number of contactors installed in the contactor matrix, it is possible to connect any load to any PEM. With IMPEC, the loss of one PEM does not lead to a load loss whereas it is the case with a classical structure (each load has its dedicated PEM).

On the Figure 9, this feature is illustrated by the examples on which two different PEM loss configuration are highlighted. As shown on the Figure 7 (right), in the nominal configuration, the load n°1 is supplied by the PEM1 & 2 and the load n°2 by the PEM4. When there is a failure of the PEM2, the PEM3 takes over the supply of the load n°1. When this is only the PEM4 which is inoperative, the load n°2 is still supplied thanks to the PEM3. The contactor matrix offers the possibility to have better matching between the PEMs use during the flight mission and the power demand fluctuations of the loads. This capability is illustrated by the example showing two cases (Figure 10). In the case n°1, the load n°1 and n°2 are consuming respectively 75 kW and 30 kW: 3 PEMs are then required for the load n°1 and only 1 PEM for the load n°2. When the load n°2 consumes its maximum power (40 kW) in the case n°2, both PEM3 and PEM4 are required. The PEM3 has been "freed" by the load n°1 because its consumption has decreased in comparison with the case n°1.

O F F

II.3 Studied use case and requirements

All the design methods developed in this paper are applied and tested on an industrial use case with a bleedless aircraft configuration. As shown on the Figure 11, IMPEC shall supply 6 different electrical loads: 3 electrical air compressors of the environmental control system (loads: ECS1, ECS2, ECS3), 1 electrical pump for the fuel tank inerting system (load: FTIS), 1 load modelling the cabin systems supplied by the IMPEC (load: CS) and 1 load representing the electrical motor ensuring the engine starting (load : ENG START). As depicted by the Figure 11, the number of PEMs is not assumed, being a degree of freedom for system design. The value of this variable will be part of the optimization process in order to reach an optimal solution in terms of IMPEC system weight. The optimal design problem comprises 2 types of decision variables:

 the number of PEMs (modules) « » ;

 the reconfiguration solution defining all the connections load/module for all the operational cases.

A reconfiguration solution is modeled by the 3-dimension matrix « » storing the binary variables « , ». These variables are defined as follows:

∀ ∈ {1, . , } , =
1, load " " is connected to the PEM " " for the case "c" 0, otherwise Eq 1

∀ ∈ {1, . , ℒ}, ∀ ∈ {1, . , },
Where « ℒ » is the number of loads to be supplied by IMPEC, and « » is the number of loading cases to be considered for the design. The loading cases represent all the operational cases that shall to be taken into account in the design [START_REF] Giraud | Load allocation problem for optimal design of aircraft electrical power[END_REF]. A loading case is defined by 4 dimensions :  a flight phase (e.g take-off);  a mode of the IMPEC electrical system (e.g loss of 2 PEMs);  a mode of the other systems (e.g loss of one ECS electrical compressor);  an external condition (e.g hot day).

As one dimension can have many values (for example, 21 flight phases are considered), several thousands of loading cases « » are considered for the industrial use case. The use of the 2 design variables introduced here before is illustrated by an example displayed on the Figure 12. The IMPEC structure is made of « = 5 ». The reconfiguration solution « » is graphically represented by a parallelepiped. The configuration of the loading case n°1 « = 1 » is highlighted.

Eventually, it must be noted that according to this configuration of the case n°1, at least 5 contactors are required between : L1-PEM1 ; L2-PEM3 ; L3-PEM2 ; L4-PEM4 ; L4-PEM5. In addition, as PEM4 and PEM5 are in parallel to supply the FTIS load, inductances are required at the outputs of both PEMs. Where « » is the overall IMPEC system weight, « » and « » are respectively the number of PEMs and their weight; « » and « » are respectively the number of inductances and their weight ; « » and « » are respectively the number of contactors and their weight ; « » is the chassis weight and « » is the heat exchanger weight.

In order to reduce the number of "part numbers" and to push the standardization as far as possible, all components of the same type have the same reference (i.e the same features: dimensions, weight). It means that there is one reference for : the IGBT modules, the capacitors, the cold plates, the inductances and the contactors.

With a view to differentiate the organic constitution of a solution and its management/control, an IMPEC solution « » is defined by the couple « = ( , ) » where :  « » defines the physical constitution of the solution : the number of PEMs « », the number of inductances « », the number of contactors « », etc.  « » represents the reconfiguration solution defined in the part III.1.

III.1.3 Design constraints

The first design constraint imposes that a consuming load shall be supplied and therefore connected to one available PEM:

∀ { , } ∈ {1, . . , ℒ} × {1, . . , }, > 0 ℎ ∑ , ≥ 1 Eq 3
Where « » is the power demand of the load « » in the loading case « ». The dual constraint specifies that a PEM shall not be connected to a non-consuming load.

A second constraint specifies that one PEM supplies only one load in a given loading case:

∀ { , } ∈ {1, . . , } × {1, . . , }, ∑ , ℒ ≤ 1
Eq 4

The last main constraint consists in verifying that the power demand of one load does not exceed the power capacity of the PEM(s) connected to it. In terms of current, this is formulated by the following equation: As highlighted in the Figure 13, the IMPEC architecture mixes components that shall be designed by considering different physical fields. For example, IGBT modules shall be defined by considering simultaneously thermal and electrical aspects. Moreover interactions exist between the different component size. For instance, the design decisions on the cold plate (e.g the fluid flow rate) will have an impact on the IGBT cooling as well as on the heat exchanger. This problem complexity requires applying a systemic approach by considering all devices and corresponding physical fields as a whole. According to the perimeter and objectives of the study, the integrated models shall have the "right" level of details, being "just enough accurate" due the complexity of the optimal design process. This methodology developed to cope with problem complexity is presented in the part IV of this paper. The « , » values (Eq 1) shall be decided in order to define a reconfiguration solution « ». The number of possibilities is very high. Let us take a reduced-size example having the following features :

« = 3 », « = 4 », « ℒ = 2 ».
The Figure 14 shows all the matrix contactor possibilities for each loading case (the number of parallelized PEMs is limited to 2). There are 12 matrix contactor choices for the case n°1; 2 for the cases n°2 and 3; 6 choices for the case n°4. In a graphical way, the reconfiguration problematic consists in selecting one square for each case. The combination of the squares form a reconfiguration solution. On the Figure 14 Any combination of configurations leads to a reconfiguration solution. Thus the size of the design space is set by 2 features : the number of cases « » and the number of potential contactor configuration per case « ( ) ». The total number of different reconfiguration solutions « », assessing the combinatorial space, is given by:

= ∏ ( ) Eq 6
Applying this formula, the total solution number of the Figure 14 example is: 12 × 2 × 2 × 6 = 288. For a problem having 5 consuming loads and 6 available PEMs, the number of potential matrix configurations for one case is : « ( ) = 2520 ». For a problem made of 50 loading cases (still not in line with several thousand loading cases mentioned in the section III.1.1), the total number of reconfiguration solutions is given by: = 2520 ≈ 10

Eq. 7

This simple example shows the highly combinatorial feature of the matrix reconfiguration problem. A combinatorial optimization approach is required in order to handle this complexity: an algorithm dedicated for the IMPEC application is presented in the section V.

Note: it is interesting to notice that this combinatorial explosion is also encountered by the other aircraft shared resource system: the IMA [START_REF] Ramsey | Integrated Modular Avionics : Less is More[END_REF]. As a consequence the IMA optimal design also requires to build methodologies based on the use of combinatorial optimization algorithms [START_REF] Annighofer | Automated selection, sizing, and mapping of Integrated Modular Avionics Modules[END_REF].

L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L2 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 L1 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3 M1 M2 M3

IV. Evaluation function : multi-physical sizing process

The multi-physical sizing process aims to compute a relevant system weight « » according to the two design variables « » and « ». As displayed on the Figure 15, a modular design methodology is intended to be developed by cutting the sizing process into different component blocks. Each block aims to design one component type: IGBT module, cold plate, etc. Two important axis are worked in order to have a computation-efficient sizing process with regards to high combinatorial aspect of the IMPEC optimal design problem:

 Axis 1: the sizing process organization: what is the "less expensive" calculation sequence and the interactions (inputs/outputs) between the different blocks?

 Axis 2 : the sizing models: what are the methods and the hypothesis to be applied in order to compute the main features of components?

State of the art

There are several different technics aiming to help the designer to organize a process. The formalisms can be different but they always manipulate generic notions:

 A decomposition of the process by entities ;  A description of the interactions or interdependencies between these entities.

The causal diagram is a graphical representation in which each entity is modelled by a block and the interactions are identified by oriented arrows. This formalism was used to organize the sizing process of a flight control actuator [START_REF] Liscouet-Hanke | A Methodology for Systems Integration in Aircraft Conceptual Design -Estimation of Required Space[END_REF]. On a causal diagram, it may be difficult to rapidly visualize the couplings between the different entities. They can be automatically detected either by a dedicated tool [START_REF] Delinchant | A Component-Based Framework for the Composition of Simulation Software Modeling Electrical System[END_REF] or by formalizing the process using another type of diagram : the N-square Diagram (N2D) [START_REF]Techniques of Functional Analysis[END_REF].

As the causal diagram does, a N2D models a process by a set of entities (blocks) and interactions (arrows). But the N2D introduces a graphical convention in order to have a better visualization of the block interactions : inputs are vertical arrows and outputs are horizontal arrows (Figure 16). Thanks to this formalism, the coupling of the Figure 16 example can be rapidly identified : between block 2 and 3 through the parameter 3 and 4.

The causal diagram and N2D can be mathematically represented by a Design Structure Matrix (DSM) [20]. The modelling is made by an adjacency matrix in which each entity is a column and a line. Thanks to the matrix formalism automated computation treatments can be applied [START_REF] Browning | Applying the Design Structure Matrix to System Decomposition and Integration Problems:A Review and New Directions[END_REF].  Intermediate calculations of component features. They are calculated by the associated component blocks but the features are injected into other blocks. For example, the IGBT current rating is derived by the "IGBT block" ant it is provided to the "Inductances block".

The algebraic cycles implying several components are broken by the introduction of the design hypothesis. The N2D formalism clearly shows that the sizing process can be calculated in a sequential way. There is no data feedback below the diagonal which would generate algebraic loops. This important feature enables to decrease significantly the computation routine time assessing the system weight « ». The sizing models represented here by different blocks are mainly formed by estimation models. As displayed on the Figure 18 these models aim to define secondary characteristics (e.g weight of an electromechanical actuator) according to primary characteristics (e.g rated torque) [START_REF] Budinger | Estimation models for the preliminary design of electro-mechanical actuators[END_REF]. Different techniques can be applied in order to link primary and secondary characteristics:

 Component catalogue database. The classical and easiest way is to use databases. For example, large IGBT module database with searching filters can be used online in order to select the components fulfilling the requirements (e.g voltage and current rating) [23]. Nevertheless if the number of product references is too low this technique is not suitable. This is particularly the case in the aerospace field for which the components are often designed for a particular application.

 Detailed design models consist in designing the whole component: choosing the geometrical dimensions, the type of material, etc. The complexity depends on the study context : from a simple routine made of analytical equations to a model made of numerical simulations coupled with optimization algorithms [START_REF] Sareni | Niching Genetic Algorithms for Optimization in Electromagnetics, II. Shape Optimization of Electrodes using the CSM[END_REF].

 Scaling laws (also called similarity models) study the effects of varying representative parameters of a component compared with a known reference. This technique is used in different physical fields : electrical [START_REF] De | Integrated optimal design for power systems of more electrical aircraft[END_REF] [26], electro-mechanical [START_REF] Budinger | Estimation models for the preliminary design of electro-mechanical actuators[END_REF] [28], thermal [START_REF] Budinger | Scaling-law-based metamodels for the sizing of actuation systems[END_REF].

 Meta models and response surface approximations (also called surrogate models) is a mathematical method enabling to model analytically a phenomena or a system according to a limited number of references. The response surface approach is often used in order to obtain a polynomial equation modelling directly a complex system. This technique is widely spread in different engineering fields : electrical [START_REF] Brisset | Optimization with experimental design : an approach using Taguchi's Methodology and Finite Element Simulations[END_REF], automatics [START_REF] Faucher | Response surface methodology for the tuning of fuzzy controller dedicated to boost rectifier with power factor correction[END_REF], electro-mechanical [START_REF] Halabi | Optimal geometrical integration of electromechanical actuators[END_REF].

IV.2.2 Summary of the different techniques in the sizing process

The Table 1 provides an outlook of the different techniques employed for the estimation models necessary for our case study. For some components such as the IGBT modules and the inductances, several techniques are required in order to build the estimation models.

Database

Detailed design model

Scaling laws

Response surface The estimation model of the inductance is a good example representative of the techniques employed in order to obtain accurate and rapid sizing blocks.

1) IGBT modules X X X 2) Cold plates X 3) Capacitors X X 4) Inductances X X X X 5) Contactors X 6) Heat exchanger X 7) Chassis X

IV.2.1.1 Detailed design

In a first step, a detailed design model is defined in order to derive the geometrical dimensions of a toroidal shape mono-phase inductance as displayed in Figure 19. A computation routine is developed; it is based on 4 items:  Item 1: 8 high level parameters (Table 2): the self-inductance value, the inductance peak current, etc.

These parameters are given by the routine user;

 Item 2: a set of geometrical relationships derived by approach mixing scaling laws and the use of industrial database [34] [35];  Item 3: a set of physical relationships modelling the inductance (magnetics and thermal fields). These equations form the physical model;  Item 4: 2 design variables defined in the Table 3.

These 2 design variables are processed by a solver in order to find a feasible solution of the problem according to the data specified by the item 1, 2 and 3. The running time for solving the problem is around 1 second.

IV.2.1.2 Response surface

In a second step, the response surface technique is used. Several runs of the detailed design model are here necessary before to enable finding an analytical relationship linking the weight « », the peak current « » and the self-inductance value « »: = 1.45. 10 . . . .

Eq. 8

with the following unit choices: « : » ; « : μ » ; « : ».

By using the Eq.8 instead of running directly the detailed design model, the weight of the inductance can be determined nearly "instantaneously" in comparison with the 1 second range of the solver. A significant computation time gain is reached. The algorithm step 1 was applied on the industrial use case studied in this paper (part II.3). On a Quad Core 2.67 GHz machine, it takes around 1 hour to carry out the step 1.

The results are displayed on the Figure 24 If these 3 numbers are seen as an image of the solution physical complexity, thus the Figure 24 shows that the « , » decrease implies the increase of the solution complexity.

VII. Conclusion

Today the More Electrical Aircraft (MEA) trend is transforming step by step the aeronautical field. Among several technological bricks required for the MEA, power electronics is one of the most important technologies. With MEA, more and more power electronics modules need to supply high power systems consuming up to 100 kW (e.g the Environmental Control System called ECS). This significant change in terms of power demand requires to design innovate concepts for supplying of electrical loads requiring power electronic interfaces.

Traditionally aircraft loads requiring power electronics are supplied by a dedicated power electronics module (PEM). For such a structure, this paper has highlighted several drawbacks: reliability aspect (the PEM loss implies the load loss), cost aspect (the number of PEM references is equivalent to the load number), the PEM low utilization rate (if the load is intermittent). The Integrated Modular Power Electronic Cabinet (IMPEC) has been designed to tackle these drawbacks. IMPEC can be seen as the power electronics twin of the well-known IMA (Integrated Modular Avionics) that is the current embedded solution for processing electronics (28V DC).

IMPEC is made of several identical (modular) PEMs that can supply a set of electrical loads. A load can be supplied by one or several PEMs at a time. A matrix contactor enables to make the association « PEMload ». IMPEC offers flexibility to the operational needs by changing the association « PEM-load » from one operational case to another one.

In this context, the present paper has described a methodology for the IMPEC optimal design. To our knowledge, this is the first time that a complete design methodology is developed for such a system. As often in aeronautical applications, the design objective consists in minimizing the overall system weight. In addition, two design variables can be processed by the optimization algorithm: the number of PEMs and the reconfiguration of the contactor matrix. The characterization of the optimal design problem highlighted two main features: the multi-physical system sizing and the combinatorial explosion mainly due to the contactor matrix reconfiguration. The mixing of these two issues within a same problem is not usual and forms a complex problematic requiring innovative approaches.

Since designer shall cope with a combinatorial explosion, an explicit and rapid sizing process was developed. The process was cut into several sizing blocks: one for each component type. By introducing design hypothesis, all loops were deleted between the different sizing blocks. Within each sizing block, simple relations such as analytical equations were used. These equations can have been established beforehand by using high demanding processing methods such as solvers or simulations.

Since the classical combinatorial optimization algorithms were not sufficiently efficient due to the problem size, a dedicated method was developed. This heuristic was made of 2 steps: 1) pre-processing enabling to define beforehand a set of physical features of the solutions (PEM number, contactor number, inductance number) 2) the successive launch of a greedy algorithm and a local search algorithm in order to define the contactor matrix reconfiguration. The performance of this innovative algorithm was validated on problems of different sizes. In particular, the heuristic was capable to detect the optimal solution nearly instantaneously on reduced-size problems whereas several hours can have been required for a branch and bound algorithm.

Additionally, the complete heuristic results showed that -except for the minimal number of PEMwithin a set of solutions having the same number of PEMs, the lightest solution is the one minimizing the PEM power. In the future, this observation could be employed in order to speed up the detection of the more promising solutions.

Figure 1 Figure 2

 12 Figure1Sharing of the power-off take extracted from the engines (Airbus A330)[START_REF] Liscouët-Hanke | A model-based methodology for integrated preliminary sizing and analysis of aircraft power system architectures[END_REF]

Figure 3

 3 Figure 3 The 2 axis of the more electrical aircraft (MEA)

Figure 4

 4 Figure 4 Two classical structures to supply an electrical machine I.4 Non-simultaneous power demand peaks for the MEAThe analysis of the power demand of most consuming loads for the MEA, such as the ECS electrical air compressors, show important fluctuations during the flight mission. This trend is illustrated by the Figure5showing that the power demand profiles are different from one load to another : the respective load power demand peaks appear in different flight phases. For example, the maximum power demand of the ECS load n°1 appears in the flight phases n°1, 2, 4 and 21 (ground phases) whereas the ECS load n°2 work a full power during the phases n°5 and 20.

Figure 5

 5 Figure 5 Power demand of the 3 electrical air compressors according to the aircraft mission phases I.5 Paper structure

Figure 6

 6 Figure 6 Simplified schematics for the Boeing B787 electrical network

Figure 7

 7 Figure 7 Structure of an IMPEC implementation (on the left) -Example of configuration of the contactor matrix with loads n°1 and n°2 consuming respectively 60 and 20 kW (on the right)

Figure 8

 8 Figure 8 Left side : electrical schematic view of two PEM supplying one load each / Right side : functional view of this PEM-load association

Figure 9

 9 Figure 9 Examples of contactor matrix configuration for different PEM losses II.2.4 IMPEC advantage n°2 : flexibility according to power demands

Figure 10

 10 Figure 10 Reconfiguration of the contactor matrix in order to fit with the load power demands

Figure 11

 11 Figure 11 Use case for the assessment of the optimization design methods

Figure 12

 12 Figure 12 Graphical representation of the design variables with a selected configuration for the loading case n°1 III.1.2 Design objective As very often in the aeronautics industry, the design objective consists in minimizing the overall weight of the system. This is here modelled by: = min( ) = min ( . + . + . + + ) Eq 2

5 III. 2

 52 ∀ { , } ∈ {1, . . , } × {1, . , ℒ},, Two sources of complexity III.2.1 Complexity source n°1 : a multi-physical system design

Figure 13

 13 Figure 13 Main IMPEC components with identification the main physical fields

  , 2 different examples « 1 » et « 2 » are colored (green and purple).

Figure 14

 14 Figure 14 Potential contactor matrix configurations for each loading case and highlighting of 2 reconfiguration solutions

Figure 15

 15 Figure 15 Highlighting of the sizing process problematic IV.1 Axis 1 : the sizing process organisation IV.1.1 State of the art

Figure 16 ,

 16 Figure 16 Input/output block convention for a N2D and N2D example made of several blocks

Figure 17

 17 Figure 17 N2D of the sizing process

Figure 18

 18 Figure 18 Input/output of an estimation model

Figure 19

 19 Figure 19 Geometrical dimensions of a toroidal shape mono-phase inductance Parameter definition Notation Parameter definition Notation Self-inductance value « » Electrical frequency « » Inductance peak current « » Relative current ripple « / » Bus voltage « » Ambient temperature « » Switching frequency « » Thermal convection coefficient « ℎ »

Figure 22

 22 Figure 22 Step 2.A : Outcome of a greedy algorithm on problem made of 3 PEMs, 2 loads for a given minimal contactor matrix

  indicating the 3 main features computed by this data preparation step: « , », « , », « , » according to the 160 different values of « , » for the industrial use case. For each « , » value, 3 different points representing the 3 main features are plotted. For a solution having the highest « , » value, the data given by the step 1 are : « , = 5 », « , = 14 », « , = 0 » (PEM paralleling is not possible).

Figure 24

 24 Figure 24 Evolution of the minimal number of PEMs, contactors and inductances

Figure 26 ForFigure 27

 2627 Figure 26 For 81 solution performance combinations -Upper part : the mean difference to the optimal solution -Lower part: standard deviation to optimal solution.

Table 1

 1 Summary of the different estimation model techniques used in the sizing process IV.2.1 Sizing model example : the inductances

Table 2

 2 High level parameter of the detailed design model of the inductance

			75 kg		70 kg		72 kg
		PEM1	PEM2	PEM3	PEM1	PEM2	PEM3	PEM1	PEM2	PEM3
	Loading
	Case n°1	L1	L2		L1	L2		L1	L2
				1			2			3
	Loading	PEM1	PEM2	PEM3
	Case n°2	L1	L2	1	76 kg		INPUTS FOR THE GREEDY ALGORITHM
									30 A
									PEM1	PEM2	PEM3
	Loading Case n°3	L1 PEM1	L2 PEM2	PEM3					L1	L2	min MAT
				1	80 kg
	Loading	PEM1	PEM2	PEM3	PEM1	PEM2	PEM3	PEM1	PEM2	PEM3
	Case n°4	L1	L2		L1	L2		L1	L2
		13 A 11 A	1			2			3
			83 kg		86 kg		81 kg

The complete methodology introduced in this paper was employed to identify the most robust classes of solutions according to different technology characteristics: IGBT power losses level; weight of contactors, inductances and capacitors. The study showed that only 4 solution classes were optimal among 160 created by the heuristic. A design of experiments, applied on these 4 solution classes, quantified their linear sensibilities according to the technology characteristics. The results have highlighted that the weight of contactors and inductances were clearly the main drivers enabling to bring weight savings in the future. Eventually, the optimization of the component sizing (in particular these commutation components) is the main perspective of this work.

Design variables

Value range Notation

Maximal induction in the tore 0 ; 0.5 « » Winding fill factor 0 ; 0.5 « »

Table 3 Design variables of the detailed design model of the inductance V. Optimization algorithm: a dedicated heuristic

Exact methods such as the branch and bound algorithm [START_REF] Clausen | Branch and Bound Algorithms -Principles and Examples[END_REF] cannot be applied due to the size of the problem (combinatorial explosion). For the same reason, meta-heuristics (Tabu search [START_REF] Hertz | Tabu search for large scale timetabling problems[END_REF], Simulated annealing [START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF], etc) usually use in combinatorial optimization problems are inefficient. Therefore a dedicated algorithm is developed for the IMPEC optimal design problem. This algorithm uses the evaluation function introduced here before each time a system weight assessment is required.

V.1 Presentation of the heuristic algorithm

As illustrated by the Figure 21, the heuristic is made of 2 successive steps.

V.1.1 Step n°1 : data preparation

The step 1 consists in a pre-processing computation enabling to define several physical features of the system for different PEM current levels « , ». This current is the maximum current level seen by the PEM semi-conductors over all the loading cases « ». The value of « All these data will be used as inputs for the heuristic step 2. These physical characteristics and in particular the contactor matrix definition enable to drastically reduce the combinatorial space for the definition of the reconfiguration solution « » (connexion of load/PEM for all the loading cases) without excluding the most promising solutions.

V. 1.2 Step n°2 : optimization

The step 2 is the core part of the optimization algorithm. It defines the design variables « , » forming the reconfiguration solution « ». From a value of « , » and for the associated set of physical features, the step 2 is made of 2 successive steps:

A intends to build a solution having the minimal number of contactors « , » by using the minimal contactor matrix « ». For each loading case « », a matrix configuration is found by the greedy algorithm in order to minimize the system weight « » as specified in the part III.1.2. The greedy algorithm is a typical method used in combinatorial optimization [START_REF] Garcia-Martinez | Tabu-enhanced iterated greedy algorithm: A case study in the quadratic multiple knapsack problem[END_REF]. The simplicity and rapidity of this algorithm enables to cope with the combinatorial explosion for the choice of the reconfiguration solution. This is qualified as a construction algorithm since the solution is built at each step of the process. In our application, it means that at each step of the algorithm, the matrix configuration of one loading case is decided. This choice is not re-challenged during the rest of the exploration of the greedy algorithm.

The Figure 22 shows the type of outcome that could produce the greedy algorithm on a simple problem made of 4 loading cases, 2 PEMs and 3 loads. In this illustrative problem the minimal contactor matrix « » is given on the left-hand side of the Figure 22. The PEM is rated as « , = 30 ». The first step of the method consists in assessing the 3 different configurations for the loading case n°1. The configuration leading to minimum weight is taken. Here it is the configuration n°2. This choice will not be re-challenged in the greedy algorithm process. Then the second iteration of the algorithm goes to the case n°2 for which there is only one choice due to « » and the loads' power demand. The same situation appears for the loading case n°3. Then the algorithm finishes with the loading case n°4 for which there are again 3 different configurations. The choice n°3 is taken since it minimizes the weight.

The greedy algorithm ends up with the solution « 1 = ( 1, 1) » (see Figure 21) . On our illustrative example, the solution « 1 » has a weight of 81 kg. The matrix reconfiguration solution « 1 » consists of the configuration highlighted in red on the Figure 22. The physical constitution « 1 » of the solution « 1 » is 3 PEMs, 4 contactors and 3 inductances (one at each PEM output -since each PEM is parallelized during the loading case n°2 and n°3).

If the solution « 1 » has the minimum number of inductances « , », the step 2 is over. Otherwise the step 2.B is launched. In our illustrative example, the minimum number of inductances is , = 2 » therefore the step 2.B shall be launched.



Step 2.B intends to modify locally « 1 » in order to obtain a solution having a minimal number of inductances. Only loading cases implying PEM paralleling are modified in order to delete one or several inductances. During this step, several solutions « » (with « ≥ 2 ») can be produced. They are the trade-off between solutions having a minimal number of inductances « , » and those having a minimal number of contactors « , ».

In the illustrative example, a possible modification brought by the step 2.B would have been to modify the configuration of the loading case n°3 as illustrated on the Figure 23. The load n°1 is now connected to the PEM3 and the load n°2 is supplied by PEM1 and PEM2. If it is the only modification brought to « 1 », the new matrix reconfiguration solution « 2 » will lead to delete the inductance ». Additionally, it is interesting to notice that even though the least complex and most complex solutions have very different physical constitutions, their weight are nearly similar.

Except for the class with minimal number of PEMs (i.e. 5 in this example), Further analysis of solutions displayed on the Figure 25 shows that, given a class of solutions (i.e. set of solutions having the number of PEM), the best one is provided by minimizing « , ». For example among the class of solutions having 7 PEMs, the lightest solution is the one having the minimal PEM power: « , = 2.0 ». PEMs (up to 7), the problems can be exactly solved by the Branch and Bound algorithm [START_REF] Clausen | Branch and Bound Algorithms -Principles and Examples[END_REF].

Although the problems can be qualified as "reduced-size", the combinatorial space can go up to « ≈ 10 » solutions.

The performance assessments show that the heuristic is capable to find the optimal solution for nearly all the problems (except one). The heuristic computation time is extremely fast (several seconds) whereas the Branch and Bound can take up to 5 hours to find the optimal solution.

 On industrial problems presented in the part II.3 but with a limited number of PEMs (up to 7), a lower bound values can be identified (a dedicated algorithm, not explained in this paper, was developed in order to find a lower bound value).

The assessments show that the heuristic is on average at 6 % of the lower bound. This means that in the most unfavorable situation, the heuristic produces a solution 6 % away from the optimum.

 On industrial problems presented in the part II.3 with a number of PEMs going from 5 to 12. A reference solution is defined by applying the rules of thumb usually applied on the industrial field.

The assessments show that the heuristic is capable to find on average a solution offering 25 % weight savings according to the reference solution.

VI. Post-optimal analysis As technologies are evolving rapidly and as there are always uncertainties on the sizing models, cartography is carried out by solving the problem with different combinations of component characteristics:  the IGBT power losses ;  the inductance weight ;  the contactor weight ;  the capacitor weight.

In addition, 3 levels (low, nominal & high) are defined for each component characteristic. For instance, if the "nominal" value of the contactor weight characteristics is « 200 / ». The "low" (pessimistic performance view) and "high" (optimistic performance view) values are respectively « 100 / » and « 300 / ».

In total, 81 (3 ) combinations are set. In order to have a synthetic view over these 81 runs, the Figure 26 is built. On each graph, the 160 solution classes are displayed.

The upper part of the graph shows for each solution class the mean difference to the optimal solution. For instance, the solution class n°53 defined by the couple: « , = 2 ; = 7 » is on average at 5 % of the optimal solution for the 81 component characteristic combinations. The lower part of graph of the Figure 26 shows the standard deviation to the optimal solution. Both graphs clearly show that the most robust solution classes are between « = 7 » and « = 9 » with a PEM power between « Finally, it is interesting to notice that over the 81 runs, only 4 solution classes are optimal. They are highlighted in red color (the number of times in which they are optimal is written just below). Thus the solution class n°1 defined by the couple: « , = 4.7 ; = 5 » is the best solution for 17 out of 81 runs (combination of component characteristics).