Sonia Cafieri
email: sonia.cafieri@enac.fr

Leo Liberti
email: liberti@lix.polytechnique.fr

Frédéric Messine
email: messine@n7.fr

Bertrand Nogarede
email: nogarede@n7.fr

Optimal Design of Electrical Machines: Mathematical Programming Formulations *

Keywords: analytical model, formulation, modeling, local optimization, inverse problem, design, electrical machine

The optimal design of electrical machines can be mathematically modeled as a (mixed-integer) nonlinear optimization problem. We investigate the impact of different mathematical formulations on the results obtained using a local optimization solver which is well-known in the engineering community: MatLab's fmincon function. Our analysis is based on six different mathematically equivalent formulations for the same problem of the design of an electrical machine without slot. Our results underline the important impact that formulation differences may have on solver performance even on a small example of design.

Introduction

The design of electromechanical actuators is known as an inverse problem, i.e. from the characteristic values given by the schedule of conditions (for example the torque), obtain the structure, the dimensions and the material compositions of the actuator constitutive parts, [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF][START_REF] Fitan | The electromagnetical actuators design problem: a general and rational approach[END_REF][START_REF] Neittaanmki | Inverse Problems and Optimal Design in Electricity and Magnetism[END_REF]. One is usually interested in performing an optimal design where a given criterion is optimized (for example, the volume of the magnet is minimized). The interest of the electromagnetic actuators design combining optimization algorithms and analytical models has in fact already been widely shown in the literature [START_REF] Messine | A Deterministic Global Optimization Algorithm for Design Problems[END_REF][START_REF] Messine | Deterministic global optimization using interval constraint propagation techniques[END_REF][START_REF] Huang | A direct approach to electrical machine performance evaluation: Torque density assessment and sizing optimisation[END_REF][START_REF] Nogarede | Optimal design of permanent-magnet machines using analytical field modeling[END_REF][START_REF] Rao | Mean field annealing (mfa) and optimal design of electromagnetic devices[END_REF][START_REF] Slemon | Modeling and design optimization of permanent magnet motors[END_REF][START_REF] Wurtz | A methodology and a tool for the computer-aided design with constraints of electrical devices[END_REF]. The inverse problems of electromechanical actuator design are more general than optimal dimensioning problems [START_REF] Fitan | The electromagnetical actuators design problem: a general and rational approach[END_REF]. These problems must be formulated as mixed-constrained optimization problems, see [START_REF] Fitan | A general analytical model of electrical permanent magnet machine dedicated to optimal design[END_REF][START_REF] Fitan | The electromagnetical actuators design problem: a general and rational approach[END_REF] for complete formulations. Optimization problems of this kind are in general quite difficult to solve to global optimality. Exact global optimization solvers can however be adapted for their solution. Methods based for example on interval analysis have been proposed in [START_REF] Fitan | The electromagnetical actuators design problem: a general and rational approach[END_REF][START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF][START_REF] Messine | Deterministic global optimization using interval constraint propagation techniques[END_REF][START_REF] Messine | A Deterministic Global Optimization Algorithm for Design Problems[END_REF]. Other recent methods for global optimization of mixed-integer problems can be applied as well. In practice, however, the picture is not as simple: recent powerful solvers are not always publically available (this is the case for interval analysis based solvers) or do not interface to a modelling environment for easy using. In practice, in the context of engineering applications, one usually prefers to resort to well known and easy to use optimization solvers [START_REF] Neittaanmki | Inverse Problems and Optimal Design in Electricity and Magnetism[END_REF][START_REF] Nogarede | Optimal design of permanent-magnet machines using analytical field modeling[END_REF][START_REF] Rao | Mean field annealing (mfa) and optimal design of electromagnetic devices[END_REF][START_REF] Wurtz | A methodology and a tool for the computer-aided design with constraints of electrical devices[END_REF]. In this sense, a very successful solver is provided by the MatLab's fmincon function [8]. It is based on Active-Set and Sequential-Quadratic Programming (SQP) methods with computation of the Hessian from Quasi-Newton techniques [START_REF] Gill | Procedures for optimization problems with a mixture of bounds and general linear constraints[END_REF][START_REF] Gill | Numerical Linear Algebra and Optimization[END_REF] and finds local optima of nonconvex nonlinear optimization problems with continuous decision variables. This solver does not represent the state of the art for nonlinear nonconvex optimization and, being a local solver, does not provide global solutions. However, it usually provides good quality (local) solutions which are often considered satisfactory in practical contexts.

Even though this kind of solver is quite easy to use, the mathematical features of the optimization problem may have an impact on the practical efficiency of the solver. It is known that the same optimization problem can be often formulated in different ways, sharing some properties such as optima and feasible regions. Furthermore, the formal description of optimization problems has an impact on the applicability and efficiency of the corresponding solution methods. Indeed, the study of reformulations is an active research area in the optimization community [START_REF] Liberti | Reformulations in mathematical programming: Definitions and systematics[END_REF][START_REF] Liberti | Reformulations in Mathematical Programming: a Computational Approach[END_REF]. The aim of this paper is to investigate the efficiency and reliability of standard local optimization solvers when handling different mathematical formulations. We illustrate the impact of such different formulations on solver performance. While providing any advancement on the global solution of engineering design problems is out of the scope of this paper, we consider a local solver which is widely used in the engineering community for optimal design, with the aim of providing guidelines for designers in practical engineering applications.

Our application testbed is a slotless electrical rotating permanent magnet machine. This example was first presented in [START_REF] Nogarede | Optimal design of permanent-magnet machines using analytical field modeling[END_REF] and also studied in a lot of papers such as [START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF][START_REF] Nogarede | Optimal design of permanent-magnet machines using analytical field modeling[END_REF][START_REF] Messine | Deterministic global optimization using interval constraint propagation techniques[END_REF][START_REF] Wurtz | A methodology and a tool for the computer-aided design with constraints of electrical devices[END_REF]. The analytical equations of the considered design problem, which come from approximations of Maxwell's equations (taken in the quasi-static mode) and mechanical considerations, can be found in [START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF][START_REF] Miller | Brushless Permanent-Magnet and Reluctance Motor Drives[END_REF][START_REF] Nogarede | Optimal design of permanent-magnet machines using analytical field modeling[END_REF]. These equations are recalled below:

Γ em = π 2λ (1 -K f ) k r βE ch ED 2 (D + E)B e (1)
E ch = AJ cu = k r EJ 2 cu (2) p = πD ∆ p (3)
K f = 1.5pβ e + E D (4)
C = πβB e 4pB iron D (5)
B e = 2l a P D ln D+2E D-2(la+e) [START_REF] Liberti | Reformulations in mathematical programming: Definitions and systematics[END_REF] where Γ em is the electromagnetic torque (from an energetic calculation, see Eq. ( 1)); D(m) is the bore diameter, λ the diameter over length ratio, E(m) the winding thickness, β the polar arc factor, k r a coefficient of occupation; the global heating up of the winding is rather roughly modeled by E ch (function of current electric loading A and J cu , see Eq. ( 2)); J cu (A/m 2 ) is the current areal density; p is the number of pole pairs, ∆ p the polar step (p is linked to ∆ p by Eq. ( 3)); K f is a semi-empiric magnetic leakage coefficient (established by numerical simulations, see Eq. ( 4)); e(m) is the thickness of the mechanical air-gap; C(m) is the thickness of yoke (Eq. ( 5) is obtained by neglecting interpolar leakages and armature reaction flux); B iron is the magnetic field in the iron; B e (T ) is the no-load magnetic radial flux density (see Eq. ( 6)); l a (m) is the thickness of the permanent magnets, P the magnetic polarization. For this study, we fix Γ em = 10N.m, P = 0.9T , k r = 0.7, B iron = 1.5T , E ch = 10 11 A/m and the polar step ∆ p = 0.1m as in [START_REF] Nogarede | Optimal design of permanent-magnet machines using analytical field modeling[END_REF]. The other parameters can vary inside the following intervals:

