Interval-based global optimization in engineering using model reformulation and constraint propagation

Engineering Applications of Artificial Intelligence 25 (2012) 404–417

Contents lists available at SciVerse ScienceDirect

Engineering Applications of Artificial Intelligence journal homepage: www.elsevier.com/locate/engappai

Interval-based global optimization in engineering using model reformulation and constraint propagation Issam Mazhoud a, Khaled Hadj-Hamou a,n, Jean Bigeon a, Ghislain Remy b a b

Grenoble-INP/UJF-Grenoble 1/CNRS, G-SCOP UMR5272, Grenoble F-38031, France LGEP CNRS/SUPELEC, 11 rue Joliot Curie, Plateau de Moulon, 91192 Gif sur Yvette, France

a r t i c l e i n f o

a b s t r a c t

Article history: Received 1 October 2010 Received in revised form 4 May 2011 Accepted 20 October 2011 Available online 21 November 2011

This paper deals with the preliminary design problem when the product is modeled as an analytic model. The analytic models based method aims to use mathematical equations to address both multiphysic and economic characteristics of a product. The proposed approach is to convert the preliminary design problem into a global constrained optimization problem. The objective is to develop powerful optimization methods enough to handle complex analytical models. We propose to adapt an approach to solve this problem based on interval analysis, constraint propagation and model reformulation. In order to understand the optimization algorithm used for engineering design problems, some basic definitions and properties of interval analysis are introduced. Then, the basic optimization algorithms for both unconstrained and constrained problems are introduced and illustrated. The next section introduces the reformulation technique as main accelerating device. An application of the reformulation device and its global optimization algorithm on the optimal design of electrical actuators is presented. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Preliminary design Analytic model Global optimization Interval analysis Constraint propagation Model reformulation

1. Introduction Global optimization is an important issue in the field of product engineering allowing to reduce costs, to improve product performance, reliability and shorten design time. The objective of this paper is to develop powerful optimization algorithms enough to handle complex engineering problems and particularly electromagnetic devices (electric actuators). Many authors worked on the design theory and proposed several descriptions of the design process. The systematic design approach proposed by Pahl et al. (2007) and described in Fig. 1 is one of the most used design theory. The main task of the design engineers is to propose an ‘‘optimal’’ design solution satisfying all the customer needs. The feasibility of a technical solution is related to the constraints satisfaction including technological, economic, environmental and legal aspects. The optimization step is performed according to one or more criteria (economic, performances, etc.). Fig. 1 shows the design steps proposed by Pahl et al. (2007) related to the design of an electric actuator. Firstly, the customer needs are expressed through the specifications, in which some n

Corresponding author. Tel.: þ33 4 76 57 48 11; fax: þ 33 4 76 57 48 95. E-mail addresses: [email protected] (I. Mazhoud), [email protected] (K. Hadj-Hamou), [email protected] (J. Bigeon), [email protected] (G. Remy). 0952-1976/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.engappai.2011.10.010

constraints are imposed. Generally, these constraints are related to the product performances, dimensions and/or the materials. Then, these needs are analyzed, adapted if necessary and validated. As the electric actuator is a well known concept, the mathematical model expressing its physical behavior is known. According to the classification of Gero (2001), the example of the electric actuator design belongs to the routine design class. The next step consists in formulating the final analytical model that takes the specifications’ constraints into account. During the preliminary design step, a first quantification of the design parameters is given using an optimization tool. The most used method to get a first design solution is the Finite Elements Method (FEM). But this complex method is becoming more and more inappropriate during the preliminary design since it is based on the device geometry. Thus, the multi-physic behavior and the economic objectives are poorly managed and the computations are time-consuming. Let us take, as example, the electromagnetic actuator design problem. This kind of product involves simultaneously the electromagnetic, the mechanical and the heat fields. Therefore, in order to describe the physical behavior of the electromagnetic actuator, it is necessary to solve complex partial differential equations: Maxwell equations for the electromagnetic behavior, the heat equations for the heat behavior and so on. Generally, design engineers make assumptions in order to integrate those equations and get consequently analytic equations that explicitly link the product performances to the

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

405

Fig. 1. Simplified steps of the design process (Pahl et al., 2007).

design parameters. Roters gives in Roters (1941), a set of analytic equations based on simplifying assumptions for the electromagnetic devices. This is all the more interesting as the customer specifications generally deal with the geometry, the materials and the performances parameters simultaneously. Unlike the FEM method, the analytic models based method aims to use mathematical equations describing the technical and the economic characteristics of the device. These models respond to the preliminary design complexity and the use of a global optimization approach aims to build an initial design solution. The final definition of the structure is validated using a CAD software and an FEM analysis of the design solution. Finally, the production and operating documents are prepared. This paper deals with the first step of the preliminary design step: computation and dimensioning. It attempts to propose an optimization tool adapted to a wide range of design problems. This tool is particularly suitable for the routine design as the analytic models are generally available. During the preliminary design step also called embodiment design, the concepts selected during the conceptual design step starts to materialize the product and the preliminary layouts are generated. These layouts are then optimized and several criteria could be used to evaluate them in order to select the most promising. Finally, the selected layout will be the object of a detailed design. Consequently, an important part of the product costs is incurred in the preliminary design step. Thus, it is important to use an optimization tool.

The selection of the solving algorithm is extremely important. This will be guided by the problem formulation and the expected results. Numerous optimization techniques could be used at the preliminary design step. These techniques may be classified according to two criteria: the operating approach (stochastic or deterministic) and the quality of the solution sought (local or global). Many authors worked on design optimization using stochastic approaches. Isfahani and Vaez-Sadeh (2007), Ahn et al. (2003) and Ho et al. (2008) applied a Genetic Algorithm to the design of engineering devices: electric motors and squeeze film damper. Pierre et al. (1995) applied a Simulated Annealing algorithm for the optimal design of a computer communication network. He and Wang (2007) proposed a Particle Swarm Optimization for constrained engineering design problems. The advantage of the stochastic algorithms is their wide range applicability and their implementation simplicity. Indeed, the objective function and the constraints are not required to be linear, convex, differentiable and continuous. Nevertheless, these algorithms do not give optimality neither feasibility proofs. In addition, they need some parameters setting which may affect the algorithm convergence. Deterministic algorithms have also been used in engineering design. Gol and Sobhi-Najafbadi (2005) and Lin and Lin (2001) used the Sequential Quadratic Programming local algorithm for the optimal design of, respectively, electromotion devices and rotor systems. The main drawback of the local methods is that

406

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

