New global optimization methods for ship design problems - CiteSeerX

Optim Eng (2009) 10: 533–555 DOI 10.1007/s11081-009-9085-3

New global optimization methods for ship design problems Emilio Fortunato Campana · Giampaolo Liuzzi · Stefano Lucidi · Daniele Peri · Veronica Piccialli · Antonio Pinto

Received: 1 January 2008 / Accepted: 6 April 2009 / Published online: 15 April 2009 © Springer Science+Business Media, LLC 2009

Abstract The aim of this paper is to solve optimal design problems for industrial applications when the objective function value requires the evaluation of expensive simulation codes and its first derivatives are not available. In order to achieve this goal we propose two new algorithms that draw inspiration from two existing approaches: a filled function based algorithm and a Particle Swarm Optimization method. In order to test the efficiency of the two proposed algorithms, we perform a numerical comparison both with the methods we drew inspiration from, and with some standard Global Optimization algorithms that are currently adopted in industrial design optimization. Finally, a realistic ship design problem, namely the reduction of the amplitude of the heave motion of a ship advancing in head seas (a problem connected to both safety and comfort), is solved using the new codes and other global and local derivative-

This work has been partially supported by the Ministero delle Infrastrutture e dei Trasporti in the framework of the research plan “Programma di Ricerca sulla Sicurezza”, Decreto 17/04/2003 G.U. n. 123 del 29/05/2003, by MIUR, FIRB 2001 Research Program Large-Scale Nonlinear Optimization and by the U.S. Office of Naval Research (NICOP grant N. 000140510617). E.F. Campana () · D. Peri · A. Pinto INSEAN—Istituto Nazionale per Studi ed Esperienze di Architettura Navale, Via di Vallerano 139, 00128 Roma, Italy e-mail: [email protected] G. Liuzzi Consiglio Nazionale delle Ricerche, Istituto di Analisi dei Sistemi ed Informatica “A. Ruberti”, Viale Manzoni 30, 00185 Roma, Italy S. Lucidi Dipartimento di Informatica e Sistemistica “A. Ruberti”, Università degli Studi di Roma “Sapienza”, Via Ariosto 25, 00185 Roma, Italy V. Piccialli Dipartimento di Ingegneria dell’Impresa, Università degli Studi di Roma “Tor Vergata”, Via del Policlinico 1, 00133 Roma, Italy

534

E.F. Campana et al.

free optimization methods. All the numerical results show the effectiveness of the two new algorithms. Keywords Nonlinear programming · Global optimization · Simulation based design

1 Introduction Simulation Based Design (SBD) of complex systems is an emerging tool to improve time-intensive industrial design which combines complex simulation codes based on fine computational meshes and algorithms for numerical optimization. In particular, in the fluid-dynamic design of marine, aeronautical, and automotive transport systems, the shape plays a key role and its detailed analysis often requires the solution of nonlinear partial differential equations (PDE), namely the NavierStokes equations, which are particularly expensive from a computational point of view in case of a realistic three dimensional geometry. The use of numerical codes able to solve this set of equations in SBD is allowed by the availability of high performance computing platforms. However, the cost of a simulation, i.e., an objective function evaluation, is CPU time consuming. A typical ship design analysis with Navier-Stokes equations takes as much as 10 hours per function evaluation for a three dimensional mesh, for example 2 × 106 nodes. Nonetheless, other applications may require more time. Under these circumstances the development of a SBD framework which combines these highly costly analysis tools and Global Optimization (GO) algorithms may appear as a paradox, but design engineers are tempted to take this direction (see Cox et al. 2001; Peri and Campana 2005), for at least three good reasons: (i) typical objective functions (e.g., total drag of cars, airplanes, ships, etc.) display, almost ever, noise, non-smoothness, and unavailability of derivatives, so that local optimizers could be trapped by local minima; (ii) the feasible design spaces of these optimization problems are very often nonconvex due to the presence of nonlinear geometrical and functional constraints that have to be enforced to prevent unrealistic results and provide a meaningful design; (iii) continuous experimental and computational activities have lead design engineers to produce near optimal configurations in many industrial fields, so that the margin for improvements is narrowing, and the probability that further improvements could come from local optimization methods is likely small. These open issues motivate the development of new derivative free global optimization algorithms which have great effectiveness and efficiency in terms of number of objective function evaluations. To this aim, we propose a new Particle Swarm Optimization (DDFPSO) method and a new filled function based approach (FILLDIR). We draw inspiration from Particle Swarm Optimization (PSO) methods because they are widely used for solving many industrial and engineering problems even though they have poor if no theoretical foundations (see for instance Venter and Sobieszczanski-Sobieski 2003, 2004; Pinto et al. 2004). On the other hand, we have

New global optimization methods for ship design problems

535

chosen filled function based algorithms because they have interesting theoretical properties (see Lucidi and Piccialli 2002; Xu et al. 2001), even though they are less used in practice. In the new algorithm DDFPSO we introduce a local search phase based on the derivative free minimization method (DF) described in Lucidi and Sciandrone (2002). As concerns FILLDIR we introduce a new filled function and we analyze its theoretical properties. Furthermore we propose a global optimization algorithm based on this new filled function, that roughly speaking is a deterministic multi-start approach. In particular, the filled function and the original function are minimized by means of the derivative free minimization method (DF) described in Lucidi and Sciandrone (2002). The new solvers are compared with a standard PSO method, with the filled function based algorithm proposed in Lucidi and Piccialli (2002), and also against standard methods which are currently used in industrial optimal design area on a set of well known test problems. The standard GO algorithms are chosen from different classes of GO methods, i.e., decomposition methods, genetic algorithms, clustering techniques and multi-start simulated annealing algorithms. Finally, we compare the enhanced GO algorithms with others global and local (DF described in Lucidi and Sciandrone (2002) and NEWUOA described in Powell (2004)) optimization methods on a real-world ship design problem, which consists in minimizing the vertical motion of the center of gravity (heave motion) of a ship advancing at fixed speed in head seas. We conclude this section by formally stating the constrained GO problem that arises from the ship design problem: min f (x) g(x) ≤ 0, l ≤ x ≤ u,

(1)