D(m) ∈ [0.01, 0.5], λ ∈ [1, 2.5], l a (m) ∈ [0.003, 0.05], E(m) ∈ [0.001, 0.05], C(m) ∈ [0.001, 0.05], β ∈ [0.8, 1], B e (T ) ∈ [0.1, 1], J cu (A/m 2 ) ∈ [10 5 , 10 7 ], K f ∈ [0.01, 0.3], e(m) ∈ [0.001, 0.005] and p ∈ {1, • • • , 10}.
The motor structure is presented in Figure 1. The addressed design problem can be formulated as an optimization problem by minimizing or maximizing a given criterion under constraints given by equations ( 1)- [START_REF] Liberti | Reformulations in mathematical programming: Definitions and systematics[END_REF].

Note that, even though we focus on a simple specific design problem for the sake of illustration, the same analysis may be generalized to more complex examples formulated as nonlinear nonconvex constrained problems.

The rest of this paper is organized as follows. In Section 2, we propose six equivalent mathematical formulations of the optimal design problem of a slotless electrical rotating permanent magnet machine. In Section 3, we first computationally compare the six proposed formulations in the continuous case, i.e., with a fixed value for p and we discuss about their efficiency when the optimization problem is solved using the solver fmincon of MatLab. We show that numerical performances are different depending on the formulation. The the effect of changing starting points on the optimization results is also discussed. We then summarize some results when p is free on three mixed-integer formulations. Some concluding remarks are given in Section 4.

Mathematical formulations

In this paper, the addressed optimization problem has the following general form:

P :                        min x∈IR n y∈IR m f (x, y) s.t. g i (x, y) ≤ 0, ∀i ∈ {1, • • • , p}, h i (x, y) = 0, ∀i ∈ {1, • • • , q}, y i = A i (x, y J i ), ∀i ∈ {1, • • • , m}, x i ≤ x i ≤ x i , ∀i ∈ {1, • • • , n}, y i ≤ y i ≤ y i , ∀i ∈ {1, • • • , m}.
where

J i ⊆ {1, • • • , m} \ {i}.
Moreover, y i depends explicitly or implicitly on x by recursive calls to A j functions and there is no cycle in the definition of y i . Hence,

y i = A i (x, y J i ) = A R i (x)
; as a vectorial notation, we use: y = A R (x). By replacing the occurrences of y i in (P) by A R i (x), we obtain the following reformulation: R :

                 min x∈IR n f (x, A R (x)) s.t. g i (x, A R (x)) ≤ 0, ∀i ∈ {1, • • • , p}, h i (x, A R (x)) = 0, ∀i ∈ {1, • • • , q}, x i ≤ x i ≤ x i , ∀i ∈ {1, • • • , n}, y i ≤ A R i (x) ≤ y i , ∀i ∈ {1, • • • , m}.
This reformulated problem (R) of (P) is obtained by removing m variables (all the variables y) and changing m equality constraints (y i = A i (x, y J i )) to 2m inequality ones (y i -A R i (x) ≤ 0 and A R i (x) -y i ≤ 0). The two optimization problems (P) and (R) are mathematically equivalent because they provide the same solution which corresponds to the global minimum of the two problems: all realizable solutions of (P) are realizable solutions of (R) and reciprocally, moreover the objective function values are equal (because y = A R (x)).

From equations of the design problem such as those defined in (1)-( 6), and by introducing the minimization of the volume of the magnets: V m = πβl a D λ (D-2e-l a ), we can generate some distinct but equivalent formulations of our design problem, as described above by problem (R). The differences concern some variables which are replaced by functions and equality constraints by inequality ones. This yields six formulations.

The first formulation comes directly from the equations ( 1)-( 6) of the design problem. Hence, we obtain:

F 1 :          min (D,λ,••• )∈D⊂I R 10 p∈{1,••• ,10} V m (D, λ, ...) = πβl a D λ (D -2e -l a ) s.t.
constraints defined by Equations ( 1)-( 6)

Formulation F1 is a nonlinear, nonconvex optimization problem, the nonlinearities arising in the analytical formulae expressing the objective function and the constraints. The variables D(m), λ, l a (m), E(m), C(m), β, B e (T ), J cu , K f and e(m) are continuous, while p is integer. Hence, the problem is a mixed-integer nonlinear program (MINLP). We remark that the problem is badly scaled, because some parameters such as J cu have large values and the others, such as D, e, . . ., are definitely small. However, this seems to have no impact when a solver like fmincon is used on all the formulations considered in this paper. We also note that the integer variable p will be fixed to a constant in the first part of Section 3 to provide a continuous nonlinear problem (NLP) and p will be an integer variable in subsection 3.4.

The following formulations are reformulations of F1 which are mathematically equivalent.

First, we consider B e as an auxiliary function (depending on D, l a , E, e) which returns

2laP D ln( D+2E D-2(la+e) )
and whose value is bounded by 0.1 and 1; in fact, in all the expressions of Formulation F1, B e could be directly replaced by

2laP D ln( D+2E D-2(la+e) )
. Thus, from F1 we remove one variable and its corresponding equality constraint, and we replace this constraint with two inequality constraints which impose lower and upper bounds on the value of B e . We obtain the following formulation:

F 2 :          min (D,λ,••• )∈D⊂I R 9 p∈{1,••• ,10} V m (D, λ, ...) = πβl a D λ (D -2e -l a )
s.t. constraints defined by Equations ( 1)-( 5)

0.1 ≤ B e (.) ≤ 1
where B e (.) is an auxiliary function and returns

2laP D ln( D+2E D-2(la+e) )
(or where B e (.) is directly replaced in Formulation F2 by its corresponding expression).

The two formulations only differ in size (number of variables and constraints). The first one, F1, has 10 continuous variables: D, λ, l a , E, C, β, J cu , K f , e, B e and 6 equality constraints, while the second one, F2, has 9 continuous variables: D, λ, l a , E, C, β, J cu , K f , e, 5 equality constraints and 2 inequality constraints involving B e . Additional constraints are given in both cases by the bounds on the variables. We remark that any solution (global minimum) of F1 is also a solution of F2 and conversely.

We now proceed to eliminate another variable and considering it as a function. Again, the equality constraint corresponding to the selected variable is replaced by two inequality constraints which impose lower and upper bounds on its value. Considering C as a function which returns πβBe 4pB iron D, we have:

F 3 :              min (D,λ,••• )∈D⊂I R 8 p∈{1,••• ,10} V m (D, λ, ...) = πβl a D λ (D -2e -l a )
s.t. constraints defined by Equations ( 1)-( 4)

0.1 ≤ B e (.) ≤ 1 0.001 ≤ C(.) ≤ 0.05
Formulation F3 has 8 continuous variables, 4 equality constraints, 4 inequality constraints and bounds on variables.

Considering now K f as a function which returns 1.5pβ e+E D , we obtain the following formulation:

F 4 :                  min (D,λ,••• )∈D⊂I R 7 p∈{1,••• ,10} V m (D, λ, ...) = πβl a D λ (D -2e -l a )
s.t. constraints defined by Equations ( 1)

-(3) 0.1 ≤ B e (.) ≤ 1 0.001 ≤ C(.) ≤ 0.05 0.01 ≤ K f (.) ≤ 0.3
Formulation F4 has 7 continuous variables, 3 equality constraints, 6 inequality constraints and bounds on variables. Another possibility to reformulate the problem is to introduce a new variable y which replace a nonlinear term appearing in the first equality constraint, and add the corresponding constraint to the formulation. Let y be equal to √ βE. We can reformulate the problem as follows:

F 5 :                              min (D,λ,••• )∈D⊂I R 8 p∈{1,••• ,10} V m (D, λ, ...) = πβl a D λ (D -2e -l a ) s.t. Γ em = π 2λ (1 -K f (.)) √ k r E ch yD 2 × (D + E)B e (.) y 2 = βE constraints defined by Equations (2),(3) 0.1 ≤ B e (.) ≤ 1 0.001 ≤ C(.) ≤ 0.05 0.01 ≤ K f (.) ≤ 0.3
Note that now the square root term appearing in the first constraint is a constant term. Note also that y has a positive value because of the bounds on variables β and E. Formulation F5 has 8 continuous variables, 4 equality constraints, 6 inequality constraints and bounds on variables.

Finally, we consider again y = √ βE and we add this variable to formulation F1, obtaining a new formulation with 11 continuous variables and 7 equality constraints:

F 6 :                  min (D,λ,••• )∈D⊂I R 11 p∈{1,••• ,10} V m (D, λ, ...) = πβl a D λ (D -2e -l a ) s.t. Γ em = π 2λ (1 -K f ) √ k r E ch yD 2 (D + E)B e E ch = AJ cu = k r EJ 2 cu y 2 = βE constraints defined by Equations (2)-(6)
These six formulations vary in size from 7 to 11 continuous variables, from 3 to 7 equality constraints and from 0 to 6 inequality constraints.

Comparing the six formulations above, we note that when the number of variables increases, the nonlinearities of the equations decrease yielding simpler optimization problems but with more variables. Since the solver performance is roughly directly proportional to both the number of variables and the number of nonlinearities in the objective and constraints, a natural trade-off situation arises. Note that this behavior is independent of the specific considered problem and the same considerations apply to other possibly more general design problems having nonlinear optimization formulations like (P).

In the following section, we numerically compare the proposed formulations and discuss about their efficiency when a local optimization solver is used to solve the optimal design problem. We also consider two other objective functions for the addressed problem:

V u = π D λ (D + E -e -l a )(2C + l a + e + l a ) (7)
P j = πρ cu D λ (D + E)E ch (8)
where V u represents the volume of the active parts and P j the losses by joule effects with ρ cu = 0.018 × 10 -6 .

Computational comparison between formulations

Different mathematical formulations of a given problem may share the same properties, such as for example feasible region and optima. However, some of them are easier to solve than others. The same optimization solver can perform differently depending on the formulation. For this reason, a considerable amount of work is often devoted to investigate efficient reformulations, see e.g. [START_REF] Liberti | Reformulations in Mathematical Programming: a Computational Approach[END_REF].

In this section, we computationally compare formulations F1 to F6 and evaluate their impact on the performance of the deterministic local optimization solver fmincon of MatLab v7. In a first series of experiments, the integer variable p is fixed to the constant value 5 (known to be optimal), so that formulations F1 to F6 become continuous optimization problems. We employ MatLab's optimset function to fix the following parameters for fmincon: (i) the maximum number of iterations and of function evaluations is fixed to 30000, (ii) the tolerance on the value of the function is fixed to 10 -10 , (iii) the tolerance on the solution point is fixed to 10 -6 , (iv) the tolerance on the constraints satisfaction is fixed to 10 -7 and (v) FunValCheck is fixed to 'on', which generates an error if values become complex during the computation (this can occur in equations ( 1)-( 6) when computing the logarithm of a negative number). In this last case (v), we have to use MatLab's exception handling mechanism (the try and catch functions). The considered tolerances correspond to standard settings in optimization solvers.

The performance of a local method such as fmincon often depends on the choice of a starting point (values initially assigned to the decision variables). This may be given by the existing design configuration in case one exists, otherwise it is very difficult to choose it appropriately. Depending on this choice, we have three possibilities: (i) the algorithm converges to a local solution; (ii) the algorithm progresses slowly towards a local optimum, but finally exceeds the maximum allowed number of iterations or function evaluations; (iii) the algorithm fails for a number of other reasons, with no meaningful answer. An aprioristic choice of starting point guaranteeing the occurrence of case (i) is usually impossible. We try and make our results "starting point independent" by implementing a simple multistart approach from 1000 randomly generated starting points (random values are generated using MatLab's rand function). We then record the percentage of starting points yielding a local and a global optimum.

To generate starting points, two strategies are possible. In the first one, random values are generated for each variable of Formulation F6 (corresponding to the formulation having the maximum number of variables) and the same values are used for all Fi, i = 1, . . . , 5. In the second one, a random value is generated for the 7 continuous variables of Formulation F4 (corresponding to the formulation having the minimum number of variables) and corresponding values are computed for problem entities which are treated as functions (in formulations F1 to F3, F5, F6). For example, in Formulation F2 the value of C and K f can be computed using auxiliary functions defined for Formulation F4 on the basis of the values of D, λ ... which are randomly generated.

Numerical comparisons using V m

In this subsection, we compare the results of our Multistart method on the six optimization problems F1 to F6 defined in Section 2. The value of p is fixed to 5.

Results are shown in Tables 1 and2 depending on whether the starting points are generated from Formulation F6 or from Formulation F4 respectively (the same base of the 1000 starting points is used for all the formulations). We compare, for each formulation, the percentage of local minima and the percentage of best local minima found, the best and the worst values of local minimum found, the best CPU time, the average CPU time, the average CPU time corresponding to successful running (i.e., providing a local minimum), and the worst CPU time. It appears that the percentage of failure is quite high. Indeed, for all the considered formulations, the percentage of local minima found is always lower than 50% and slightly higher than this percentage in one case only. On the last line of the tables (x0 = mid), we report the local minima which are obtained (or not) using the middle point of the variable bounds taken as a starting point; '-' denotes no local minima. In Table 1, only two formulations provide a local minimum using the middle of the bounds as a starting point compared to Table 2 where four formulations achieve the convergence. Nevertheless, all these local minima are far from the best one (one order of magnitude of difference). This confirms that a lot of care must be paid in the choice of a good starting point. When some values are computed starting from values assigned to variables on which they depend (Table 2), the percentage of local minima found is higher. In particular, the formulations with more variables, F1 and F6, provide the highest percentage of local minima found. This remark yields to show that the formulations with more variables F1 and F6, where the nonlinearity decrease, associated with starting points made from components which are partially randomly generated and partially computed from other components using auxilary functions, provide the most efficient Multistart technique: about 50% of chance to provide a local minimum compared to less than 45% for all the 10 others (see Tables 1 and2), and about 5% to find the best local minima compared to less than 4.6% for the other formulations excepted F5; note that a minimum is considered equal to the best one if their corresponding optimal value are numerically close.

CPU times are in general practically acceptable (lower than 30 seconds). The average CPU time is however increased by the time spent in unsuccessful runs, where convergence is not achieved because of a worse choice of the starting point. We note that there is no impact on the CPU-time if the formulation has more variables than another one. In Tables 3 and4 we provide the best found solutions for each formulation. We remark that Formulation F4 produces the minimal value for the objective function V m ; this is understandable because F4 is the most compact formulation. We note also that all computed objective function values are different except for F2 and F3 in Table 3, this minimum point was also found using formulations F1, F2 and F6 in Table 4. Thus, as we found 8 different minima, we can conclude that there exist a lot of solutions which provide close values for the objective function V m . This kind of Multistart method could be adapted to provide a large set of minimal points. Note that Formulation F4 provides the same results in the two tables 3 and 4 because we take exactly the same starting point in both cases, as it is the most compact formulation, so there are no components of the starting point that can be deduced from others.

In [START_REF] Nogarede | Optimal design of permanent-magnet machines using analytical field modeling[END_REF], the considered design problem with formulation F1 was first modeled and solved using a local search algorithm. These first solutions were clearly local minima because in [START_REF] Wurtz | A methodology and a tool for the computer-aided design with constraints of electrical devices[END_REF] new better solutions were found. In [START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF], a deterministic global optimization solver was used for the first time (with a restriction to 4 digits on the accuracy of parameters). The objective function values in Table 4 are the same as the one corresponding to the global solution provided in [START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF]. Thus, some global minimizers are found here using a multistart technique with the fmincon procedure. The same remarks could be done for the results involving the objective functions V a and P j presented in next sections.

More numerical comparisons using V u and P j

In order to summarize all the computations that we performed, we present in the following and until the end of the paper, only the two most important formulations: F1, which comes directly from equations ( 1)-( 6), and F4, which is the most compact one. Moreover, two distinct techniques are applied to F1 to generate its starting point: F1(rnd) indicates that all the components of the starting point are randomly generate and F1(cmp) indicates that the starting point is partially randomly generated and some components are computed from others. In the following tables, the best and worst values of local minima and the best and worst CPU-times used for one local minimization step of the Multistart are reported as in the previous subsection.

In this subsection, two different objective function are minimized: V u and P j defined in [START_REF] Liberti | Reformulations in Mathematical Programming: a Computational Approach[END_REF] and (8).

In Table 5, results using the Multistart method with the objective function V u (7) are given. Similarly, in Table 6, results for the objective function P j (8) are presented. The best found local are given in Tables 7 for V u and 8 for P j .

We observe that results are similar to the ones of the previous subsection. This seems to confirm that the Formulation F1 associated with partially randomized starting points (denoted by F1(cmp)) provides the best results in terms of percentage to converge to a local solution, while the percentage to find the best local minimum are still low for V m but quite high (about 20%) for P j . In these two cases too, the minimum points are all different.

These remarks show that the objective function does not have a deep impact in the considered optimal design problems. Note that only Formulation F1(cmp) provides a local minimum when the middle of the bounds on the variables is used as a starting point (and some components are computed from the others). Note also that these two local minima are not the best minima found (in Table 5, this local minimum is close to the worst case). 3) and of percentage of times when the best local minimum is found (Figures 4 and5), using totally or partially randomized starting points.