they do not give global optimality proof. Many global deterministic algorithms exist to solve problems with particular structures: convex optimization for convex functions, simplex algorithm for linear programs, quadratic programming for quadratic functions. The engineering problems we solve belong to a large and complex class of optimization problems with non-linear and non-convex functions. Messine (2004) proposed a global deterministic algorithm to solve general constrained optimization problems called Interval Branch and Bound Algorithm (IBBA). This algorithm provides a global optimality proof but it is difficult to find exact global optima in a reasonable runtime. The purpose of this paper is to adapt the IBBA algorithm to the preliminary design context and enough to handle complex analytic models. An analytic model consists of two main parts: a physico-mathematical model and the user specifications. (1) A physico-mathematical model: It describes the physical behavior of the device. This model is usually made by modeling experts from several engineering fields (mechanics, heat, electrical, etc.), and is more or less independent of the manufacturing context. It is defined as a set of analytic equations related to several physical laws. Consequently, they are tight as well as necessary to satisfy and it is better not to approximate them. A physico-mathematical model has several characteristics: – it involves generally 10 up to 1000 parameters and usually up to 20 degrees of freedom, – the mathematical formulation of the equations can either be explicit, implicit and/or including functional equations (ordinary differential equations, integrals, MinMax, etc.), – it is usually multi-physic, it can incorporate several engineering fields (mechanics, heat, electromagnetic, etc.), – the analytic equations are often non-linear, having different forms according to the field: polynomial, ln, exp, cos, sin, etc., – the mathematical expressions are mostly continuous but may involve some discrete variables so as to make the problem highly combinatorial, – the convexity of the equations is not known and is hardly determinable. In many product design cases, the physico-mathematical expressions have the form pi ¼hi(P), where pi is a product performance parameter and P the vector of the design parameters.

(2) Specifications: They define the user needs. Generally, they are established by the customer or the marketing department. The specifications define the user requests in terms of production, performances and costs including: – one or more objective functions (minimization or maximization), for example cost, device power, device volume, etc., – parameters variation domains representing the bound constraints for the parameters, – user new constraints: fixing some performance and design parameters, adding new constraints related to the manufacturing context, etc., As described in Fig. 2, the model that should be solved is the aggregation of the physico-mathematical model and the specifications. Experts in computer science and applied mathematics combine these two parts to draw up the final design model. In this paper, the approach is to convert the preliminary design model into a global constrained optimization problem. This paper only deals with the reformulation and the optimization steps. We consider that the analytical model related to a specific design application is available. Often, we use to reconsider the specifications. Therefore, it is essential to be able to modify the model quickly and to get solutions within a reasonable time. The final model takes into account the objective function, the constraints and the variables (parameters) and can be formulated as a general global optimization problem: min f ðxÞ

x A X D Rn

( subject to :

g k ðxÞ r 0,

8k A f1,. . .,pg

hk ðxÞ ¼ 0,

8k A fp þ1,. . .,p þ qg

where X is the hypercube defining the bound constraints for the variables, x is a realization of all the variables belonging to X, f is the objective function, G ¼{gk, kA{1,y,p}} is the set of the inequality constraints and H¼{hk, kA{pþ 1,y,p þq}} is the set of the equality ones. An inequality constraint can be reformulated into an equality by adding a slack variable (xi ¼gi(x)) such that (xi r0 ) xiA]  Ny0]). In addition, we assume that the constraints are explicit, meaning that they have the form pi ¼ hi(P). Then, the optimization

Design Engineers

Analytical Model Multi-physic Economic

Modeling Experts

Analysis + Reformulation

Design Engineers

Preliminary design model

Specifications = Constraints + Objective User Customer Fig. 2. The approach to build the preliminary design model.

Deterministic optimization

Definitive layout

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

407

In order to understand the optimization algorithm, some basic definitions and properties of interval analysis are introduced in the second section. In the third one, the basic optimization algorithms for both unconstrained and constrained problems are presented and illustrated. The fourth section introduces the reformulation technique as a main accelerating device. Finally, the last section shows an application of the reformulation device coupled with the global optimization algorithm on the optimal design of electrical actuators.

2. Basics of interval analysis

Fig. 3. The preliminary design model of the electromagnetic actuator.

problem handled in this paper can be formulated as follows: 8 < min n f ðxÞ xAX DR

: pi ¼ hi ðPÞ,

8j A f1,. . .,qg

Fig. 3 shows an example of the analytical model of an electric actuator. This model has 6 analytic constraints and 17 variables. The constraints and the objective function are non-linear and non-convex. In the specifications, some variables are fixed and for the remaining ones, bound constraints are given. This model is detailed in the Section 5. From now on, we consider that the objective function and the constraints are expressively available and continuously differentiable. The models with implicit or differential equations systems are not handled here. This paper propose to adapt the IBBA algorithm and to solve a large class of engineering problems. Within this framework it is possible: – to achieve consistent enclosure for the global optimum and all its optimizers within a precision frame given by the user, – to display reasonable runtime (reduced number of model evaluations) in solving complex engineering preliminary design problems, This algorithm is based on interval analysis, constraint propagation and model reformulation: – Interval-based branch and bound algorithm: It is a global optimization algorithm based on interval computation. It allows to enclose all the global optima or to give a nonexistence proof if there is none (Hansen and Walster, 2004), – Contractors: The contractors are the interval extension of the constraint propagation operators. A contraction as an accelerating device is a procedure that reduces the search domain by discarding non-feasible regions, – Model reformulation: The main contribution of this paper is a reformulation step introduced as an accelerating device and made the problem-solving more powerful. The reformulation takes advantage of the intrinsic characteristics of many preliminary design problems. The main idea is to reduce the complexity of the model so as to reduce the number of model evaluations. This can be reached by reducing the number of the model parameters. As this algorithm is a general optimization method, it could be applied to the innovative or creative design as long as we have the analytical model.