where x, l, u ∈ Rn , f : Rn → R and g : Rn → Rm are twice continuously differentiable functions, even though we assume that the derivatives are not available.

2 A new filled function approach In this section we propose a new filled function algorithm to solve problem (1). This problem has a twofold difficulty. On the one hand, the objective function presents multiple local minima and among them the global minimum point is sought. On the other hand, the presence of nonlinear constraints substantially increases the difficulty of the problem by introducing a further requirement on the solution points. As a first step toward the definition of an algorithm based on a filled function, we initially ignore all the constraints g(x) ≤ 0 and l ≤ x ≤ u, thus we refer to the problem min f (x).

x∈Rn

(2)

We require the following assumption to hold true in this section. Assumption 1 For any α ∈ R, the level set Lf (α) = {x ∈ Rn : f (x) ≤ α} is compact.

536

E.F. Campana et al.

This is a minimal requirement that is in general satisfied by any optimal design problem. The main idea of an algorithm based on a filled function is that, once determined a stationary point x1∗ of problem (2), a perturbation of the objective function, called filled function, is built such that its stationary points satisfy the relation f (x) < f (x1∗ ). In this way, by minimizing the filled function it is possible to find a stationary point of the filled function with a value of the original objective function lower than f (x1∗ ). More formally, a scheme of a general algorithm based on a filled function can be stated as follows: Algorithm FILLED Data: x0 ∈ Rn . Set k ← 0. Step 1. Compute xk∗ by applying a local minimization algorithm to the solution of problem (2) starting from xk . Step 2. Define a filled function Q(x, xk∗ ). Step 3. Choose a point x¯ and find y ∈ arg minn Q(x, xk∗ ) x∈R

by using a local minimization algorithm starting from x. ¯ Step 4. If f (y) < f (xk∗ ) then set xk+1 ← y, k ← k + 1 and goto Step 1. Otherwise goto Step 3. End Algorithm We do not specify a stopping criterion, since it can vary in practice from the application. For example, it can depend on the maximum allowed number of function evaluations, on the CPU time or on the numbers of failures at Step 4. In the following we introduce a new filled function which has different properties with respect to the one proposed in Lucidi and Piccialli (2002). In particular, the filled function in Lucidi and Piccialli (2002) has the favorable property to be globally convexized, but the price that must be paid is the existence of a known local minimum point worse than xk∗ that has a large basin of attraction, so that the local minimization in step 3 is often attracted by this local minimum. The filled function we introduce here is not globally convexized, but does not have stationary points at all where the objective function is higher than f (xk∗ ), and the numerical experiments prove that the new filled function is most reliable in practice to locate the global minimum point. 2.1 A new filled function In this subsection we introduce a new filled function which has the following expression: x − x˜ ∗ 2 Q(x, x˜ ∗ ) = exp − + 1 − exp (−τ [f (x) − f (x˜ ∗ ) + ]), (3) γ2 where x˜ ∗ is a known stationary point of the original objective function f (x), γ > 0 is a constant, and > 0 and τ ≥ 1 are real parameters.


537

This filled function is constituted by two terms; the first one, exp(−x − x˜ ∗ 2 / makes point x˜ ∗ a local maximum of the filled function Q(x, x˜ ∗ ) (drawing our inspiration from Renpu 1990 and Xu et al. 2001). The second term, 1 − exp(τ [f (x) − f (x˜ ∗ ) + ]), filters the stationary points of f (x) which have objective values greater or equal to f (x˜ ∗ ) and ensures that, for right values of the parameters, Q(x, x˜ ∗ ) has a local minimum point in a point with lower objective value than f (x˜ ∗ ). More in detail, it is possible to prove the following theoretical result: γ 2 ),

Proposition 2.1 There exists a τ¯ > 0 such that for all τ ≥ τ¯ the filled function has the following properties: (i) The point x˜ ∗ is an isolated local maximizer of the filled function Q(x, x˜ ∗ ). (ii) Q(x, x˜ ∗ ) has no unconstrained stationary point in {x ∈ Lf (f (x0 )) : f (x) ≥ f (x˜ ∗ )} except x˜ ∗ . (iii) If x˜ ∗ is not a global minimum of f (x) and satisfies the condition 0 < < f (x˜ ∗ ) − f (x ∗ ),

(4)

where x ∗ is a global minimum of f (x), then all the global minimum points xˇ of the filled function Q(x, x˜ ∗ ) over Lf (f (x0 )) belong to the region {x ∈ Lf (f (x0 )) : f (x) < f (x˜ ∗ )}. Proof First of all we note that the gradient of Q(x, x˜ ∗ ) has the following expression: ∇Q(x, x˜ ∗ ) = −2

(x − x˜ ∗ ) x − x˜ ∗ 2 exp − γ2 γ2

+ τ ∇f (x) exp (−τ (f (x) − f (x˜ ∗ ) + )).

(5)

We begin by proving point (i). Since the point x˜ ∗ is a stationary point of problem (2), it satisfies ∇f (x˜ ∗ ) = 0. Therefore, (5) implies that ∇Q(x˜ ∗ , x˜ ∗ ) = 0,

(6)

and hence x˜ ∗ is a stationary point of Q(x, x˜ ∗ ). Moreover, we have that ∇ 2 Q(x˜ ∗ , x˜ ∗ ) = −

2 I + τ ∇ 2 f (x˜ ∗ ) exp (−τ ) γ2

(7)

and hence it follows for all y ∈ Rn 2 y T ∇ 2 Q(x˜ ∗ , x˜ ∗ )y = y T − 2 I + τ ∇ 2 f (x˜ ∗ ) exp (−τ ) y γ 2 y2 + τy T ∇ 2 f (x˜ ∗ )y exp (−τ ) γ2 2 2 ∗ ≤ − 2 + τ λmax (∇ f (x˜ )) exp (−τ ) y2 , γ

=−

(8)

538

E.F. Campana et al.