From Figure 2, we can see that Formulation F4 (which is the most compact one) provides the best results for V m and V u , and Formulation F3 provides the best results for P j . Comparing the six formulations on each of the considered objective functions, it appears that close results were obtained for formulations F2, F3, F4 and F5, while formulations F1 and F6 provide a percentage which is about 3% lower. Contrarily on Figure 3 (when some components of the starting points are computed from others), F1 and F6 appear clearly to be the most efficient formulations for our Multistart method for all the three objective functions, providing a percentage about 5% higher than all other formulations.

Concerning the percentage of best found local minima, from Figures 4 and5 we can see that similar results are obtained for all the formulations and all the objective functions.

Computational comparison with p free

In the considered optimal design problem, the number of pole pairs p has an integer value and can be considered as a variable of the problem. In fact, algorithms such as the one implemented in fmincon of MatLab, need continuous formulations and twice differentiable functions. In some papers of the literature (see [START_REF] Nogarede | Optimal design of permanent-magnet machines using analytical field modeling[END_REF][START_REF] Wurtz | A methodology and a tool for the computer-aided design with constraints of electrical devices[END_REF]), p is considered as a continuous variable and the user has to convert the obtained real value of p into an integer one. Previously, we considered p as a fixed parameter (p = 5). Some ways to deal with the integer variable p could be the following: (i) p is converted into a real variable p ∈ [1, 10] and its optimal value is rounded into an integer one by the user (all the other variables must be adjusted); (ii) p is treated as a real variable but an equality constraint (p -1) In this study, we focus on two formulations only, the first one where all entities are considered as variables of the problem F1 and the one where there is the minimum number of variables F4. Results with F1 are again obtained using random generated values for all variables as well as random values for the 7 variables of F4 and values computed from this 7 values for the other variables (as we did in subsection 3.2). Since, as noticed in the previous section, no relevant differences in the behavior of the considered formulations can be observed using a different objective function like V u or P j , we focus in the following only on the minimization of the volume V m of the magnet.

× (p - 2) × • • • × (p -
In Figures 6 and7, we compare the formulations in terms of percentage of found local minima and percentage of best found local minima respectively, considering the different choices of p discussed above. Thus, these figures summarize the results obtained via the multistart method associated with fmincon. From Figure 6, one can see that Formulation F1 with recomputed values for the starting point appears to be the most efficient to find the largest number of local solutions, with the different strategies to deal with p. Thus, as the values of the percentage of best found local solution are low, we can have more confidence on the local solutions found by this formulation in order to provide the global solution.

It is also evident, from Figure 6 and specially from Figure 7, that considering p as a variable, as expected the problem becomes more difficult to solve and consequently the percentages of solutions found and best found solution are low. The fourth proposed strategy, based on a lower-degree polynomial constraint, allows to obtain higher percentages, though very low, when p is treated as an integer variable.

In Figure 7, it is very clear that, using Formulation F1 with some computed values for the starting point, a higher number of starting points allows to obtain the best found local solution. Thus, as the values of the percentage of best found local solution are low (except for p ∈ {1, • • • , 10} but the minimal value for F1(cmp) is the lowest one), we can have more confidence on the local solutions found by this formulation in order to provide the global solution. We have presented in this paper six mathematical formulations of an optimal design problem. We have shown that, even though the formulations are mathematically equivalent, their numerical performances are different when an optimization solver is used. We mainly discussed the impact of reformulations on a classical local solver based on Active-Set/SQP/Quasi-Newton algorithm, thus showing that the designer must take care about the formulation of the problem in order to make more efficient the use of this kind of local algorithms, such as the one proposed by MatLab in fmincon.

First, we have discussed the solutions found using the six formulations, considering the problem as a continuous one (fixing the value of the integer variable present in the formulation), in a multistart setting where the solution is computed at each step by fmincon. The first formulation F1, with 10 continuous variables, provides the best percentage of local solutions found, specially when the starting point is computed by giving random values to some of its components and computing the other values from these ones. Second, we have discussed about strategies to deal with integer variables, such as the number of pole pairs. Again Formulation F1 (with starting points partially generated) appears to be the most efficient in providing the best local solution.

These remarks could be extended to more general design problems when a local solver is used: formulations of type (P), with starting points only depending on variables x, seem to be more easier to solve and provide better results than their possible reformulations of type (R).

Figure 1 :

 1 Figure 1: Structure of the considered permanent magnet machine

Figures 2 , 3 ,

 23 Figures 2, 3, 4 and 5 summarize for the six proposed formulations and the three objective functions the results obtained in terms of percentage of convergence to a local minimum (Figures2 and 3) and of percentage of times when the best local minimum is found (Figures4 and 5), using totally or partially randomized starting points.From Figure2, we can see that Formulation F4 (which is the most compact one) provides the best results for V m and V u , and Formulation F3 provides the best results for P j . Comparing the six formulations on each of the considered objective functions, it appears that close results were obtained for formulations F2, F3, F4 and F5, while formulations F1 and F6 provide a percentage which is about 3% lower. Contrarily on Figure3(when some components of the starting points are computed from others), F1 and F6 appear clearly to be the most efficient formulations for our Multistart method for all the three objective functions, providing a percentage about 5% higher than all other formulations.Concerning the percentage of best found local minima, from Figures4 and 5we can see that similar results are obtained for all the formulations and all the objective functions.

Figure 2 :F6Figure 3 :F6Figure 4 :F6Figure 5 :

 2345 Figure 2: Percentages of found local minima: totally randomized case

  10) = 0 is added; (iii) we first solve case (i) and following the optimal value of p (denoted p * ) we add the constraint (p -([p * ] -1)) × (p -[p * ]) × (p -([p * ] + 1)) = 0, where [p * ] indicates that the closest integer number of p * is taken into account.

Figure 6 :

 6 Figure 6: Percentages of found local minima

Figure 7 :

 7 Figure 7: Percentages of best found local minima

Table 1 :

 1 Numerical results by minimizing V m in a Multistart method: random case

		F1	F2	F3	F4	F5	F6
	% local min	40.4%	44.2%	43.7%	44.6%	44.0%	41.6%
	% best min	3.7%	3.3%	2.8%	4.6%	4.8%	3.4%
	best min value	7.3589e-5 7.3586e-5 7.3586e-5 7.3568e-5 7.3597e-5 7.3568e-5
	worst local min	1.1443e-3 1.1215e-3 1.1220e-3 1.0526e-3 1.0976e-3 1.1065e-3
	best CPU-time	0.02s	0.00s	0.00s	0.00s	0.00s	0.00s
	avrg CPU-time	0.50s	0.24s	0.27s	0.18s	0.17s	0.42s
	avrg-scf CPU-time	0.12s	0.14s	0.13s	0.14s	0.14s	0.13s
	worst CPU-time	20.88s	27.61s	33.25s	23.02s	30.86s	23.94s
	x0 = mid	-	8.4162e-4 8.4163e-4	-	-	-

Table 2 :

 2 Numerical results by minimizing V m in a Multistart method: par-

	tially random case
		F1	F2	F3	F4	F5	F6
	% local min	49.7%	45.2%	44.7%	44.6%	45.5%	50.2%
	% best min	4.8%	2.9%	3.1%	4.6%	5.0%	5.5%
	best min value	7.3586e-5 7.3586e-5 7.3598e-5 7.3568e-5 7.3606e-5 7.3586e-5
	worst local min	1.1512e-3 1.0382e-3 1.0327e-3 1.0526e-3 1.0077e-3 1.0261e-3
	best CPU-time	0.02s		0.02s	0.02s	0.02s	0.02s
	avrg CPU-time	0.75s	0.17s	0.23s	0.18s	0.16s	0.61s
	avrg-scf CPU-time	0.12s	0.13s	0.13s	0.14s	0.14s	0.13s
	worst CPU-time	20.00s	26.66s	30.91s	23.09s	30.84s	12.98s
	x0 = mid	8.0291e-4 7.7114e-4 7.7112e-4	-	-	8.7127e-4