The global optimization algorithm is based on the branch and bound principle. At each iteration, the search domain (a box) is bisected, and the regions that cannot contain the global optimum are discarded. In order to make possible to discard a box, it is essential to evaluate the objective function and the constraints over that box. Rather than evaluating them point by point, we use the interval arithmetic that permits to give lower and upper bounds of a function over a box. The interval arithmetic has been introduced by Moore in 1966 as an approach to overcome rounding errors (Moore, 1966). Since then, many applications and improvements in global optimization have been proposed taking advantages of the increasing processing power of computers (Ratschek and Voller, 1991; Kearfott, 1992; Hansen and Walster, 2004; Messine, 2004). Before explaining the optimization algorithm, let us give some definitions. Definition (Interval). An interval is a connected ([1y2][[3y5] is not an interval), closed ([1y2[is not an interval) subset of R. The set of all intervals of R will be denoted R. Definition (Box). A box is an interval vector defined by the Cartesian product of n intervals: the vector ([1y2],[4y5]) is a bi-dimensional box. An n dimensions box is an element of IRn . 2.1. Elementary operations Consider two intervals X and Y. We define X \ Y ¼ fx=x A X4x A Yg,

X [ Y ¼ fx=x A X3x A Yg,

X\Y ¼ fx=x A X4x= 2Yg,

X  Y ¼ fðx,yÞ=x A X4y A Yg

Definition (Range). Let f be a function from D DRn in Rm , let X A ID . We call image by the function f of X the following subset: range (f,X)¼ {f(x)/xAX}. Instead of defining arithmetic operations on fixed numbers, the interval analysis extends these operations to intervals. All the elementary operators have been redefined for interval. Definition (Elementary operators). Let X and Y be two intervals and 3 A f þ ,,  ,=g: X3Y ¼ fx3y=x A X,y A Yg. Thereby ½a. . .b þ ½c. . .d ¼ ½a þc. . .b þ d ½a. . .b½c. . .d ¼ ½ad. . .bc ½a. . .b½c. . .d ¼ ½minfac,ad,bc,bdg. . .maxfac,ad,bc,bdg ½a. . .b=½c. . .d ¼ ½a. . .b½1=d. . .1=c if 0 2 = ½c. . .d Note that subtraction (respectively, division) is not the inverse operation of addition (respectively, multiplication). Extended interval arithmetics have been proposed to allow division by an interval containing 0 and to allow open intervals and infinite bounds (Hansen, 1980). Other extended arithmetics have been proposed, such as Affine arithmetics (Messine, 2002), symbolic

408

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

interval arithmetics (Jaulin and Chabert, 2010). Elementary functions (exp, ln, power, cos, sin, etc.) have also been extended to intervals. We also define the following notions for intervals and boxes. Let i be an interval [ayb]: – lef tð½a. . .bÞ ¼ a, rightð½a. . .bÞ ¼ b: – mid(i)¼((aþ b)/2) the middle (center) of the interval [ayb]. For a box, the vector middle of the box intervals. For example, for a box X ¼([1y2],[4y5],[4y10]), mid(X)¼ (3/2,9/2,7). – W(i)¼(b a) the width of the interval [ayb]. For a box, it is defined by the width of the widest interval. For example, for a box X ¼[1y2],[4y5],[4y10]), w(X)¼max(1,1,6) ¼6. – The bisection operator allows to explore smaller parts of the search domain. This operator split a box B through a specific direction in order to generate two boxes B1 and B2 such that B¼B1[B2. The bisection of an interval [ayb] 1,1,61,1,6 generates the two equal intervals [ay((aþb)/2)] and [((aþb)/2)yb]. The most common bisection rule is that which defines the splitting direction as the largest edge of the box. The bisection of the box ([1y2],[4y5],[4y10]) by its largest side generates the two boxes ([1y2],[4y5],[4y7]) and ([1y2],[4y5],[7y10]). Nevertheless, there are many other splitting rules such as the bisection following Baumann center (Baumann, 1988), multi-section techniques (Casado et al., 2000). In the next section, we introduce the extension to intervals of more complex functions: the inclusion functions. 2.2. Inclusion functions The constraints addressed in this paper are a combination of elementary operations. In order to handle these constraints, extensions to intervals of the common functions are introduced. The inclusion functions allow extending any usual function to intervals. This property is particularly useful in constraint propagation. Definition (Range). Let f be a function from D D Rn in Rm , let X A ID. We call image of X by f the following subset:   rangeðf ,XÞ ¼ ff ðxÞ=x A Xg : ¼ minf ðxÞ. . .maxf ðxÞ xAX

xAX

Definition (Inclusion function). Let f be a function from D D Rn in Rm . Let F be another function from ID in IRm . The function F is called inclusion function of f in D if and only if: rangeðf ,XÞ DFðXÞ: The closer an inclusion function is to range(f,X), the better it is. As illustrated in Fig. 4, f([X]) is the range of f as it is the exact image of [X] by f. [f] is an inclusion function and [f]n the minimal inclusion function. Knowing that the extensions to intervals of the elementary operations are by definition inclusion functions for the elementary operators, a trivial inclusion function can be built. It is obtained by exchanging the usual elementary operators by their interval extensions. This inclusion function is called the natural inclusion function.

Fig. 4. Inclusion functions [f] and [f]n(minimal) (Jaulin et al., 2001).

Definition (Natural extension). Let f : Rn -R a continuous function. The interval function F : IRn -IR having the same expression as f, in which the classic operators are replaced by their interval extensions, is called the natural interval extension of f. Let us introduce a simple function and compute the image of the interval [.5y2.0]: f :

R-R x/xx2

The natural extension F1 of f has the same expression but using the interval operations:

IR-IR

F1ð½0:5. . .2:0Þ ¼ ½0:5. . .2:0½0:5. . .2:02 X/XX 2 ¼ ½3:5. . .1:75

F1 :

Now, if we use the same function under other forms: F2 :

IR-IR X/Xð1XÞ2

F2ð½0:5. . .2:0Þ ¼ ½2:0. . .1:0

F3 :

IR-IR X/0:25ðX0:5Þ2

F3ð½0:5. . .2:0Þ ¼ ½2:0. . .0:25 ¼ rangeðf ,XÞ

Using the inclusion function F3 we obtain the smallest enclosure of f on the interval [.5y2.0]. So, the quality of an interval evaluation of a function is highly dependent on the form of the inclusion function. The more precise the enclosure of the function is the more effective the optimization algorithm will be, as it deals with smaller search boxes. For the interval evaluation of an expression, each occurrence of the same variable is processed as a different variable. The dependency between the occurrences of a variable in an expression is lost. In F3 the interval variable X occurs only once, allowing no overestimation of F(X). Therefore, other inclusion functions have been introduced: Taylor 1st and 2nd orders (Hansen and Walster, 2004), optimal centered forms (Baumann, 1988), and Kite (Vinko et al., 2004). Definition (Taylor first order extension). Let f : D D Rn -Rm a multivariate continuously differentiable function, X A ID,GðXÞ the gradient enclosure of f and y any point belonging to the box X, for example the middle: range(f,X)Df(y)þG(X)(X y), thus, the function: F t : ID-Rm ,X/f ðyÞ þ GðXÞðXyÞ is an interval extension of f and is called the first order Taylor extension. The Taylor extension is a particular case of centered forms. This extension is particularly interesting since it allows the linearization of any function. This property is very useful in the constraint propagation step. For more details about interval computation, see Moore (1966) and Kearfott, (1996).

3. Interval-based global optimization The global optimization algorithms used in this paper are based on the branch and bound principle; the branching is done by bisecting the search domain (box), and the bound by evaluating the model and discarding boxes that do not contain the optimum. In the case of minimization problems, the Interval Branch and Bound Algorithm (IBBA) maintains an upper bound f~ of the global minimum xn and updates it throughout its execution. If we consider a box X such that f ðXÞ 4 f~ , meaning that the lower bound of the interval f(X) is greater than the current upper bound f~ : lef tðf ðXÞÞ 4 f~ , this box obviously does not contain the global minimum and is consequently discarded. Other tests have been introduced as accelerating devices allowing discarding more boxes. These accelerating tests depends on whether the problem is constrained or not. In addition to the upper bound, the optimization algorithm maintains throughout the iterations two

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

lists of boxes. One list contains the boxes that have to be evaluated (workList). The other one contains the boxes that have not been discarded and have a maximum width E fixed by the user, so that they will not be split (resultList). The algorithm takes an initial box as input and gives as outputs a list of boxes that contain the optima (if there is at least one) and an upper bound of the global optimum which is necessarily a feasible point. To illustrate the interval-based branch and bound principles, the first introduced algorithm deals with the unconstrained case: min n f ðxÞ xAX DR

3.1. Unconstrained IBBA The first algorithms for unconstrained optimization problems have been developed by Hansen (1979, 1980). Three accelerating devices can be used in order to eliminate boxes that cannot contain the global optima: bound test, monotony test and convexity test, as illustrated in Fig. 5. The interval I1 cannot contain the optimum as f is strictly monotonous on it. On interval I2, f is concave, so it can be discarded. I3 is discarded as its objective function value is greater than the actual upper bound. To illustrate the unconstrained problem, we use the Six-hump camel back function: min

x A ½1000...10002 D R2

f 1 ðxÞ ¼ 4x21 2:1x41 þ 13x61 þ x1 x2 4x22 þ 4x42

f1 contains 2 global minima and 4 local minima which are very close. We tested the IBBA algorithm on this function and we compared the results with and without using accelerating tests (monotony and convexity tests). The benchmark review shows that the upper bound of the global minimum obtained is (  1.031628453489). The global optimizers have been enclosed in two symmetrical boxes with a precision on the boxes width of ewidth ¼ 1:0e8: ! ½8:98420184e02. . .8:98420112e02 Opt 1 ¼ ½ þ 7:12656401e01. . . þ 7:12656408e01 ! ½ þ 8:98420112e02. . . þ 8:98420184e02 Opt 2 ¼ ½7:12656408e01. . .7:12656401e01 ( 8x A IRn ,

409

As illustrated in Fig. 6, during the algorithm execution using accelerating tests, 159 boxes have been discarded using the bound test (gray), 204 using the monotony test (white), 10 using the convexity test (black) and 2 boxes containing the global optimizers (dotted boxes). The results summarized in Table 1 show that the accelerating tests are very efficient in the unconstrained cases, even with many local minima. The lower bound obtained using accelerating tests corresponds to the global minimum with a 1.0e  8 precision. The monotony and convexity tests are not appropriate in the case of constrained problems. However, other techniques have been developed in order to take constraints into account.

3.2. Constrained IBBA The accelerating devices that will be developed for the constrained case are derived from the constraint programming, as an adaptation of the constraint propagation techniques. Constraint satisfaction problems have been used and adapted to solve multiple engineering design problems such as in Aldanondo et al. (2008), Yvars (2009) and Qureshi et al. (2010). Definition (Continuous CSP). A continuous CSP (Constraint Satisfaction Problem) is a triplet X,D,C defined by a set of n variables X ¼{x1,x2,y,xn}, a set of the variables interval domains D ¼{D1,D2,y,Dn} such that x1 A D1 ,x2 A D2 ,. . .,xn A Dn and a set of m constraints on the variables C ¼{C1,C2,y,Cm}. For each constraint, rc defines the set of all points of Rn that satisfy the constraint C. The interval extension of constraint propagation is called Contractor, an operator that allows discarding parts of the search space that do not contain feasible solutions (Jaulin et al., 2001). Definition (Contractor). A contractor, also called narrowing operator, for a constraint C is a function.C X : IRn -IRn . that satisfies the two following properties:

C X ðxÞ D xðcontractanceÞ C X ðxÞ \ rc ¼ x \ rc ðcompletenessÞ, where rc is the set of feasible points

Fig. 5. Accelerating devices for unconstrained optimization problems.

410

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

Fig. 6. Paving of the search domain in the global optimum neighborhood.

Table 1 Impact of accelerating devices on algorithm efficiency. f1

Upper bound

Without accelerating devices With accelerating devices

 1.031297687959 1 464 111

Iterations #

 1.031628453489 749

CPU time (s) 17.40 .15

For example, considering the atomic constraint C : x1 ¼ x2 þx3 and the intervals D1,D2 and D3, respectively, for x1 ,x2 and x3, a possible narrowing operator can be computed as follows: 8 > < D1 ¼ D1 \ ðD2 þD3 Þ D2 ¼ D2 \ ðD1 D3 Þ > : D ¼ D \ ðD D Þ 3

3

1

Compared to the unconstrained IBBA, three main adjustments are made to implement a constrained IBBA (also called Branch and Prune):

2

A complete library has been implemented in Java in order to compute narrowing operators for all atomic constraints (plus/ minus, mult/div, exp/ln, sin/arcsin, cos/arccos, sqr/sqrt, abs). For atomic and linear constraints, a contraction function can be formulated easily because each variable can be expressed with the other ones. This task is much harder for a general engineering design constraint, as they are generally non-linear and nonconvex. For more complex constraints, precisely non-linear constraints, two solutions have been proposed: 1. A method based on linearization using centered forms such as Taylor first order extension (Hansen and Walster, 2004). For a non-linear constraint C(x)¼0 defined on a box X, it is possible to express the Taylor first order interval extension: 8y A X,

where Xi the ith component of X,C 0i ðXÞ the ith component of G(X). Therefore, to use this contraction technique, it is necessary to compute all the first order derivatives of the constraints. Hand computing of all the derivatives is not an easy task for real engineering design problems. To handle this issue, techniques and softwares for automatic differentiation are available, such as Matlab and Pro@Design (DPT, 2010). 2. A method based on the calculus tree named Hull consistency (Benhamou et al., 1999). The main idea is to create intermediate constraints in order to transform a complex constraint into a set of elementary constraints (atomic constraints). This technique has been tested and shown good performances (Messine, 2004).

CðXÞ DCðyÞ þ GðXÞðXyÞ

where G(X) is the gradient’s enclosure of C. Then, we can express each variable Xi using the others: ! P CðyÞ nj¼ 1,j a i Gj ðXÞðX j Y j Þ X i ¼ X i \ yi þ C 0i ðXÞ

– each box taken from the workList is systematically pruned (contracted), – during the updating step, a feasibility test is added: the upper bound must be feasible, – the stopping test contains an additional condition ðf~ binf o ebound Þ. This test reflects a requirement on the quality of the solution. Through this test, we are certain that the global minimum is greater than ðf~ ebound Þ. The algorithm has a branch and bound structure. Inside, the interval analysis is used to evaluate enclosure functions over the boxes, and the contractors allow to reduce the search domain (Algorithm 1). Algorithm 1. Constrained global optimization algorithm 1: Initialize X’the initial box in which the global minimum is sought. 2: Initialize f~ ’ þ 1 3: Initialize the list workList’ð þ 1,XÞ : the structure that will contain the boxes that remain to be explored. The elements in the workList have the form: (lower bound of f, Box).

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

4: Initialize the list resultList:empty while (workList not empty) 5: take an element ðvj ,X j Þ from workList 6: bisect Xj into two boxes: X j1 and X j2 for i from 1 to 2 7: contract the box by constraint propagation techniques 8: compute the lower bound of f in X ji : binf if ðb 4 f~ Þ inf

9: discard Xji else 10: C’Center(Xji) if (f ðCÞ o f~ and C is feasible) 11: f~ ’f ðCÞ end if if ½wðX ji Þ o ewidth and f~ binf o ebound  12: resultList’resultList þðbinf ,X ji Þ else 13: workList’workList þðbinf ,X ji Þ end if end if end for end while Complete interval analysis library and optimization algorithms have been implemented in Java version 1.6.0_19. The tests have been accomplished on a 2.6 GHz Intel core i5 with 8Go-RAM. We have realized tests to evaluate the contribution of the contractors on two mathematical models. They are relatively simple examples but they will permit us to evaluate the effectiveness of the contractors and to identify some problems. The test case 2 has two variables allowing us to show the efficiency of the acceleration devices on a 2D graphic Test case 1 min f ðxÞ ¼ x3 þ ðx1 þ x2 þ x3 Þx1 x4 ( x21 þx22 þx23 þ x24 ¼ 40 , X A ½1 . . . 54 subject to : x1 x2 x3 x4 Z 25

411

feasible solutions. Four boxes enclosing the global optimum have been selected (in gray). The global optimum corresponds to the last upper bound update: X n ¼ ð0:84089641,2:99999999Þ. For more details and recent advances about contractors see Benhamou et al. (1999), Chabert and Jaulin (2009) and about constrained interval-based optimization see Ratschek and Rokne (1988), Byrne and Bogle (1996), and Hansen and Walster (2004). Despite the efficiency of the accelerating devices, especially the constraint propagation, updating the upper bound of the global optimum remains difficult. In fact, even using the contractors, IBBA was able to find only 4 feasible solutions for the test case 1 (more than 2 millions iterations) and only 1 feasible solution for the test case 2 (74 iterations). Currently, the upper bound is updated by testing whether the middle of the current box is feasible (Algorithm 1). It is very difficult to find a feasible solution when instantiating the variables independently of the constraint especially if the constraints are very tight, which is particularly true for equality ones. Knowing that the variables are strongly linked to each others in the equality constraints, it would be interesting to use this information in order to reduce the number of variables to be determined. If we note nv the number of variables and nc the number of equality constraints, supposed explicit and independent, the problem can be summarized in the determination of (nv nc) variables, the remaining variables are directly computed using the constraints. These independent variables are called the degrees of freedom. This model simplification can be reached through a reformulation step.

4. Model reformulation Frequently, the constraints of a particular optimization problem may be modeled by several equivalent CSPs. Performance

Test case 2 min f ðxÞ ¼ 12x1 7x2 þx22 subject to : 2x41 þ2x2 ¼ 0, X A ½0 . . . 2  ½0 . . . 3

The numerical results are summarized in Table 2. Without the use of contractors, neither the first nor the second test case have been solved. This is explained by the number of boxes in the worklist that increases exponentially and exceeds the size of the heap memory (512Mo). In addition, no updates have been done as no feasible point has been found. Consequently, no boxes can be discarded. In the other side, the contractors reduced significantly the search space which permits to easily find feasible points and consequently update the upper bound. As illustrated in Fig. 7 concerning test case 2, the upper bound of the global minimum obtained is (f(Xn)¼ 22.09075666) with a precision of ewidth ¼1.0e  5. During the execution of the algorithm using a contractor accelerating device, 116 boxes have been discarded using the bound test (in black), only one upper bound update and an important part of the search area discarded using constraint propagation (white parts). The white parts are not necessarily boxes as they are parts of boxes that do not contain

and behavior of an optimization solver may change drastically even with equivalent CSPs. To choose the best model for a particular problem, it is important to understand the underlying optimization and constraint propagation algorithms. Some modeling techniques are commonly adopted to improve the accuracy and efficiency of the interval-based optimization solvers. A fundamental problem of interval analysis is the dependency. In the interval arithmetic evaluation of an interval expression, each occurrence of the same variable is processed as a different variable. The dependency between the different occurrences of a variable in an expression is lost. Some expressions may be rewritten into equivalent ones that minimize this problem. For polynomial expressions, we factorize as much as possible to obtain functions with mono variable occurrences (Section 2.2). For the example function f ðxÞ ¼ xx2 , instead of using the constraint inclusion function F1(X)¼X X2 we use the optimal inclusion function F3(X)¼.25(X .5)2 ¼range(f,X). The natural inclusion function F1 overwraps the f(x) function. For

Table 2 Numerical results of test cases 1 and 2 with and without propagation. Test case 1

Without contractors With contractors

Test case 2

Iterations #

CPU time (s)

Optimum

Iterations #

CPU time (s)

Optimum

Not solved 2 144 374

17.78

17.014017

Not solved 74

0.12

 22.090756

412

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

3 Optimizers

2.5 zoom

X2

2

1.5

1

0.5

Bound test Constraint propagation Optimizers

0 0

0.2 0.4 0.6

0.8

1 X1

1.2 1.4

1.6 1.8

2

Fig. 7. Test case 2: search domain and constraint propagation.

more general constraint structure, we use better interval extensions: mean value forms, Taylor forms (first and second orders), etc. Continuous global optimization algorithms rely on interval techniques for dealing with numerical errors. A consequence of these numerical errors is the amplification of the variable domains and poor pruning results. Two major sources of numerical errors are operations with large numbers and operands with different magnitudes. Making some variable substitutions and scaling both objective function and constraint system may allow us to improve the algorithm performance. For example, instead of using constraint ð1010 x2 þx þ 1010 ¼ 0Þ, consider (x¼1010y) and use the constraint ðy2 þy þ 1 ¼ 0Þ. Constraint optimization solvers rely on the efficiency of the branch and prune algorithms to enforce consistency on the CSP variables. Precision and efficiency may be improved if the number of variables is reduced. Sometimes a set of constraints may be rewritten and reformulated into an equivalent system with fewer variables. The idea of reformulating the model came from a weak point of the interval branch and prune algorithms. In fact, there are two mechanisms that strongly impact the efficiency and performance of the optimization algorithm: the contracting step and the bound updating step. After testing the interval-based optimization algorithm, enhancements can be made by improving the upper bound updating step. The reformulation aims at stimulating this step, thus allowing the algorithm to discard more parts of the initial search box with the bound test. In order to understand the degrees of freedom based reformulation, let us look at the structure of engineering design models. Let P ¼ ðP 1 ,P2 ,. . .,P n Þ be a box representing the variables involved in the engineering model. We assume that the constraints have the following form: 8 Pi ¼ f i ðPÞ > > > > < Pj ¼ f ðPÞ j ... > > > > : P ¼ f ðPÞ k k

For instance, a preliminary design model could have the following structure: 8 minf ðPÞ > > > P A P0 > > > > < sc : P1 ¼ f 1 ðPÞ P ¼ ðP 1 ,P 2 ,P3 ,P 4 ,P5 ,P 6 Þ P2 ¼ f 2 ðPÞ > > > > P3 ¼ f 3 ðPÞ > > > : P ¼ f ðPÞ 4

4

In the classic algorithm, the upper bound of the global minimum is updated by testing if the center of the current box respects the equality constraints. In fact, every parameter is instantiated independently by the middle of its corresponding interval. This permits to ensure the bound constraints, and only afterwards, the equality constraints are tested. With the reformulation, the approach will be reversed. We will first ensure the satisfaction of the equality constraints, which are very tight, and then we will test if the variable bound constraints are satisfied, which are easier to fulfill. Since this paper deals with explicit systems, the parameters can be classified into two sets: Input parameters and Output parameters. In the reformulation step of this example illustrated in Fig. 8, the outputs are P1, P2, P3 and P4. The inputs are the remaining parameters P5 and P6. Indeed, for a given input set, the outputs can be cascade-computed via the constraints. For this, the constraints must be addressed in a specific order. For this example, we consider that the suitable order is f1, f2, f3 then f4. Thus, f1 is necessarily a function of inputs only, f2 is function of inputs and output of f1, f3 is function of inputs, output of f1 and output of f2, etc. To summarize, the reformulation consists in two main phases, each one is composed itself of two smaller steps. The first phase is to reduce the problems related intrinsically to the use of interval computation. This could be done by factorizing the constraints in order to reduce as much as possible the variables multi-occurrence, and by normalizing the variables sizes so as to obtain

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

413

Fig. 8. Illustration of the reformulation step.

homogeneous intervals and thus reducing the rounding errors and avoiding splitting the search space by the same variables (the largest ones). The second phase is to reduce the model variables by treating only the degrees of freedom. This could be done by classifying the constraints in the appropriate order that allows computing the output parameters using the input parameters. Also, it is necessary to ensure that the outputs vary in large intervals in order to make the feasibility test easier. Fig. 8 summarizes this phase. We denote by I the set of the input parameters and by O the set of the output parameters. The set O(fi) represents the output parameters of the fi function, for instance O(f1)¼P1. Let nv be the total number of variables (design parameters), and nc the total number of constraints. As the constraints are independent, the problem has (nv  nc) degrees of freedom, which corresponds to the input parameters. Since the outputs could be directly deducted from the inputs, the optimization problem can be reduced to only determining the optimal inputs. Consequently, the algorithm could process boxes of only (nv nc) dimensions. It should be noted that the reformulated problem is theoretically the same in terms of final results as the initial problem, but much simpler. According to the classification of the mathematical programming reformulations proposed by Liberti (2009), this reformulation is an opt-reformulation as all the optimums are preserved. This step considerably reduces the complexity of the problem and increases the algorithm performances. This can be explained by two main reasons: – processing smaller boxes mechanically reduces the complexity of the problem, as we can reach better precision in less computing time. This would obviously make the convergence of the algorithm much faster, since we reach the epsilon precision rapidly, – the bound updating step takes advantage of this reformulation. In fact, processing only the degrees of freedom makes the step of the upper bound updating much easier. Two steps of the constrained optimization algorithm previously presented must be revised in order to take into consideration the reformulation step: the bisection and the bound updating steps. The bisection is processed only on the input parameters. At each iteration, the bisection is performed by the largest side among the input parameters. For example, if the current box is B ¼([5y9],[1y5],[2y8],[3y4],[2y4],[1y2]), the resulting boxes of the bisection step are the two boxes B1 and B2 as the input parameters are P5 and P6: B1 ¼ ð½5. . .9,½1. . .5,½2. . .8,½3. . .4,½2. . .3,½1. . .2Þ B2 ¼ ð½5. . .9,½1. . .5,½. . .8,½3. . .4,½3. . .4,½1. . .2Þ The bound updating step is no longer done by the center of the whole box. It is however done by taking only the center of the input parameters box, and computing the corresponding output parameters. This step ensures that the resulting output parameters satisfy all the analytical constraints except the bound

Table 3 Reformulation of test cases 1 and 2. Test case 1 minf (ðxÞ ¼ x3 þ ðx1 þ x2 þ x3 Þx1 x4 x21 þ x22 þ x23 þ x24 ¼ 40 st : x1 x2 x3 x4 Z 25

Test case 2 minf ðxÞ ¼ 12x1 7x2 þ x22 st : 2x41 þ 2x2 ¼ 0 X A ½0,2  ½0,3

X A ½1,54

Reformulated Test case 1 min(f ðxÞ ¼ q x3ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ ðx1 þ x2 þ x3 Þx1 x4 st :

x1 ¼

40x22 x23 x24

x5 ¼ x1 x2 x3 x4

Reformulated Test case 2 min f ðxÞ ¼ 12x1 7x2 þ x22 st : x2 ¼ 2x41 þ 2 X A ½0,2  ½0,3

X A ½1,54  ½25, þ 1½

constraints. Indeed, we only have to check if these resulting output parameters belong to the search space. For the test case 2 introduced earlier, if we reformulate the problem as follows: min f ðxÞ ¼ 12x1 þðx2 7=2Þ2 49=4 subject to : x2 ¼ 2x41 þ2 X A ½0. . .2  ½0. . .3 Mathematically, this problem remains the same. The objective function f(x) has been rewritten into equivalent that minimizes the dependency problem, range(f(x)) in this case. The output parameter is x2, the input is x1. If the current box is ([0y1],[0y2.5]), take the middle of the input variable x1 ¼ 0:5 and compute the output variable x2 ¼ 2ð0:5Þ4 þ 2 ¼ 2:125. It is obvious that the vector (.500,2.125) is an instantiation of the variables that respects the equality constraint as it is computed in light of it. However, we have to verify that the output variable x2 respect its bound constraint as the input variable x1 respects by construction its bounds. Here we have to verify that 2.125 belongs to [0y3]. First, we tested the new algorithm on test cases 1 and 2. We added the reformulation step and made a comparison with the previous results (Tables 3 and 4). We noted marked improvement in the number of iterations and in processing time as in Table 5. In fact, we noted a gain of about 75% for the test case 1 and 40% for the test case 2. The obtained optimums are almost the same, but we noted an improvement in the number of updates that passed from 2 to 43 for the test case 1 and from 1 to 23 for the test case 2.

5. Application to the optimal design of electromagnetic actuator 5.1. Model presentation The efficiency of the reformulation has been tested on a more realistic design problem: a rotating machine with magnetic effects (Fig. 9). The detailed description of the problem and the

414

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

constraints are available (Kone et al., 1993). The problem consists on dimensioning a slot-less rectangular waveform permanent magnet machine called MAPSE. Many aspects are present in this problem and we can find constraint parameters from several engineering fields. The MAPSE design analytical model, established from simplified electromagnetic modeling, corresponds to various equations ensuring the flux conservation in the magnetic circuit, the electromechanical conversion and the overall heating up due to the losses in the stator by Joule effect (Messine et al., 1998). The electromagnetic torque, calculated by integrating into the whole zone the product of no-load magnetic flux density by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi current, is given byGem ¼ ðp=2lÞð1K f Þ kr bEch ED2 ðD þ EÞBe , where D is the bore diameter, E is the winding thickness, kr is the fitting factor, l is the magnet machine form factor (ratio of D to the length of iron L), l ¼D/L, b is the polar arc factor, Kf is the coefficient of the magnet leakage resulting from the relative importance of the magnetic air gap inherent in the slot-less structure, Be is the amplitude of the no-load-flux density at the level of the bore diameter. The machine global winding heating

the magnetic leakage coefficient Kf and the machine geometric dimensions has been established using the finite difference method: Kf ffi1.5pb((eþE)/D). The no-load magnetic field Be in the air gap, supposed exclusively radial, is given by: Be ¼ ð2la M=ðD log½ðD þ 2EÞ=ðD2ðla þ eÞÞÞÞ, where la is the thickness of the permanent magnets and M their magnetic polarization. An expression of the thickness C of the magnetic field in the iron Bfer is given. The inter polar leakage and armature reaction flux are assumed negligible allowing the conservation of the flux emitted by the poles: C ¼(pbBe/4pBfer)D. The number of pole pairs p is directly proportional to the bore diameter D and inversely proportional to the pitch pole Dp: p ¼ pD/Dp. This structural model is multi-physic. In fact, it combines constraints from several engineering fields: electrical, magnetic, mechanical and heat. It consists exclusively of equality constraints which are non-linear and explicit. The model parameters and their bound constraints (including user specifications) are summarized in Tables 6 and 7.

Ech is defined by the relation: Ech ¼ kr EJ 2cu , where Jcu represents the current areal density and krE is the current electric loading (kr is the coil winding filling factor and E the winding thickness). If e indicates the mechanical air gap, an empirical relation between

Table 6 The interval parameters of the design model.

Table 4 Classification of the parameters into inputs and outputs.

Inputs (freedom degrees) Outputs

Test case 1

Test case 2

x2,x3,x4 x1,x5

x1 x2

Variable parameters

Intervals

The bore diameter, D (m) The magnetic field in the air gap, Be (T) A semi-empiric magnetic leakage coefficient, Kf The current areal density, J cu ðA=m2 Þ The thickness of the mechanical air gap, e (m) The thickness of the permanent magnets, la (m) The winding thickness, E (m) The thickness of yoke, C (m) The polar arc factor, b Length, l

[.001,.5] [.1,1.0] [.01,.3] [105,107] [.001..005] [.001,.05] [.001,.05] [.001,.05] [.8,1.0] [1,2.5]

Table 5 The gain obtained by adding the reformulation step. Test case 1

Without reformulation With reformulation Gain (%)

Test case 2

Iterations #

CPU time (s)

Optimum

Iterations #

CPU time (s)

Optimum

2 144 374 1 148 081 46

17.780 7.142 60

17.014017 17.014017 0

74 43 42

0.123 0.115 4

 22.090756  22.090755 0

×/D D

la

C

Stator yoke Winding Mechanical gap Magnets Rotor yoke

C

e E

L Fig. 9. MAPSE bi-dimensional structure (Kone et al., 1993).

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

The user specifications consist in dimensioning a machine capable of developing a Gem ¼ 10.0 N m electromagnetic torque using rare earth-type magnets (M ¼ 0:9 T) and laminations that can stand a flux density of Bfer ¼1.5 T. Because of the form of the used conductors and insulators, the coil winding is fixed at 70% (kr ¼.7) filled by copper such that the global winding heating up is Ech ¼1011 A/m. The number of pole pairs is fixed to p ¼4 and the pitch pole must be in the region of 100 mm (Dp ¼.1 m). Because of manufacturing constraints, the mechanical air gap must be more than or equal to 1 mm (eZ.001 m). Inter polar leakages should be limited to 30% of the total flux generated by the inductor (Kf r.3). In this case study, the machine must be sized in order to minimize the active part volume Vu: Min V u ¼ p

D

l

ðD þ Eela Þð2C þ E þe þ la Þ

The structural design model that should be solved is defined as aggregation of the physico-mathematical model and the user specifications. These two parts are mixed to draw up the final design model corresponding to a global constrained optimization problem. The model consists of 17 parameters including 7 fixed parameters. Consequently, we have to determine 10 parameters. The

Table 7 The fixed parameters of the design model. Fixed parameters

Values

The The The The The The

0.7 1.5 1011 10 0.9 0.1

coil winding filling factor, kr magnetic field in the iron, Bfer (T) machine warm-up, Ech (A/m) electromagnetic torque, Gem (N m) magnetic polarization, M (T) polar step, Dp (m)

415

model has 6 independent equality constraints. Therefore it has 4 degrees of freedom (input parameters) (Fig. 10). After the reformulation step, the global optimization model has the following form: D Min V u ¼ p ðD þ Eela Þð2C þ Eþ e þla Þ 8l p Dp > > p >D¼ q > ffiffiffiffiffi > > > > J cu ¼ kErchE > > > > > > Be ¼ 2la M  < D þ 2E D log D2ðl a þ eÞ subject to > eþE > K ¼ 1:5p b > f D > > pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > l ¼ 2Gpem ð1K f Þ kr bEch ED2 ðD þEÞBe > > > > > pbBe > : C ¼ 4pB D f er In the reformulated model, we can see that the first and the second constraints are function of only inputs, the third is function of the inputs and the outputs of the previous equations and so on. This order allows us to compute iteratively the outputs for given inputs. In addition, we made sure that the outputs vary in interval so as to make easier the feasibility test. Within this reformulation, the structural model contains 4 degrees of freedom (input variables): e, la, E and a# . The remaining variables (D, Jcu, Be, Kf, and C) are output variables. In order to evaluate the impact of the reformulation, we introduced a feasibility coefficient a. The feasibility coefficient a is defined experimentally by calculating the percentage of feasible solutions among 1 000 000 randomly generated solutions from the search space. We summarize the obtained results in Table 8. We noted marked improvement in the coefficient a in the reformulated model from almost 0% to 5.1121%. The reformulated model remains mathematically the same, and as a result the constraints have the same mathematical properties. In fact, we have the same number of linear inequality constraints (LInC), non-linear inequality constraints (NLInC), linear equality constraints (LEC) and non-linear equality constraints (NLEC). Combined with the constraint propagation techniques, the reformulation impact would be more significant. In the next section, the test results using the IBBA algorithm will be detailed. 5.2. Numerical results

Fig. 10. The classification of the parameters into input and output.

Table 8 The impact of the reformulation on the feasibility coefficient. Problem

Var #

Obj-f

LInC #

NLInC #

LEC #

NLEC #

a (%)

Actuator

17

2

0

0

6

0.000

Reformulated actuator

17

Nonlinear Nonlinear

2

0

0

6

5.112

The electromagnetic actuator optimization model is solved using the IBBA based on constraint propagation. We have compared the two approaches: without and with reformulation. The results are summarized in Table 9. The table shows the definitive upper bound of the global minimum, the corresponding CPU time and iteration number. Also, it shows the number of updates and feasible solutions found during the algorithm execution. These results show that the updating of the upper bound is much easier with reformulation; fewer iterations and less execution time. It should be noted that, without reformulation, a tolerance of Et ¼ 1:0e3 has been added to the equality constraints in order to make the update step faster. Consequently, an equality constraint pi ¼f(P) is replaced by pi A ½f ðPÞEt . . .f ðPÞ þ Et . As said earlier, these equality constraints are the expression of the

Table 9 Main results.

a (%) Non reformulated model Reformulated model

0 5.11

Solution (1e  4) 6.0743 6.0736

Relaxation 10 0

3

Iterations #

CPU time (s)

Total updates #

Feasible solutions found #

60 151 649 470 960

1123 13

12 34

124 434 705 461

416

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

physical laws and thus their relaxation is not recommended. With the reformulation, no relaxation has been made. The solutions found respect scrupulously the equality constraints. Figs. 11 and 12 show that, without reformulation, the first update of the upper bound of the global minimum is done after

14 142 iterations, despite the constraints relaxation. On the other hand, the reformulation makes this step much easier. From the very fourth iteration, the algorithm has succeeded in finding a feasible solution. The update step has succeeded, almost instantly, in bringing the upper bound very close to the global minimum.

Fig. 11. Upper bound updates without reformulation: ] iterations.

Fig. 12. Upper bound updates with reformulation: ] iterations.

I. Mazhoud et al. / Engineering Applications of Artificial Intelligence 25 (2012) 404–417

The best solution is found after 19 781 585 iterations without using the reformulation. With the reformulated model, the best solution, which is slightly better and more accurate, is found after 153 067 iterations. This increase in the algorithm convergence is explained by the bound test that permits to discard more boxes since the upper bound is updated more often. We move from 12 updates (in 60 151 649 iterations) to 34 updates (in 470 960 iterations). The upper bound is updated more often simply because the algorithm finds feasible solutions easily. The number of feasible solutions found increases from 124 434 to 705 461.

6. Conclusion Deterministic methods are very interesting in the way that they guarantee the global optimum’s enclosure. Their main drawback is that they may be time-consuming. The purpose of this paper is to introduce an adapted optimization procedure that permits to: – find and enclose the global optimum and all its optimizers with a precision fixed by the user, – display reasonable runtime in solving engineering preliminary design problems. The interval analysis guarantees the enclosure of the global optimum and all its optimizers. This allows us to reach the first goal. The accelerating devices have helped us to reach the second goal. By adding a reformulation step, which came from using particular characteristics of many preliminary design models, the upper bound updating of the global minimum has become much easier. We succeeded to diminish notably the number of model evaluations, which constitute the greater part of the total runtime, from millions to thousands. In addition, no relaxation is required to make the update of the upper bound possible. This is all the more important as the analytical model consists of strong constraints. Many tests have been conducted in order to illustrate the efficiency of the algorithms and the reformulation step. Among these tests, we have validated our solution on the electromagnetic actuator design model. The obtained results correspond to the fixed targets. Naturally, many other improvements have to be made. The next steps revolve around two axes: – Theory: Improve the bound updating step and the quality of the inclusion functions used in the contractors. Also, many papers show the efficiency of the affine arithmetic and algorithms have been proposed (Messine, 2002), – Implementation: Currently every design model must be implemented separately and ‘‘manually’’ with its corresponding contractors. In order to make the algorithm introduced in this paper more practical, it is essential to automate these operations. Actually, we are testing our algorithm on other real engineering design problems such as a motorized hatch for the automotive industry, a heatsink for the laptop industry and so on. Also, we are working on global stochastic methods, such as PSO or Ant Colony algorithm, on real design problems. The work is in progress and, at this stage, we have only preliminary results. Our main purpose is to put our examples in benchmarks in order to make comparisons with other research teams working on the same subject.

417

References Ahn, Y.K., Kim, Y.C., Yang, B.S., 2003. Optimal design of a squeeze film damper using an enhanced genetic algorithm. J. Mech. Sci. Technol. 17 (12), 1938–1948. Aldanondo, M., Vareilles, E., Hadj-Hamou, K., Gaborit, P., 2008. Aiding design with constraints: an extension of quad trees in order to deal with piecewise functions. Int. J. Comput. Integr. Manuf. 21 (4), 353–365. Baumann, E., 1988. Optimal centered forms. BIT 28 (1), 80–87. Benhamou, F., Goualard, F., Granvilliers, L., Puget, J.-F., 1999. Revising Hull and Box consistency. In: Proceedings of the International Conference on Logic Programming, pp. 230–244. Byrne, R.P., Bogle, I.D.L., 1996. An accelerated interval method for global optimization. Comput. Chem. Eng. 20 (1), S49–S54. Casado, L.G., Garcia, I., Csendes, T., 2000. A new multisection technique in interval methods for global optimization. Computing 65 (3), 263–269. Chabert, G., Jaulin, L., 2009. Hull consistency under monotonicity. Lect. Notes Comput. Sci., 188–195. DPT, 2010. Design processing technologies: /http://designprocessing.free.frS. Gero, J.S., 2001. Mass consumption of creative designs. In: Proceedings of the International Conference on Engineering Design ICED’01, Glasgow. Gol, O., Sobhi-Najafbadi, B., 2005. Deterministic techniques for design optimization of electromotion devices with multiple objectives. J. Mater. Process. Technol. 161, 282–287. Hansen, E., 1979. Global optimization using interval analysis—the one-dimensional case. J. Optim. Theory Appl. 29 (3), 331–344. Hansen, E., 1980. Global optimization using interval analysis—the multidimensional case. Numer. Math. 34 (3), 247–270. Hansen, E., Walster, G.W., 2004. Global Optimization Using Interval Analysis. Marcel Dekker, New York. He, Q., Wang, L., 2007. An effective co-evolutionary particle swarm optimization for constrained engineering design problems. Eng. Appl. Artif. Intell. 20 (1), 89–99. Ho, W., Ho, G.T., Ji, P., Lau, H.C., 2008. A hybrid genetic algorithm for the multidepot vehicle routing problem. Eng. Appl. Artif. Intell. 21, 548–557. Isfahani, A.H., Vaez-Sadeh, S., 2007. Design optimization of a linear permanent magnet synchronous motor for extra low force pulsations. Energy Convers. Manage. 48, 443–449. Jaulin, L., Chabert, G., 2010. Resolution of nonlinear interval problems using symbolic interval arithmetic. Eng. Appl. Artif. Intell. 23 (6), 1035–1040. Jaulin, L., Kieffer, M., Didrit, O., Walter, E., 2001. Applied Interval Analysis with Examples in Parameter and State Estimation, Robust Control and Robotics. Springer-Verlag. Kearfott, R.B., 1992. An interval branch and bound algorithm for bound constrained optimization problems. J. Global Optim. 2 (3), 259–280. Kearfott, R.B., 1996. Interval computations: introduction, uses and resources. Euromath Bull. 2 (1), 95–112. Kone, A.D., Nogarede, B., Lajoie Mazenc, M., 1993. Le dimensionnement des actionneurs electriques: un proble me de progarammation non line´aire. J. Phys. 3, 285–301. Liberti, L., 2009. Reformulation in mathematical programming: definitions and systematics. RAIRO Oper. Res. 43 (1), 55–85. Lin, Y.H., Lin, S.C., 2001. Optimal weight design of rotor systems with oil-film bearings subjected to frequency constraints. Finite Elem. Anal. Des. 37, 777–798. Messine, F., 2002. Extensions of affine arithmetic: application to unconstrained global optimization. J. Universal Comput. Sci. 8 (11), 992–1015. Messine, F., 2004. Deterministic global optimization using interval constraint propagation techniques. RAIRO Oper. Res. 38 (4), 277–294. Messine, F., Nogarede, B., Lagouanelle, J.L., 1998. Optimal design of electromechanical actuators: a new method based on global optimization. IEEE Trans. Magn. 34 (1), 299–308. Moore, R.E., 1966. Interval Analysis. Prentice-Hall, Englewood Cliffs. Pahl, G., Beitz, W., Fddhusen, J., Grote, K.H., 2007. Engineering Design : A Systematic Approach, third ed Springer-Verlag. Pierre, S., Hyppolite, M.A., Bourjolly, J.M., Dioume, O., 1995. Topological design of computer communication networks using simulated annealing. Eng. Appl. Artif. Intell. 8 (1), 61–69. Qureshi, J., Dantan, J.-Y., Bruyere, J., Bigot, R., 2010. Set based robust design of mechanical systems using the quantifier constraint satisfaction algorithm. Eng. Appl. Artif. Intell. 23 (7), 1173–1186. Ratschek, H., Rokne, J., 1988. New Computer Methods for Global Optimization. Ellis Horwood, Chichester. Ratschek, H., Voller, R.L., 1991. What can interval analysis do for global optimization. J. Global Optim. 1 (2), 111–130. Roters, H.C., 1941. Electromagnetic Devices. John Wiley & Sons Inc.. Vinko, T., Lagouanelle, J.-L., Csendes, T., 2004. A new inclusion function for optimization: Kite—the one dimensional case. J. Global Optim. 30 (4), 435–456. Yvars, P.A., 2009. A CSP approach for the network of product lifecycle constraints consistency in a collaborative design context. Eng. Appl. Artif. Intell. 22 (6), 961–970.