where λmax (∇ 2 f (x˜ ∗ )) is the maximum eigenvalue of the matrix ∇ 2 f (x˜ ∗ ). The above inequality implies that there exists a τ1 > 0 such that, for all τ ≥ τ1 , the Hessian matrix ∇ 2 Q(x˜ ∗ , x˜ ∗ ) is negative definite. Therefore the point x˜ ∗ is an isolated local maximizer of Q(x, x˜ ∗ ) for all τ ≥ τ1 . As for point (ii), recalling the expression (5) of the gradient of Q(x, x˜ ∗ ) we note that if there would exists an unconstrained stationary point xˆ ∈ Lf (f (x0 )) of ˆ ≥ f (x˜ ∗ ), it must satisfy Q(x, x˜ ∗ ) such that xˆ = x˜ ∗ and f (x) xˆ − x˜ ∗ xˆ − x˜ ∗ 2 = τ ∇f (x) ˆ exp(−τ (f (x) ˆ − f (x˜ ∗ ) + )). 2 exp − γ2 γ2

(9)

Point (i) implies the existence of > 0 such that xˆ − x˜ ∗ > . By Assumption 1 the compactness of Lf (f (x0 )) follows. Therefore there exist two constants D and L such that xˆ − x˜ ∗ ≤ D and ∇f (x) ≤ L for all x ∈ Lf (f (x0 )). This implies the following estimates for the two sides of (9): xˆ − x˜ ∗ D2 xˆ − x˜ ∗ 2 2 2 exp − 2 ≤ 2 exp − γ γ γ2 γ2

(10)

τ ∇f (x) ˆ exp(−τ (f (x) ˆ − f (x˜ ∗ ) + )) ≤ τ L exp(−τ ).

(11)

and

Therefore (10) and (11) imply that there exists a τ2 ≥ τ1 such that for all τ ≥ τ2 condition (9) does not hold. Finally, we prove point (iii). The compactness of Lf (f (x0 )) implies the existence of a global minimum point xˇ ∈ Lf (f (x0 )) of Q(x, x˜ ∗ ). For all xˆ ∈ ∂Lf (f (x0 )) we have that f (x) ˆ = f (x0 ) ≥ f (x˜ ∗ ) since x˜ ∗ is obtained by a local minimization within the level set Lf (f (x0 )). Then xˆ − x˜ ∗ 2 + 1 − exp(−τ [f (x) ˆ − f (x˜ ∗ ) + ]) > 0. Q(x, ˆ x˜ ∗ ) = exp − γ2 Let x ∗ be a global minimum point of f (x). Then (4) implies that there exists a τ3 such that for all τ ≥ τ3 Q(x ∗ , x˜ ∗ ) < 0, from which we have Q(x, ˇ x˜ ∗ ) ≤ Q(x ∗ , x˜ ∗ ) < 0,

(12)

which implies that xˇ is in the strict interior of Lf (f (x0 )). Moreover, we have Q(x˜ ∗ , x˜ ∗ ) = 2 − exp (−τ ) > 0, and hence (12) implies that xˇ ∈ LQ (Q(x˜ ∗ , x˜ ∗ )). Therefore, since by point (i) we know that xˇ ∈ {x ∈ Lf (f (x0 )) : f (x) ≥ f (x˜ ∗ )}, and we proved that it is in the interior of Lf (f (x0 )), we get that xˇ ∈ {x ∈ Lf (f (x0 )) : f (x) < f (x˜ ∗ )} Thus, choosing τ¯ ≥ max{τ2 , τ3 }, the result follows.


539

The properties so far described hold for the unconstrained problem (2). The problem we are interested in, however, is problem (1) which presents some nonlinear constraints g(x) ≤ 0 and some box constraints l ≤ x ≤ u. In order to handle the constraints and obtaining a problem of form (2), we use an ∞ −penalty function. The resulting problem is minn p (x) = f (x) +

x∈R

+

1 max{0, g1 (x), . . . , gm (x)}

1 max{0, l1 − x1 , . . . , ln − xn , x1 − u1 , . . . , xn − un },

(13)

where > 0 is a penalty parameter which should be chosen sufficiently small in order for problem (13) being equivalent to problem (1). We note that function p (x) is non-smooth on the boundary of the set {x ∈ [l, u] : g(x) ≤ 0}, but two times continuously differentiable in its interior, where p (x) = f (x). Hence, under Assumption 1, p (x) is twice continuously differentiable in a sufficiently small neighborhood of every strictly feasible point x ∗ . Moreover, we note that since f (x) satisfies Assumption 1, then all the level sets of p (x) are compact as well. 2.2 A deterministic generation of starting points A critical aspect of algorithm FILLED resides in Step 3. Namely, depending on the starting point x¯ it could happen that point y produced by Step 3 does not satisfy condition f (y) < f (xk∗ ), because the sequence of points produced by the local minimization algorithm leaves the level set Lf (f (x0 )). Therefore, it could be necessary to reiterate Step 3 by choosing a different starting point x. ¯ The algorithm based on a filled function proposed in Lucidi and Piccialli (2002) solves this problem by extracting the needed starting points from a pseudo-random sequence. In this paper we propose a different selection rule based on a deterministic generation of points. In particular, we adopt a generation strategy, named DIRGEN, inspired by the DIRECT algorithm proposed in Jones et al. (1993). Whenever a new starting point is needed a single iteration of a DIRECT-type algorithm is performed. This amounts to subdividing one of the largest hyperintervals of the current partition of the feasible set using the same partitioning strategy as used by DIRECT. In this way it is guaranteed that the distance between successively generated starting points decreases quite uniformly. We stress the difference between the original DIRECT algorithm and the point generator DIRGEN, which consists in that DIRGEN does not perform any objective function based hyperinterval selection. In the limit, if applied an infinite number of times, DIRGEN produces a dense set of point in the box [l, u]. The algorithm which implements the above strategy and uses the new filled function applied to problem (13) is the following: Algorithm FILLDIR Data: l ≤ x0 ≤ u, > 0. Set k ← 0, k−1 ← ∅. Step 1. Compute xk∗ by applying a local minimization algorithm to the solution of problem (13) starting from xk .

540

E.F. Campana et al.

Step 2. Define a filled function x − x˜ ∗ 2 Q(x, xk∗ ) = exp − + 1 − exp(−τ [p (x) − p (x˜ ∗ ) + ]). γ2 Step 3. Generate a set of points k using DIRGEN. Set k ← k ∪ k−1 . Repeat. ¯ and find Choose a point x¯ ∈ k , set k ← k \ {x} y ∈ arg minn Q(x, xk∗ ) x∈R