Table 3 :

 3 Comparison between local minima: completely random case

	Best local min	F1	F2	F3	F4	F5	F6
	D * =	0.159155	0.159155	0.159155	0.159155	0.159155	0.159155
	λ * =	2.5	2.5	2.5	2.5	2.5	2.5
	l * a =	0.003	0.003	0.003	0.003	0.003	0.003
	E * =	0.002549	0.003419	0.003419	0.002897	0.003542	0.002997
	C * =	0.005100	0.004556	0.004556	0.004853	0.004496	0.004790
	β * =	0.8	0.8	0.8	0.8	0.8	0.8
	B * e =	0.382466	0.341724	0.341724	0.363953	0.337209	0.359256
	J * cu =	7485729.14 6464050.67 6464050.68 7021833.60 6350604.11 6904057.59
	K * f =	0.149734	0.182640	0.182640	0.163657	0.186849	0.167446
	e * =	0.001422	0.001426	0.001426	0.001444	0.001414	0.001445
	V * m	7.3589e-5	7.3586e-5	7.3586e-5	7.3568e-5	7.3597e-5	7.3568e-5

Table 4 :

 4 Comparison between local minima: partially random case

	Best local min	F1	F2	F3	F4	F5	F6
	D * =	0.159155	0.159155	0.159155	0.159155	0.159155	0.159155
	λ * =	2.5	2.5	2.5	2.5	2.5	2.5
	l * a =	0.003	0.003	0.003	0.003	0.003	0.003
	E * =	0.003419	0.003419	0.003549	0.002897	0.003633	0.003419
	C * =	0.004556	0.004556	0.004493	0.004853	0.004454	0.004556
	β * =	0.8	0.8	0.8	0.8	0.8	0.8
	B * e =	0.341724	0.341724	0.336951	0.363953	0.334045	0.341724
	J * cu =	6464050.65 6464050.67 6344113.56 7021833.60 6271081.34 6464050.68
	K * f =	0.182640	0.182640	0.187094	0.163657	0.189877	0.182640
	e * =	0.001426	0.001426	0.001413	0.001444	0.001404	0.001426
	V * m	7.3586e-5	7.3586e-5	7.3598e-5	7.3568e-5	7.3606e-5	7.3586e-5

Table 5 :

 5 Numerical results by minimizing V u in a Multistart method

		F1(rnd) F1(cmp)	F4
	% local min	41.0%	49.8%	45.1%
	% best min	2.4%	2.2%	2.4%
	best min value	5.3632e-4 5.3631e-4 5.3635e-4
	worst local min	1.7000e-3 1.6494e-3 1.7543e-3
	best CPU-time	0.02s	0.00s	0.02s
	avrg CPU-time	0.44s	0.73s	0.13s
	avrg-scf CPU-time	0.13s	0.12s	0.20s
	worst CPU-time	19.00s	12.03s	24.09s
	x0 = mid	-	1.5392e-3	-

Table 6 :

 6 Numerical results by minimizing P j in a Multistart method

		F1(rnd)	F1(cmp)	F4
	% local min	44.5%	50.2%	46.1%
	% best min	19.7%	19.6%	21.6%
	best min value	5.7810e+1 5.7810e+1 5.7810e+1
	worst local min	8.5743e+1 7.0680e+1 6.9911e+1
	best CPU-time	0.02s	0.00s	0.02s
	avrg CPU-time	0.28s	0.55s	0.09s
	avrg-scf CPU-time	0.10s	0.09s	0.12s
	worst CPU-time	12.41s	20.45s	12.42s
	x0 = mid	-	58.057	-

Table 7 :

 7 Comparison between local minima using V u

	Best local min	F1(rnd)	F1(cmp)	F4
	D * =	0.159155	0.159155	0.159155
	λ * =	2.5	2.5	2.5
	l * a =	0.003	0.003	0.003
	E * =	0.003373	0.003419	0.003419
	C * =	0.004580	0.004556	0.004556
	β * =	0.8	0.8	0.8
	B * e =	0.343463	0.341724	0.341724
	J * cu =	6507727.58 6464050.67 6464050.67
	K * f =	0.181054	0.182640	0.182640
	e * =	0.001429	0.001426	0.001426
	V * u	5.3632e-4	5.3631e-4	5.3635e-4

Table 8 :

 8 Comparison between local minima using P j

	Best local min F1(rnd) F1(cmp)	F4
	D * =	0.159155 0.159155 0.159155
	λ * =	2.5	2.5	2.5
	l * a =	0.047767 0.042697 0.040992
	E * =	0.001429 0.001429 0.001429
	C * =	0.007328 0.007250 0.007870
	β * =	0.845289 0.803578 0.870069
	B * e =	0.520176 0.541294 0.542699
	J * cu =	10000000 10000000 10000000
	K * f =	0.181863 0.193636 0.227066
	e * =	0.003137 0.003685 0.004109
	P * j	57.8101	57.8101	57.8101
	3.3 Comparison between the six formulations and the
	three objective functions

The second and third author were partially supported by grant ANR 07-JCJC-0151 "ARS".


Sonia Cafieri
email: sonia.cafieri@enac.fr

Pierre Hansen
email: pierre.hansen@gerad.ca

Leo Liberti
email: liberti@lix.polytechnique.fr

Hierarchical clustering for the identification of communities in networks

Keywords: clustering, network, graph, community, modularity, hierarchical heuristic

Introduction

The analysis of networks and in particular the identification of communities, or clusters, is a topic of active research and attracts an increasing attention in the operations research as well as the physics communities. Complex systems arising in a variety of fields can be represented as networks, or graphs, where the set of vertices is given by the entities under study and the edges represent relations holding for pairs of vertices. A typical example is given by social networks, modeling interactions among people. Other real-life applications include communicatons networks, such as the World Wide Web, and transportation networks, representing movements of people or goods.