by using a local minimization algorithm starting from x. ¯ Until p (y) < p (xk∗ ) or k = ∅. Step 4. If p (y) < p (xk∗ ) then set xk+1 ← y, k ← k + 1 and goto Step 1 ∗ ← xk∗ , k ← k + 1 and goto Step 3. else set xk+1 End Algorithm In our practical implementation, we set the constant γ = 1, the parameters τ = 100, ρ = 10−3 , = 10−3 . Moreover, since our final aim is to solve a problem where the derivatives are not available, we use as a local minimization tool the unconstrained derivative free minimization method DF proposed in Lucidi and Sciandrone (2002), that is the same used in the algorithm FILLED proposed in Lucidi and Piccialli (2002) we compare with.

3 Particle swarm optimization The Particle Swarm Optimization (PSO) algorithm is a recent addition to the list of global search methods (Russell et al. 2001). Since it was originally introduced in Kennedy and Eberhart (1995), the PSO algorithm has been studied by many authors (see Kennedy and Spears 1998; Parsopoulos and Vrahatis 2002; Shi and Eberhart 1998a, 1998b), and successfully applied to solve several engineering problems (see Pinto et al. 2004; Shi and Eberhart 2001; Venter and Sobieszczanski-Sobieski 2003, 2004). In the following, the standard PSO method is described and a new version is then introduced and tested. 3.1 Description of algorithm PSO The swarm strategy simulates the social behavior of a set of individuals (particles) which share information among them while exploring the design variables space. In PSO method each particle has its own memory to remember the best places that it has visited, whereas the swarm has a global memory to remember the best place ever visited. Moreover, each particle has an adaptable velocity to move itself in the design space. According to these principles, each particle investigates the search space analyzing its own flying experience and that of the other members of the swarm. The PSO algorithm is composed of the following four steps:


541

Step 0. (Initialization) Distribute a set of particles x0i inside the design space with random distribution and random initial velocities. Set k = 1. Step 1. (Compute velocity) Calculate a velocity vector v i for each particle, using the particle’s memory and the knowledge gained by the swarm according to i i b i vki = χ[wk vk−1 + c1 r1 (p i − xk−1 ) + c2 r2 (pk−1 − xk−1 )],

(14)

where χ is a constriction factor, w is called inertia weight, c1 and c2 are positive constants, r1 and r2 are random numbers equally distributed between b 0 and 1, p i is the best position found by the particle i and pk−1 is the best position found by the swarm up to iteration k − 1. Step 2. (Update position) Update the position of each particle, x i , using the velocity vector and previous position i xki = xk−1 + vki .

(15)

Step 3. (Check convergence) Set k = k + 1. Go to Step 1 and repeat until convergence. Literature reports that fine tuning of the parameters in (14) is crucial for the optimization process, and that the final solution and the calculation time are strictly linked to the parameters setting. The inertia weight w regulates the trade-off between the global (wide-ranging) and local (nearby) exploration abilities of the swarm. A large inertia weight facilitates global exploration (searching new areas), while a small one facilitates local exploration. A suitable value for w usually provides balance between global and local exploration abilities and consequently results in a reduction of the number of iterations required to locate the optimal solution. Experimental results indicate that it is better to initially set the inertia to a large value, in order to promote global exploration of the design variables space and gradually decrease it to get a more refined solution (Venter and Sobieszczanski-Sobieski 2003). For these reasons, an initial value for w is set and the decrease rate is calculated by wk = wk−1 gw ,

(16)

where wk is the new value for the inertia weight, wk−1 is the previous one and gw is a constant chosen between 0 and 1. In Venter and Sobieszczanski-Sobieski (2003) it is suggested to use 0.35 < w < 1.4, and gw = 0.975. The constriction factor χ is used alternatively to w to limit the maximum velocity. The major difference between the two is that while the inertia w is employed to control the impact of the previous history of velocities on the current one, χ offers to the user the chance to select the search resolution. Quantitatively, if box constraints are given, χ takes a value equal to a fraction of the characteristic dimension of the box. The constants c1 and c2 are called cognitive and social parameter, respectively. The cognitive parameter indicates how much confidence the particle has in itself, while the social parameter indicates how much confidence it has in the swarm. Proper fine-tuning of these coefficients may result in faster convergence and alleviation of

542

E.F. Campana et al.

local minima. In basic PSO algorithm (Kennedy and Eberhart 1995) the authors propose c1 = c2 = 2, so that the mean of stochastic multipliers of (14) is 1. In Venter and Sobieszczanski-Sobieski (2003) different values for the two coefficients are used. In particular, c1 = 1.5 and c2 = 2.5 work well in their examples. In Campana et al. (2006a, 2006b) a more rigorous analysis is carried out. In these papers a generalized PSO iteration is described by means of a dynamic linear system, whose properties are analyzed. The influence of the particles’ starting points and the use of deterministic or stochastic parameters are investigated and some partial convergence results are given. In particular, the PSO parameters are selected by imposing that the particles trajectories are confined in a suitable compact set. 3.2 A new PSO algorithm: DDFPSO This section focuses on the development of an enhanced PSO version, named Deterministic Derivative-Free Particle Swarm Optimization (DDFPSO). In the following all the modifications are described. (i) Initialization of the swarm: The use of GO algorithms based on expensive analysis tools imposes a substantial reduction of the particles number. Instead of using a random distribution, a deterministic one is proposed. In particular, at Step 0 the initial swarm is built with one particle at the center of each face of the ndimensional hyper-cube which represents the design space. As a consequence, the total number of particles is 2n. Moreover, the initial velocity vector, defined in (14), is set equal to 0. (ii) Boundary search phase: According to (14) and (15) each particle is attracted by the others. As a consequence, they cannot escape from the dashed region indicated in Fig.1 unless the inertia term is sufficiently strong, i.e., the corners of the design space cannot be reached by any particle of the swarm. In order to force the search along the boundaries of the design space we introduce a threshold value that limits the d-component of the particle’s speed orthogonal to the face where the particle has been initially placed: |vkd | ≤

d − xd | |xmax min , gd