Given a clustering criterion, the problem of community detection in networks can be formulated as an optimization problem. The current mainstream approach is based on the criterion proposed by Newman and Girvan in [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. These authors introduced the concept of modularity for a partition of a network, defined as the sum for all communities of the difference between the fraction of edges they contain and the expected fraction of edges they would contain if all edges were drawn at random, keeping the same degree distribution. Maximizing modularity gives an optimal partition with its optimal number of clusters.

Community detection based on modularity maximization is currently done with hierarchical as well as with partitioning heuristics, hybrids and, in a few papers, exact algorithms. See [START_REF] Fortunato | Community detection in graphs[END_REF] for an in-depth survey and [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF] for recently proposed exact algorithms. Partitioning schemes aim at finding a single partition or possibly several partitions into given numbers of clusters. They are based on a variety of approaches, including genetic search and simulated annealing. Hierarchical heuristics aim at finding a set of nested partitions. They are in principle devised for finding a hierarchy of partitions implicit in the given network when it corresponds to some situation where hierarchy is observed or postulated. This is often the case, for instance, in social organization and evolutionary processes. Hierarchical heuristics can be further divided into agglomerative and divisive ones. Given a graph G = (V, E) with |V | = n, agglomerative heuristics proceed from an initial partition with n communities each containing a single entity and iteratively merge the pair of entities for which merging increases most the objective function (e.g., modularity). Divisive heuristics proceed from an initial partition containing all entities and iteratively divide a community into two in such a way that the increase in the objective function value is the largest possible, or the decrease in the objective value is the smallest possible. Mergings or bipartitions can be ended once they do not improve the objective function value anymore. Results can be presented on a dendrogram, which displays mergings or divisions of communities.

In this work, we consider the case of hierarchical networks and propose a divisive heuristic which is locally optimal, in the sense that each of the successive bipartitions is done in a provably optimal way. The bipartition subproblem is expressed as a quadratic mixed-integer program with a convex relaxation. To that effect, first we write the modularity Q as a function, for each community, of its number of inner edges and of the sum of degrees of its vertices :

Q = ∑ s m s m - d s 2m 2 , (1)
where m = |E| and m s and d s denote respectively the number of edges and the sum of degrees of the vertices of community s. Here s ∈ {1, 2} since we aim to find a bipartition. We express the sum of degrees d 2 of vertices belonging to the second community as a function of the sum of degrees d 1 of vertices belonging to the first one :

d 2 = d t -d 1
, where d t is the sum of degrees in the community to be bipartitioned. Hence, Q can be rewritten as :

Q = m 1 + m 2 m - d 2 1 4m 2 - d 2 2 4m 2 = m 1 + m 2 m - d 2 1 2m 2 - d 2 t 4m 2 + d t d 1 2m 2 . (2)
We introduce binary variables X r1 , X r2 and Y i1 to identify to which community each vertex and each edge belongs. For r = 1, 2, . . . m and s = 1, 2, X rs is equal to 1 if edge r belongs to community s and 0 otherwise, and for i = 1, 2, . . . n Y i1 is equal to 1 if vertex i belongs to the first of the two communities of the bipartition. We then impose constraints to ensure consistency, i.e. that any edge r = {v i , v j } with end vertices indiced by i and j can only belong to community s if both of its end vertices belong also to that community, and express m s and d 1 in terms of the considered variables : m s = ∑ r X rs ,

d 1 = ∑ i∈V 1 k i Y i1
. We obtain a quadratic convex mixed-integer program which can be solved by CPLEX.

We use this bipartitioning method for the splitting step in a hierarchical divisive scheme. Hence, our divisive heuristic is based on bipartitions which are done at each step in an optimal way. We present a comparison with the spectral-based hierarchical divisive heuristic of Newman [START_REF] Newman | Modularity and Community Structure in Networks[END_REF] and with the hierarchical agglomerative heuristic of Clauset et al. [START_REF] Clauset | Finding community structure in very large networks[END_REF] and we show that the proposed locallyoptimal divisive heuristic gives better results.


Daniel Aloise
email: daniel.aloise@gerad.ca

Sonia Cafieri
email: sonia.cafieri@enac.fr

Gilles Caporossi
email: gilles.caporossi@gerad.ca

Pierre Hansen
email: pierre.hansen@gerad.ca

Leo Liberti
email: liberti@lix.polytechnique.fr

Sylvain Perron
email: sylvain.perron@gerad.ca

Algorithms for network modularity maximization

Keywords: clustering, network, graph, community, modularity, row and column generation

Introduction

Networks are often used to represent complex systems arising in a variety of fields. Social networks model interactions among people. Telecommunication networks model communications betwen them, such as in the World Wide Web. Transportation networks model movements of goods and passengers. Biological networks model interactions between organisms, such as in food networks. A network (or graph) G = (V, E) is composed of a set of vertices, representing the entities of the system under study, and a set of edges joining pairs of vertices and representing a relation holding for such pairs. Identifying communities, or clusters, in complex networks is a topic of particular interest and is currently a very active research domain. A clustering criterion is chosen, in terms of presence or absence of edges, as well as a structure for the clusters. The corresponding mathematical programming problem is then solved exactly or, more often, heuristically. Clustering algorithms and heuristics can be divided into partitioning algorithms, which aim at finding the best partition into a (often but not always) given number of clusters, and hierarchical algorithms, which lead to a set of nested partitions and can be further divided into agglomerative and divisive ones. Algorithms for community identification need a precise definition of community. A very successful definition is that proposed by Newman [START_REF] Newman | Finding and evaluating community structure in networks[END_REF], which leads to the introduction of the concept of modularity for a partition of a network. It is defined as the sum for all communities of the difference between the fraction of edges they contain and the expected fraction of edges if they are placed at random, keeping the same degree distribution. A few algorithms and many heuristics have been developed to maximize modularity. Heuristics rely upon agglomerative and divisive hierarchical clustering, the latter including a spectral approach, and partitioning schemes based on simulated annealing, genetic search, and a variety of other approaches. These heuristics can solve approximately very large instances with up to hundred or thousand entities or even more. However, they do not have either an a priori performance guarantee (e.g. in the form of finding always a solution with a value which is at least a given percentage of the optimal one), nor an a posteriori performance guarantee (e.g. that the obtained solution is at least a computable percentage of the optimal one). By contrast, algorithms provide an optimal solution together with the proof of its optimality. Very few papers propose exact algorithms for modularity maximization. Moreover, they can only solve small instances (with about a hundred entities) in reasonable time. We believe, by analogy with other problems such as the traveling salesman problem, this performance can, and probably will, be much increased. For that purpose, the present paper reviews and compares the algorithms yet proposed. It also presents a new stabilised column generation approach which raised the size of instances solved exactly from 115 to 512 entities.

Row and column generation algorithms

As shown by Brandes et al. [START_REF] Brandes | On modularity clustering[END_REF], modularity maximization in NP-complete, and can be expressed, when there are n entities, as a clique partitioning problem with O(n 2 ) 0-1 variables expressing that pairs of entities belong to the same cluster or not, and O(n 3 ) constraints expressing that if entities i and j belong to the same cluster and entities j and k belong to the same cluster, then entities i and k must also belong to the same cluster. These constraints are numerous and can be added by batches of unsatisfied ones, as already noted by Grötschel and Wakabayashi [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF]. Such row-generation approach can be easily implemented with the "lazy constraints" feature of CPLEX.

Alternately, one can solve this clique partitioning problem by column generation. The auxiliary problem is then an unconstrained quadratic 0-1 program with a 100% dense matrix. It can be solved by a Variable Neighborood Search (VNS) heuristic as long as a column with positive reduced cost can be found and with an exact branch and bound algorithm when this is no more the case.

A different approach was recently proposed by Xu, Tsoka and Papageorgiou [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF]. These authors introduced indicator variables for entities and for edges of G = (V, E) to belong to the first, second,. . .cluster. When G is sparse and the number of clusters is small, the model has less variables and constraints than the clique partitioning one.

Again, a column generation approach can be applied, replacing the difficult problem of finding simultaneously all clusters by a sequence of problems of finding one improved cluster at a time. Column generation algorithms are known to have slow convergence, therefore stabilization is applied [START_REF] Du Merle | Stabilized column generation[END_REF] or, more precisely, focussing. For that purpose, an initial heuristic solution is found by a VNS or other heuristic, from which a plausible interval for the optimal values of the dual variables of the master problem are deduced by sensitivity analysis. Then departures from these intervals are penalized, in a decreasing way.

The four exact algorithms described are compared on a series of problems from the literature. Each of them except that one of Xu et al. is best in some case. The new column generation algorithm outperforms the others for the larger instances.


Alberto Costa
email: costa@lix.polytechnique.fr

Sonia Cafieri
email: sonia.cafieri@enac.fr

Pierre Hansen
email: pierre.hansen@gerad.ca

Reformulation of a locally optimal heuristic for modularity maximization

Keywords: binary decomposition, clustering, modularity maximization, reformulation

Introduction

A network, or graph, G = (V, E) consists of a set of vertices V = {1, . . . , n} and a set of edges E = {1, . . . , m} connecting vertices. One of the most studied problems in the field of complex systems is to find communities, or clusters, in networks. A community consists of a subset S of the vertices of V where inner edges connecting pairs of vertices of S are more dense than cut edges connecting vertices of S to vertices of V \S. Many criteria have been proposed to evaluate partitions of V into communities. The best known of them appears to be the modularity, defined as follows by Newman and Girvan [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]:

Q = c Q c = c m c m - D c 2 4m 2 , ( 1
)
where Q c is the modularity of the cluster c, m c is the number of edges with both end vertices within the cluster c, D c is the sum of the degrees of the vertices in the cluster c, and m is the number of edges of the whole network. The modularity is the difference between the fraction of edges within communities and the expected fraction of such edges in a random graph having the same distribution of degrees than the graph under study. In order to find a good partition into communities for a given network, according to Newman and Girvan one should maximize its modularity. This is a strongly NP-hard problem [START_REF] Brandes | On Modularity Clustering[END_REF].

A few exact algorithms [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF][START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF][START_REF] Xu | Finding community structures in complex networks using mixed integer optimisation[END_REF]] and many heuristics have been proposed for network modularity maximization. They consist in divisive and agglomerative hierarchical clustering approaches [START_REF] Clauset | Finding and evaluating community structure in very large networks[END_REF][START_REF] Newman | Modularity and community structure in networks[END_REF], as well as exact or approximate partitioning ones. In this paper, we focus on a recent locally optimal heuristic based on a hierarchical divisive approach [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF]. We propose several ways to reformulate the model of [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF] in order to accelerate the resolution by reducing efficiently the number of variables and constraints. Computation results are reported for a series of real-world problems from the literature in which the different reformulations are compared. It appears that computing times are very substantially reduced.

Initial model

The model used in the framework of the hierarchical divisive heuristic proposed in [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF] to split a cluster (V c , E c ) into two clusters maximizing the modularity, and based on the one proposed in [START_REF] Xu | Finding community structures in complex networks using mixed integer optimisation[END_REF], is the following:

max 1 m m 1 + m 2 - 1 2m D 1 2 + D c 2 2 -D 1 D c (2) s.t. X i,j,1 ≤ Y i ∀(v i , v j ) ∈ E c (3) X i,j,1 ≤ Y j ∀(v i , v j ) ∈ E c (4) X i,j,2 ≤ 1 -Y i ∀(v i , v j ) ∈ E c (5) X i,j,2 ≤ 1 -Y j ∀(v i , v j ) ∈ E c ( 6
)
m s = (vi,vj )∈Ec X i,j,s ∀s ∈ {1, 2} (7)
D 1 = vi∈Vc k i Y i,1 (8)
Y i ∈ {0, 1} ∀v i ∈ V c ( 9
)
X i,j,s ≥ 0 ∀(v i , v j ) ∈ E c , ∀s ∈ {1, 2}, (10)
where the variable X i,j,s is equal to 1 if the edge (v i , v j ) is inside the community s (i.e., both vertices v i and v j are inside the community s) and 0 otherwise, Y i is equal to 1 if the vertex v i is inside the community 1, and 0 otherwise, and k i is the degree of the vertex v i ; note that D c is a parameter, and it is known before solving the problem.

Reformulations

Power of two reformulation

The heuristic proposed in [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF] works by recursively splitting a cluster into two clusters in an optimal way (in the sense that the computed bipartition corresponds to the best possible modularity). The model is a quadratic integer programming one, with a convex relaxation. The only non-linear term is D 1 2 . The usual Branch-and-Bound approach implemented in CPLEX [START_REF]ILOG CPLEX 12.2 User's Manual[END_REF] is to relax the integrality constraints, solve the continuous quadratic program obtained and then branch. Alternately, one may linearize D 1 2 by replacing it with its expansion in power of two, as proposed for mixed-integer quadratic programming in [START_REF] Billionnet | Extending the QCR method to general mixed-integer programs[END_REF]:

D 1 = t i=0 2 i a i , a i ∈ {0, 1}. (11)
Therefore, the term D 1 2 in (2) can be written as:

D 1 2 = t l=0 2 l a l • t h=0 2 h a h = t l=0 t h=0 2 l+h a l a h = t l=0 t h=0 2 l+h R lh = t l=0 2 2l a l + t l=0 h<l 2 l+h+1 R lh , ( 12
) where R lh is the linearization variable for a l a h ; hence, we have to adjoin the following constraints to our model: R lh ≥ a l + a h -1, ∀l ∈ {0, . . . , t}, ∀h ∈ {0, . . . , l -1} R lh ≥ 0, ∀l ∈ {0, . . . , t}, ∀h ∈ {0, . . . , l -1}.

To estimate t, recall that the maximum value which can be assumed by D 1 is the sum of the degrees of all the vertices in the current cluster, that is D c . Moreover, from (11) the maximum possible value for D 1 is 2 t+1 -1. Hence, t can be computed as:

2 t+1 -1 ≥ D c ⇒ t = ⌈log 2 (D c + 1) -1⌉. ( 13
)

Change of variables

The model of [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF] uses variables assigning edges or vertices to a specific community. When bipartitioning, as there are only two communities to be determined at each iteration, one can use other variables S i,j , associated with the fact that the two end vertices v i and v j of an edge belong to the same cluster or not (i.e., S i,j = 1 if Y i = Y j , and 0 otherwise). This leads to the following reformulation:

max 1 m   (vi,vj )∈Ec (2S i,j -Y i -Y j ) + |E c | - 1 2m D 1 2 + D c 2 2 -D 1 D c   (14) s.t. S i,j ≤ Y i ∀(v i , v j ) ∈ E c (15) S i,j ≤ Y j ∀(v i , v j ) ∈ E c ( 16
)
D 1 = vi∈Vc k i Y i ( 17
)
Y i ∈ {0, 1} ∀v i ∈ V c . ( 18
)

Symmetry breaking

To avoid considering twice equivalent solutions, one fixes a vertex to belong to the first (or second) community. It appears that the vertex with largest degree is a good choice.

Compact model

Applying all the reformulations presented in the previous sections leads to the following compact model:

max 1 m   (vi,vj )∈Ec (2S i,j -Y i -Y j ) + |E c | - 1 2m t l=0 2 2l a l + t l=0 h<l 2 l+h+1 R l,h + D c 2 2 -D 1 D c   (19) s.t. S i,j ≤ Y i ∀(v i , v j ) ∈ E c (20) S i,j ≤ Y j ∀(v i , v j ) ∈ E c (21) R l,h ≥ a l + a h -1 ∀l ≤ t, ∀h < l (22) R l,h ≥ 0 ∀l ≤ t, ∀h < l (23) D 1 = t l=0
2 l a l (24)

D 1 = vi∈Vc k i Y v ( 25
)
Y g = 0, g = arg max{k i , ∀v i ∈ V c } (26) Y i ∈ {0, 1} ∀v i ∈ V c ( 27
)
a l ∈ {0, 1} ∀l ≤ t. ( 28

Results

Table 1 presents the comparison of computing times for the initial model and the final one. Results have been obtained on a 2.4GHz Intel Xeon CPU of a computer with 24 GB RAM running Linux and CPLEX 12.2 [START_REF]ILOG CPLEX 12.2 User's Manual[END_REF]. M denotes the number of clusters, and Q the modularity; computing times are in seconds. Note that slight discrepancies may arise in the values of M and Q; they are due to the fact that optimal bipartitions are not necessarily unique. It appears that the computing time is reduced by a factor of 2 to over 265.

Network

)

  This model has |V c |+t+1 binary variables, |E c |+ t 2 +t 2 +1 continuous variables and 2|E c |+t 2 +t+3 constraints, while the initial model has |V c | binary variables, 2|E c | + 3 continuous variables and 6|E c | + 3 constraints.

  TAB. 1: Results obtained with the hierarchical divisive heuristic using respectively the original formulation and the compact reformulation.

					Initial model		Compact model
		n	m	M	Q	time	M	Q	time
	Karate	34	78	4 0.4188	0.32	4 0.4188	0.16
	Dolphins	62	159	4 0.5265	1.45	4 0.5265	0.65
	Les misérables	77	254	8 0.5468	4.47	8 0.5468	0.67
	A00 main	83	135	7 0.5281	0.71	7 0.5281	0.37
	P53 protein	104	226	7 0.5284	16.82	7 0.5284	1.55
	Political books	105	441	4 0.5263	16.74	5 0.5244	2.66
	Football	115	613 10 0.6009	238.47	10 0.6009 82.21
	A01 main	249	635 15 0.6288	563.41	15 0.6288 38.12
	USAir97	332 2126 8 0.3596 113545.00 8 0.3596 428.40
	Netscience main 379	914 20 0.8470	11.83	20 0.8470	5.24
	S838	512	819 15 0.8166	24.48	15 0.8166	6.40
	Power	4941 6594 40 0.9394	3952.72	41 0.9396 567.07


Implementing an interval computation library for OCaml on x86/amd64 architectures