g d < 1.0.

(17)

This threshold value is initially set as a fraction of the box’s dimension in the d-th direction. This limit is progressively relaxed during the iterations. By applying this normal speed limiter it is possible to explore the corners of the design space and find the correct global minimum, as reported in the example shown in Fig. 1. The dashed (inner) region of the feasible space is obtained by connecting the initial positions of the particles. Without the normal speed limiter (a), particles are strongly attracted each other and tend to be confined inside the dashed region, hence failing to locate the global optimum. The use of the normal speed limiter (b) allows the particles to explore the design space near the boundary too and to find the global optimum. (iii) Suppression of random coefficients: In the new version PSO is modified according to a deterministic flavor. Indeed, we decide to fix in (14) the parameters


543

Fig. 1 The path of the swarm particles for the n = 2 Griewank test function, with and without the normal speed limiter introduced with (17)

r1 and r2 equal to 1, thus we eliminate the random factor introduced by these two coefficients. In this way, we transform a pure stochastic method into a deterministic one (Pinto et al. 2004). The motivation stems from the use of PSO in combination with CPU time expensive numerical simulations used to obtain objective function and constraints information: a stochastic approach would require repeated runs which might require simply too much computing time for real industrial applications. (iv) Particles with violated constraints: The original PSO algorithm is defined only for unconstrained optimization problems. For this reason, when we are dealing with constrained problems we set equal to 0 the inertia weight term of particles with violated constraints (Venter and Sobieszczanski-Sobieski 2003; Pinto et al. 2004). As a consequence, (14) is substituted by i b i ) + c2 (pk−1 − xk−1 )]. vki = χ[c1 (p i − xk−1

(18)

In most cases the new velocity vector will point back to the feasible region. This feature of the DDFPSO algorithm is particularly useful in optimization problems with a non-convex feasible design space, as for example the ship design problem that will be discussed in the following sections (Pinto et al. 2004). (v) Convergence criterion for the global search phase: In the original algorithm there was no stopping criterion. For real applications, however, a maximum number of iterations is always fixed. We here introduce an heuristic convergence criterion for the PSO phase, used in the algorithm to switch from the global to the local search phase. For the convergence of the global phase we focus onto the identification of all the particles which fall into the same basin of attraction. In order to recognize a basin of attraction, the stored information about the previous steps of the algorithm, i.e., the visited locations and the corresponding objective function values, are necessary. We monitor two quantities:

544

E.F. Campana et al.

– At each iteration k, all the particles of the swarm are analyzed: if a descent trend (19) is shown for a given number of consecutive steps, the ith particle is marked as attracted by a basin. i )) · (f (x i ) − f (x i )) ≥ 0, (f (xki ) − f (xk−1 k−1 k−2 (19) i ) < 0. f (xki ) − f (xk−1 – Then the distance among the attracted particles is performed. These particles are assumed to be in the same basin if the maximum distance among them is a given fraction of the maximum distance among all the initial positions. – Furthermore, an average radius of the distance among these attracted particles is computed and the center of the basin is estimated. This radius is finally used to check the remaining (non attracted) particles of the swarm. Given their distances from the center of the basin and the remaining steps before the maximum iteration is reached, if their current speed is already too slow to bring them back to the basin where the current optimum is located, they are abandoned. When all the swarm particles are (i) in the basin or (ii) abandoned, the global search phase is ended. (vi) Local search phase: It has been observed that PSO is generally fast in the identification of the attraction basins but it is quite slow to converge. The strategy adopted here is to use a two phase global-local approach. Thus, we introduce the following new step to perform a local refinement of the solution: Step 4. (Local search) Starting from the point with lowest objective function value, we perform a local minimization with the derivative free linesearch method DF proposed in Lucidi and Sciandrone (2002).

4 Numerical results on a set of test problems The aim of this section is twofold. On the one hand we want to test whether the new methods (FILLDIR and DDFPSO) represent a significant improvement with respect to the existing ones. On the other hand, we compare the new solvers against a selection of well-known global optimization methods. First of all, a description of the test problems used in the numerical experience and comparison is in order. A set of 31 well-known test problems was selected (see Lucidi and Piccioni 1989), with dimensions ranging from n = 2 through n = 20. Table 1 reports the problem name, dimension, and the value of the known global minimum fglob . In an effort to obtain meaningful results, the same stopping criterion for all the methods was used. In particular, the iterations are stopped either when a point x¯ with f (x) ¯ − fglob ≤ 10−4 max{1, |fglob |}

(20)

has been located or when the number of objective function evaluations has exceeded a prefixed and dimension-dependent amount maxnf. Thus, two sets of results are

New global optimization methods for ship design problems Table 1 Test problems

Problem

545 n

fglob −1.03162

Six humps camel back

2

Treccani

2

Quartic

2

−0.3523

Schubert

2

−186.7309

Schubert pen.1

2

−186.7309

Schubert pen.2

2

−186.7309

Cosine mixture

2

−0.2

Cosine mixture

4

Shekel 5 loc. minima

4

−10.1532


4

−10.4029


4

−10.5364

Exponential

2, 4

−1.0

Hartman

3

−3.8627

Hartman

6

−3.3223

5n loc. minima

2, 5, 10, 20

0.0

10n loc. minima

2, 5, 10, 20

0.0

0.0

−0.4

15n loc. minima

2, 5, 10, 20

0.0

Griewank

2, 5, 10, 20

0.0

obtained by setting maxnf to 100n and 1000n, respectively. Furthermore, when a method makes use of some randomness, we report (for every problem) the average behavior on 100 runs. The results are presented showing figures which report the success rates and the average number of function evaluations, nf∗ , needed to get convergence on problems with the same dimension. Further, these values are normalized with respect to the prefixed maximum number of objective function evaluations. In this way we obtain values ranging from 0 to 1. To take into account those cases where condition (20) is not met, we simply sum up the maximum number of function evaluations allowed for that specific run. As concerns the selection of existing off-the-shelf global optimization solvers, we consider only methods which do not require any derivative knowledge on the problem and do not use any technique to numerically estimate first order derivatives. This derivative-free feature is an essential requirement in that derivatives can be neither computed nor approximated for the ship design problem that is our ultimate task. In particular, we consider methods for global optimization which are used in the optimal design area and which belong to the following categories: 1. Decomposition methods. Namely those methods which recursively subdivide the search space into smaller regions that, hopefully, should be finer around global minima and coarser away from them. 2. Evolutionary and genetic algorithms. The methods of this class try to solve the problem by exploiting the genetic behavior of a set of individuals which, through mutations and breeding, try to achieve the fittest individual ever possible which should correspond to the global minimum of the problem.

546

E.F. Campana et al.

3. Clustering techniques. The aim of the methods belonging to this class is that of locating all the local minima of the problem. In order to accomplish this task, they look for the basins of attraction of every local minimum. In particular, an initial random sample of points is divided into clusters each one containing points which are in the same basin. In this way, a local optimization started from one point for every cluster should be able to locate the corresponding different local minima with the best one being the global minimum point. 4. Multi-start algorithms. Methods of this class perform many local minimizations starting from different starting points. The starting points can be chosen either at random or according to some probabilistic criterion. Another aspect which characterizes multi-start algorithms is the particular local search method adopted which should depend on the problem to be solved in such a way to better exploit all the information available. In an attempt to get a significant numerical comparison with the three new methods, we selected one method from every class choosing it among those which are believed to be particularly efficient. Further, we privileged those algorithms for which a Fortran implementation is publicly available and, we do not consider any software designed for specially-structured problems. We first collected some codes available through the web page of Neumaier1 and in particular DIRECT (see Jones et al. 1993) and DE (see Price and Storn 1997), which won the third place at the 1st International Contest on Evolutionary Computation on a real-valued function testsuite. DE proved to be the best genetic algorithm approach. The first two places of the contest were won by non-GA algorithms. Moreover, we have chosen GLOBAL (see Boender et al. 1982) using the derivative free local optimizer UNIRANDI which is considered to be a robust implementation of a clustering-type algorithm. Finally, we use DDFSA as an example of a multi-start algorithm. This algorithm has been recently proposed in Liuzzi et al. (2004) for the solution of an optimal design problem. DDFSA adopts a simulated annealing criterion to generate suitable starting points for the local searches which are carried out in a distributed fashion by means of an efficient derivative-free procedure, as proposed in Lucidi and Sciandrone (2002). In Figs. 2 and 3 we report, for every solver, the number of function evaluations and the success rates, respectively, both when maxnf = 100n and maxnf = 1000n. Whichever value the maximum number of function evaluations takes, the methods which are the best are FILLDIR and DDFPSO; as concerns the other methods the rating changes quite a lot, in the sense that some methods (as DDFSA and DIRECT) become more and more competitive as the number of function evaluations increases. This behavior is due to the intrinsic distributed structure of these algorithms. It emerges that FILLDIR and DDFPSO require uniformly a lower number of function evaluations to converge. 5 A ship design problem: hull form optimization for seakeeping Aim of the design optimization problem solved in the following is to improve the seakeeping performance of a ship, i.e., the behavior of a ship in a seaway, a well 1 http://solon.cma.univie.ac.at/ neum/.


547

Fig. 2 Average number of function evaluations necessary to converge, normalized with maxnf (100n top, 1000n bottom). It should be noted that for some of the algorithms it was impossible to enforce exactly the maximum number of evaluations equal to maxnf

recognized issue in ship design. It is indeed a common experience that ship performance worsens in rough weather. Passengers and crew may be seasick, and to carry out any task on board takes longer than it does in calm weather. Operability of the ship may dramatically decrease as the vertical acceleration increases and eventually, for severe ship’s motions, the ship may capsize. Additionally, massive structural failures may occur. These are often the cause of ship’s sinking, even of large bulk carriers

548

E.F. Campana et al.

Fig. 3 Number of successes when maxnf is equal to 100n and 1000n

longer than three football pitches, wide as a six lane freeway and with a gross tonnage of 100,000. Broadly speaking, the ship designer takes into account the ability of the ship to achieve its mission in rough weather conditions. But while it is possible to reduce ship’s roll motions by applying some additional devices (i.e., roll fins), vertical motions may only be varied by changing the shape of the hull and hence the mass distribution. Therefore, the problem considered here is one of shape optimization.


549

5.1 The mathematical formulation of the seakeeping problem A frequent request from shipyards and navies is to limit the vertical motion and acceleration while the ship is advancing in irregular seas. The prediction of this behavior is typically obtained by studying the problem in the frequency domain under two major assumptions: it is assumed that (i) the irregular sea waves can be represented as a sum of a simple sine waves with specified amplitudes from selected spectral densities and random phases with uniform distributions and (ii) the response of the ship to irregular waves can be obtained as the sum of the ship responses to the individual sine components. In the following, the mathematical formulation for computing the ship motions is briefly summarized, whereas for more details the interested reader is referred to Salvesen et al. (1970) and Newman (1977). Consider a ship advancing at a constant forward speed U in regular sine waves and assume that the resulting motions are linear and harmonic. The ship is assumed to be a rigid body with 6 degrees of freedom and the motion’s amplitudes are indicated as ζj (j = 1, . . . , 6).2 Neglecting viscous effects, a velocity potential (x, y, z, t) can be introduced and the problem is formulated in terms of the potential flow theory. is then separated into two parts: the steady contribution due to the forward speed motion (−U x + φS (x, y, z)) and the time-dependent part associated with the incident wave system and unsteady body motions: (x, y, z, t) = −U x + φS (x, y, z) + φT (x, y, z)eiωt , where ω is the encounter frequency of the waves in the moving reference frame. The exact boundary conditions are then linearized considering only small oscillatory motions and decomposing the time-dependent potential as follows: φT = φI + φD +

6

ζj φj ,

j =1

where φI is the incident wave potential, φD the diffraction potential and φj the contribution to the velocity potential from the j th mode of motion. According to this analysis, the following linearized boundary conditions have to be satisfied • on the hull surface (at the mean position): ∂ [−U x + φS ] = 0, ∂n

∂φD ∂φI + = 0, ∂n ∂n

∂φj = iωnj + U mj = 0, ∂n

• on the free surface (z = 0): ∂ 2 ∂ 2 φS ∂φS ∂ = 0, iω − U (φI , φD ) = 0, + + g ∂n ∂x ∂z ∂x 2 ∂ 2 ∂ φj = 0. iω − U φj + g ∂x ∂z

U2

2 The index j = 1, . . . , 6 refers to the 6 degrees of freedom: surge, sway, heave, roll, pitch and yaw, respec-

tively.

550

E.F. Campana et al.

Finally, each of the potentials φS , φI and φD must satisfy the Laplace equation in the fluid domain and some radiation conditions at infinity. Once the problem is solved, the hydrodynamic forces can be obtained by integrating the pressure on the hull, given by the linearized, time-dependent Bernoulli’s equation: ∂ p = −ρ iω − U φT eiωt − ρg(ζ3 + ζ4 y − ζ5 x)eiωt . ∂x The numerical solution of the continuous problem described so far may be obtained via a Boundary Element Method in the frequency domain (e.g., Meyers et al. 1981). 5.2 Description of the design problem The objective function selected for this application is the peak of the Response Amplitude Operator (RAO) for the heave motion, namely the vertical motion of the center of gravity, when the ship is advancing at constant speed in head seas. For a given ship the RAO is defined as the square of the amplitude of the regular wave transfer function at each frequency of the incoming sea, i.e., the transfer function of the ship for all the six motions under any incoming wave of assigned relative direction. RAO is frequently used by the ship designers as a measure of the behavior of a ship when operating at sea. In the design phase of a ferry or a cruiser ship, emphasis will be placed upon stability, to ensure the comfort of the passengers. For a naval warship the focus will be on the maximum value of vertical motions and accelerations. In this paper, the analysis will be carried out for a containership while it is advancing at the speed of 16 knots in head seas. The minimum of the response is searched for nondimensional frequencies higher than 0.4 and for a ship speed of Fr = 0.198 where the Froude number Fr is defined as the ratio of inertial to gravitational forces. The S175 containership (Fig. 5) is selected as test case for the design problem. It has been adopted by the ITTC3 Seakeeping Committee as a benchmark test. The seakeeping performance of the S175 containership is numerically evaluated with a simulation code based on the strip theory, whose mathematical background is given in the previous section (see Newman 1977). With this approach the hull is represented as a series of two-dimensional slices and the model for the fluid is based on potential flow assumptions. Simulations based on this approach for computing heave motion for conventional mono-hulls at conventional speeds are of sufficient accuracy for design purposes. The output is the transfer function of the ship for all the six rigid body motions under any incoming wave of assigned relative direction. Hull shape parameterization is performed according to Peri et al. (2001) and a Béziér patch is superimposed to the original hull shape. Six control points of the patch are used as design variables to modify the hull shape. Consequently, we can assume that the shape of the ship hull depends on a design variable vector x = (x1 , x2 , . . . , x6 )T . The advantages of this parameterization technique are its great flexibility in terms of construction of new shapes and simplicity of implementation. The patch acts in the y 3 ITTC stands for Int. Towing Tank Conf., an international organization of the naval hydrodynamic com-

munity.


551

(transversal) direction only. As a consequence, the keel line is left unchanged during the optimization process. In order to come up with a problem having the characteristics of a real design case, we impose some geometrical constraints, which must be taken into account during the optimization process. Three different constraints are adopted, in particular we define a range of variation for the displacement, (x), and for the beam, B(x): 2398 ≤ (x) ≤ 2460,

(21)

25 ≤ B(x) ≤ 26,

(22)

where (x) is measured in tons and B(x) in meters. Finally, in order to avoid an unreliable hull form we impose a maximum range of variation for the design variables, xi : −20.0 ≤ xi ≤ 20.0,

i = 1, . . . , 6.

(23)

The imposition of the nonlinear constraints (21) and (22) creates a nonconvex feasible design space (see Pinto et al. 2004). Whenever the objective function must be evaluated on an infeasible point, we simply assign to it the dummy value +∞. Moreover, we mention that also for the ship design problem we use a maximum number of objective function evaluations as stopping criterion. As a result, two runs have been performed for each algorithm setting maxnf equal to 100n and 1000n. 5.3 Optimization of the S175 containership We performed the shape optimization of containership S175 using the enhanced versions, presented in previous sections, and DIRECT, which has proved to be the best among off-the-shelf codes. Values of the objective function, the improvement with respect the initial objective function value, and the number of function evaluations required by each method are reported in Table 2 for the sake of comparison. All the algorithms show good outcomes which are very similar to each other. In particular, when maxnf = 100n all the codes perform quite similarly except for FILLDIR where the improvement in the objective function value is lower than the others. In this setting, DIRECT achieves the best result. On the other hand, when maxnf = 1000n, all the methods find a better solution with FILLDIR having the biggest improvement. In this case, FILLDIR is also able to find the best solution. When maxnf = 1000n, we note that DDFPSO stops before the function evaluation limit is reached just because an inner stopping condition gets satisfied. DIRECT, which does not have an inner stopping condition other than the maximum function evaluations, also stops prematurely in 2345 evaluations. This is due to the constrained nature of the hull optimization problem. In particular, on this problem DIRECT generates nearly 90000 points which is the maximum allowed by the code of which only the reported 2345 happen to be feasible. Moreover, as concerns FILLDIR, the method does not have an inner stopping condition so that it continues iterating until the 1000n = 6000 function evaluations stopping condition is met. Nonetheless, it

552 Table 2 Comparative results for the RAO’s heave motion peak optimization (maxnf = 100n and maxnf = 1000n)

E.F. Campana et al. Algorithm maxnf = 100n fmin

maxnf = 1000n

improv. (%) nf

fmin

improv. (%) nf

FILLDIR 0.8796 33.36

320 0.8572 35.06

DDFPSO

0.8648 34.48

600 0.8645 34.50

982 714

DIRECT

0.8636 34.57

601 0.8627 34.64

2345

Fig. 4 Comparison of the heave motion RAOs of the original and optimized hulls (maxnf = 100n), where the minimum has been searched for non-dimensional frequencies higher than 0.4 only

finds the best solution of the problem after only 982 iterations. More in particular, it should be noted that DDFPSO cannot manage to significantly improve the objective function value when maxnf = 6000. Indeed, DDFPSO stops with 714 function evaluations returning a point with an objective function value which is higher than the value obtained by DIRECT when maxnf = 600. In Fig. 4 we report the RAOs associated to the various hull forms obtained with 600 objective function evaluations. The shape of all curves is significantly better than the starting one. In fact, the initial heave motion peak of the RAO has been damped. Furthermore, since the optima achieved by DIRECT and DDFPSO are very close, also their RAOs are similar. Figure 5 shows the original and the optimized hull forms of the S175 containership and the corresponding design variables values in the case of maxnf = 100n. As a consequence of the small differences in the numerical results among DIRECT and DDFPSO (i.e., design variable values and objective function optima) their corresponding optimized hull shapes are similar.


553

Fig. 5 S175 hull shape and optimum design variables comparison (maxnf = 100n)

6 Conclusions In this paper we propose new enhanced versions of a filled function approach and of a Particle Swarm Optimization method, with the aim of improving their effectiveness and efficiency. FILLED and PSO are enhanced in various ways. In particular, in DDFPSO we introduce a local search phase based on the derivative free minimization method DF. As concerns FILLED, which already uses the local minimization method DF, we introduce a new filled function and we analyze its theoretical properties. The new solvers are compared with respect to their standard versions on a set of well-known test problems. The results show that the new versions are more effective and efficient with respect to the original ones. Moreover, we compare the enhanced methods also against some codes chosen from different classes of GO methods and currently used to solve industrial optimization problems. Finally, we solve a realistic ship design problem using the enhanced methods and DIRECT. FILLDIR proves to be the best among off-the-shelf codes. The heave mo-

554

E.F. Campana et al.

tion peak of the Response Amplitude Operator of a containership is selected as objective function for the design problem. All the new GO methods proposed here and DIRECT code are able to find satisfactory designs which significantly improve the performance of the original geometry.

References Boender CGE, Rinnooy Kan AHG, Timmer GT, Stougie L (1982) A stochastic method for global optimization. Math Program 22:125–140 Campana EF, Fasano G, Pinto A (2006a) Dynamic system analysis and parameter selection in PSO, for costly applications. INSEAN Tech Rep 2006-019 Campana EF, Fasano G, Pinto A (2006b) Dynamic system analysis and initial particles position in particle swarm optimization. In: IEEE Symposium on swarm intelligence, Indianapolis (USA) Cox SE, Haftka RT, Baker CA, Grossman B, Mason WH, Watson LT (2001) A comparison of global optimization methods for the design of a high-speed civil transport. J Glob Optim 21:415–433 Jones DR, Perttunen CD, Stuckman BE (1993) Lipschitzian optimization without the Lipschitz constant. J Optim Theory Appl 79(1):157–181 Kennedy J, Eberhart R (1995) Particle swarm optimization. In: Proceedings of the IEEE international conference on neural networks, vol 4. IEEE Service Center, Piscataway, pp 1942–1948 Kennedy J, Spears WM (1998) Matching algorithms to problems: an experimental test of the particle swarm and some genetic algorithms on the multimodal problem generator. In: Proceedings of the 1998 IEEE international conference on Evolutionary Computation, Anchorage, Alaska Liuzzi G, Lucidi S, Piccialli V, Sotgiu A (2004) A magnetic resonance device designed via global optimization techniques. Math Program Ser B 101:339–364 Lucidi S, Piccialli V (2002) New classes of globally convexized filled functions. J Glob Optim 24:219–236 Lucidi S, Piccioni M (1989) Random tunneling by means of acceptance-rejection sampling for global optimization. J Optim Theory Appl 62(2):255–279 Lucidi S, Sciandrone M (2002) On the global convergence of derivative free methods for unconstrained optimization. SIAM J Optim 13:97–116 Meyers WG, Applebee TR, Baitis AE (1981) User’s manual for the standard ship motion program, SMP. David Taylor internal report DTNSRDC/SPD-0936-01 Newman JN (1977) Marine hydrodynamics. MIT Press, Cambridge Parsopoulos KE, Vrahatis MN (2002) Recent approaches to global optimization problems through particle swarm optimization. Nat Comput 1:235–306 Peri D, Campana EF (2005) High fidelity models in global optimization. In: Global optimization and constraint satisfaction. Lecture notes in computer science, vol 3478. Springer, Berlin, p 112 Peri D, Rossetti M, Campana EF (2001) Design optimization of ship hulls via CFD techniques. J Ship Res 45(2):140–149 Pinto A, Peri D, Campana EF (2004) Global optimization algorithms in naval hydrodynamics. Ship Technol Res 51(3):123–133 Powell D (2004) The NEWUOA software for unconstrained optimization without derivatives. In: 40th workshop on large scale nonlinear optimization, Erice, Italy Price K, Storn R (1997) Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces. J Glob Optim 11:341–359 Renpu G (1990) A filled function method for finding a global minimizer of a function of several variables. Math Program 46:191–204 Russell C, Eberhart RC, Shi Y (2001) Particle swarm optimization: developments, applications, and resources. In: Proceedings of the congress on evolutionary computation, pp 81–86 Salvesen N, Tuck EO, Faltinsen O (1970) Ship motions and sea loads. Trans SNAME 135:250–287 Shi Y, Eberhart RC (1998a) Parameter selection in particle swarm optimization. In: The 7th annual conference on evolutionary programming, San Diego, USA Shi Y, Eberhart RC (1998b) A modified particle swarm optimizer. In: Proceedings of the IEEE international conference on evolutionary computation, Anchorage, Alaska


555

Shi Y, Eberhart RC (2001) Particle swarm optimization: developments, applications and resources. In: IEEE congress on evolutionary computation, pp 27–30 Venter G, Sobieszczanski-Sobieski J (2003) Particle swarm optimization. AIAA J 41(8):1583–1589 Venter G, Sobieszczanski-Sobieski J (2004) Multidisciplinary optimization of a transport aircraft wing using particle swarm optimization. Struct Multidiscip Optim 26(1–2):121–131 Xu Z, Huang HX, Pardalos PM, Xu C (2001) Filled functions for unconstrained global optimization. J Glob Optim 20:49–65