Jean-Marc Alliot1 and Jean-Baptiste Gotteland 1,2 and Charlie Vanaret 1,2 and Nicolas Durand 1,2 and David Gianazza 1,2 Abstract. In this paper we present two implementation of interval arithmetics for Ocaml on x86/amd64 architectures. The first one is simply a binding to the classical MPFI/MPFR library. It provides access to multi-precision floating point arithmetic and multi-precision floating point interval arithmetic. The second implementation has been natively written in assembly language for low-level functions and in ocaml for higher-level functions ans is as fast as classical C or C++ implementations of interval arithmetic.

Fundamentals of interval arithmetic

Interval arithmetic has been used in computer science and numerical computations for years [START_REF] Moore | Interval Analysis[END_REF]. Its main goal was to create computing environments where the exact value of a computed result lies with certainty within an interval, which might be paramount for some critical applications. Floating point units (FPU) only work with a fixed size for the mantissa of the operands, and numerical errors are unavoidable. For example, The IEEE-754 standard requires that every FPU is able to make any elementary computation in four rounding modes: toward +∞ (upper rounding), toward 0, toward -∞ (lower rounding) and toward the nearest representable number (nearest rounding). It is thus quite easy to implement the most basic functions of interval arithmetic with proper roundings of both endpoints. The addition of two intervals is for example defined by:

[a, b] + [c, d] = [a + low c, b +up d].
Multiplication and division are slightly more complex tasks that require to take into account the position of the intervals relatively to 0. For example, if a < 0 < b, we have

1/[a, b] = [-∞, 1/a] ∪ [1/b, +∞].
As the second member of the equation has to be reduced to a single interval, we simply have:

1/[a, b] = [-∞, +∞]
There is thus a loss of information, which is unavoidable. The usual realvalued functions (cos, sin, arcsin, log, etc.) can also be extended to interval arithmetic. The extension is trivial for monotonic functions (such as exp(x)), as the image of interval [a, b] by exp is simply [exp low (a), exp up (b)]. Computing interval extensions of periodic functions such as sin(x) or cos(x) is a much more complex task. When elementary functions are defined, more complex functions can be computed by composition.

Existing implementations

There are numerous interval arithmetic implementations with various bindings to different languages such as Profil/BIAS (a C++-class library developed in 1993 at the Hamburg University of Technology), a template class for interval arithmetic in the Boost C++ libraries, Gaol (a C++ interval arithmetic library that implements operators for interval constraint programming), MPFI [START_REF] Revol | Motivations for an arbitrary precision interval arithmetic and the MPFI library[END_REF] (a multi-precision interval arithmetic library for C or C++), and the SUN interval arithmetic implementation for Fortran 95 or C++ [START_REF]C++ Interval Arithmetic programming manual[END_REF]. There is a proposal to have interval arithmetic integrated to the standard C++ language [START_REF] Brönnimann | A proposal to add interval arithmetic to the c++ standard library[END_REF] and an IEEE interval standard is currently under development.

Bindings to the MPFR/MPFI libraries

We chose to develop bindings to the MPFR/MPFI libraries since these libraries provide a unique feature: the possibility to work with an arbitrary precision. This feature is extremely valuable when precision (and not time) is paramount. The MPFR library is a multiprecision arithmetic library which implements a very large number of functions on arbitrary-precision floating-point numbers. We only implemented the bindings to MPFR functions which are needed when using the MPFI bindings (more complete MPFR bindings are available in the APRON library). The MPFI library uses the MPFR library to provide arbitrary-precision floating-point interval arithmetic. OCaml bindings exist for almost all MPFI functions with almost no syntactic sugar. Thus, an interval is an opaque structure that must be allocated with the init function prior to use. After allocation, it can be used and modified in place. The semantic of the functions is thus not functional whatsoever. Marshalling of these objects is not possible yet.

The MPFR/MPFI library is widely used, but is quite slow due to its versatility. We thus decided to implement a faster library limited to double-precision floating-point interval arithmetic. There were two possibilities: developing bindings to an existing library, or implementing the library by ourselves. From an overview of the existing implementations, we concluded that developing bindings to these libraries (usually in C++) would be tedious and would require the use of an extra layer that would slow us down. We chose to directly implement the library in assembly language for the lower-level functions, and in OCaml for the higher-level functions. The latter implements the "logical" part of the function, while the former implements the elementary function computations for both rounding modes. The target architecture was the Intel processor family (x86/amd64), on the three main operating systems (Windows, Linux and Mac OS X). We decided to use a functional semantics for all functions: implemented operators such as $+, $-, $*, $/ extend standard operators to interval arithmetic. This allows us to write let a = b $+ c with a, b and c being interval objects. Implementations choices regarding inclusion functions and throwing exceptions are based on the following properties: Let f be a function and F its extension to interval arithmetic, then

• ∀x ∈ [a, b], if f (x) is defined then f (x) ∈ F ([a, b]) • if {f (x) | x ∈ [a, b]} is empty then F ([a, b]) raises an exception
All the elementary operations with both rounding modes are available via the Fpu module of our library, while the interval function extensions are available in the Interval module.

We had to deal with some unexpected problems. The x87 is supposed to return the nearest value, the upper and lower bounds for each elementary operation, but this is not always the case: some functions such as cos, sin or tan are not properly implemented everywhere. For example, we computed cos(a), with a = atan2_low 1. 0. and with the following cosine implementations:

1. the MPFI library (with 128-bit precision), 2. the x87 in round-toward -∞ mode, 3. the x87 in nearest mode (default value for the C and OCaml libraries on 32-bit Linux), 4. the x87 in round-toward +∞ mode, 5. the SSE2 implementation (default value for the C and OCaml libraries on 64-bit Linux):

We got the following results: cos x87low (a) < cosx87(a) = cos x87high (a) < cosMP F I (a) < cosSSE2(a), so the upper bound (4) computed by the x87 is clearly incorrect, as it is lower than the correct value computed by the MPFI library. The value computed by the SSE2 ( 5) is much more precise than the one computed by the x87. However, it is unfortunately impossible to get upper and lower bound values with the SSE2 implementation, and we have no other choice but to use the x87 for computing these (sometimes incorrect) bounds. The problem here is that the value computed by the standard C-lib (or OCaml) cos(x) function does not always lie within the lower bound/upper bound interval returned by the x87 functions. This can be very prejudicial when executing branch and bound algorithms where the mid-value is expected to lie within the lower/upper interval. We solved this issue by rewriting quite efficiently the trigonometric functions in assembly language. With our new implementation, the lower and upper bounds are properly set and they are always lower (resp. higher) than the values computed by the standard cos(x) functions on 32 and 64-bit architectures.

Values returned by the standard (C-lib or OCaml) cos(a), sin(a) or tan(a) functions remain different on 32 and 64-bit architectures. In order to obtain the same behavior on both architectures, the fcos, fsin or ftan functions from module Fpu can be used. They always return the same values on all architectures, and they can transparently replace standard functions by using the Fpu_rename or Fpu_rename_all modules. From the above figure and table, we can compare our native implementation with some of the fastest implementations of interval arithmetic. We can observe that on most functions, our library is on par with Filib. Moreover, the native implementation suffers from its functional semantics, as it creates at each of the 10 6 operation a new interval object, while the classical libraries (including MPFI) create interval objects once and for all before the computation and modify them in place. So the observed performance is indeed excellent.

Performance issues

Figure 1

 1 Figure 1 from [3] presents some comparison results of classical interval libraries, including MPFI. The scale is logarithmic.

Figure 1 .

 1 Figure 1. Comparison of standard interval libraries

Table 1

 1 is a comparison of MPFI with our native implementation on 10 6 operations (times are in seconds). The third column is the logarithm of the ratio of the execution times of both implementations (log 10 (M/N )).

	Op	MPFI	Nat	log(M/N)	Op	MPFI	Nat	log(M/N)
	+	0.124	0.076	0.21	log	2.696	0.096	1.44
	-	0.172	0.068	0.40	exp	4.568	0.224	1.31
	*	0.148	0.088	0.23	cos	3.544	0.136	1.42
	/	0.240	0.088	0.43	sin	3.868	0.136	1.45

Table 1 .

 1 Comparison of MPFI bindings and native implementation

Institut de Recherche en Informatique de Toulouse, name.surname@irit.fr

Laboratoire "Mathématiques Appliquées et Informatique" de l'ENAC