Oscar Brito Augusto

Fouad Bennis
email: fouad.bennis@ec-nantes.fr

Stephane Caro
email: stephane.caro@irccyn.ec-nantes.fr

Oscar Brito

A NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS

Keywords: multi-objective optimization, multi-attribute decision making, engineering design optimization

Many engineering sectors are challenged by multi-objective optimization problems. Even if the idea behind these problems is simple and well established, the implementation of any procedure to solve them is not a trivial task. The use of evolutionary algorithms to find candidate solutions is widespread. Usually they supply a discrete picture of the non-dominated solutions, a Pareto set. Although it is very interesting to know the non-dominated solutions, an additional criterion is needed to select one solution to be deployed. To better support the design process, this paper presents a new method of solving non-linear multi-objective optimization problems by adding a control function that will guide the optimization process over the Pareto set that does not need to be found explicitly. The proposed methodology differs from the classical methods that combine the objective functions in a single scale, and is based on a unique run of non-linear single-objective optimizers.

INTRODUCTION

Many engineering design problems are multi-objective in nature as they often involve more than one design goal to be optimized. These design goals potentially impose conflicting requirements on the technical and economic performance of a given system design. To study the trade-offs between these conflicting design objectives and to explore design options, an optimization problem with multiple objectives has to be formulated. To facilitate the exploration of different multi-objective formulations, we can incorporate Pareto optimality concepts into optimization algorithms that require the designer's involvement as a decision maker. Edgeworth (1881) was the first academic to define an optimization problem involving multiple criteria. The problem, elaborated in the context of two consumers called P and π can be defined

A NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS

as: "it is required to find a point (x, y) such that, in whatever direction we take an infinitely step, P and π do not increase together, but, while one increases, the other decreases".

A few years later, in 1896, [START_REF] Pareto | Manual of Political Economy[END_REF] establishes the optimum for n consumers: "we will say that members of a collectivity enjoy maximum ophelimity in a certain position when it is impossible to find a way of moving from that position very slightly in such a manner that the ophelimity enjoyed by each of the individuals of that collectivity increases or decreases. That is to say, any small displacement in departing from that position necessarily has the effect of increasing the ophelimity which certain individuals enjoy, and decreasing that which other enjoy, being agreeable to some and disagreeable to others".

Since then, many researchers have been dedicated to developing methods to solve this kind of problem. Interestingly, solutions for problems with multiple objectives, also called multi-criteria optimization or vector optimization, are treated as Pareto optimal solutions or Pareto front, although as Stadler observed (1988), they should be treated as Edgeworth-Pareto solutions.

On the other hand, thanks to the evolution of computers, the optimization of large scale problems has become a common task in engineering design processes. Computers, ever faster, with higher capacity of data storage and Internet connection capability, are revolutionizing the traditional ways of developing engineering designs. Currently, engineers can rely on a wide range of design alternatives and optimization methods that allow systematic choices between alternatives when they are based on measurable criteria. Used properly, these methods can improve or even generate the final solution of a design, as shown in [START_REF] Giassi | Multi-Objective Optimization Using Asynchronous Distributed Applications[END_REF], which held multi-objective optimization of a vessel hull shape with collaborative design method where the computational tools for the analysis of alternatives used in the cycles of optimization were distributed in three different sites. With this approach they demonstrate that if the optimization is conducted in a distributed way it becomes a powerful tool for the concurrent design.

Rao et al. (2007) designed a disk for high-pressure turbine aircraft engines, treating the problem as multi-objective by considering the manufacturing cost and fatigue life as functions to be optimized. [START_REF] Bouyer | The Multiobjective Optimization of a Prismatic Drive[END_REF] performed the optimal design of a mechanical transmission system, Slide-O-Cam, which turns rotary into translation motion and replaces the conventional rack-pinion systems, a problem with four objective functions. In both instances, the optimal solutions were presented by their Pareto fronts.

Plenty of methods for solving design problems with more than one objective are found in scientific and technical literature. However, since the multi-objective problems have no single solution but a set of points called optimal solutions, the question that arises is how to choose one alternative among those contained in the set of solutions. The answer is not simple unless another criterion is added to help choose the alternative that should be deployed.

In this paper we propose a different alternative to solve the problem of optimizing a set of objective functions that will be divided into two groups: one containing only one function that will be called control function and one in which the p functions will be called the performance functions group. With this strategy, the optimization problem is formulated as single-objective, optimizing the control function over the Pareto set that would result from the optimization problem constituted by the group of performance functions if this multi-objective optimization problem had been solved. However, since this task is not performed, the Pareto front is not an output of the optimization process and the computational task is significantly smaller than the traditional multi-objective optimization algorithms.

Applying this strategy to a multi-objective optimization problem requires the use of an algorithm that can solve a single-objective problem and work with the functions involved. The result will be a single solution. However, the methodology differs significantly from other classical algorithms as there is no need to know a value function relating the objectives that articulates de designer preferences.

After this brief introduction, the remaining of the paper is organized into six sections. The first formulates a general multi-objective optimization problem and defines the nature of optimal solutions from the Pareto perspective and the necessary conditions to be met. The second presents a brief literature review on multi-objective optimization problem algorithms. In the sequence, it presents the main contribution of this work, a new strategy to solve the optimization problem with multiple objectives. Next, it goes on to the implementation of the strategy to solve three problems with increasing levels of complexity, and finally presents the conclusions on the work done. average of the i th decision variable included in the optimal set X decision variable vector X * non-dominated solution of a multi-objective optimization problem X inf , X sup lower and upper bounds of the design space X extended extended vector for unknowns with X, α, λ and μ α i weighting factor for i th objective function gradient in KKT condition α vector of α is ω i weighting factor for the i th objective function ω * i normalizing weighting factor for the i th objective function ε j upper bound for the j th objective function in the ε-restricted method λ j weighting factor for j th inequality constraint gradient in KKT condition λ vector of λ js σ i standard deviation for the i th decision variable included in the optimal set μ i weighting factor for i th equality constraint gradient in KKT condition μ vector of μ is ∇ gradient operator

Nomenclature

MULTI-OBJECTIVE OPTIMIZATION PROBLEM

Multi-objective optimization problems (MOOP) can be defined by the following equations:

minimize: f(X) (1a)
subject to:

g i (X) ≤ 0, i = 1, 2, . . . m. (1b) h j (X) = 0, j = 1, 2, . . . , l. (1c) X inf ≤ X ≤ X sup (1d)
where

f(X) = f 1 , f 2 , f 3 , . . . , f k T : X → R k
is a vector with the values of objective functions to be minimized. X is the vector containing the design variables, also called decision variables, defined in the design space R n . X inf and

Pesquisa Operacional

X sup are respectively the lower and upper bounds of the design variables. g i (X) represents the i th inequality constraint function and h j (X) the j th equality constraint function. The three equations (1b)-(1d), define the region of feasible solutions, S, in the design space R n . The constraints g i (X) are of type "less than or equal" functions in view of the fact that "greater or equal" functions may be converted to the first type if they are multiplied by -1. Similarly, the problem is the "minimization" of the functions f i (X), given that functions "maximization" can be transformed into the former by multiplying them by -1.

Pareto optimal

The notion of "optimum" in solving multi-objective optimization problems is known as "Pareto optimal." A solution is said to be Pareto optimal if there is no way to improve one objective without worsening at least one other, i.e., the feasible point X * ∈ S is Pareto optimal if there is no other feasible point X ∈ S such that ∀i, j with i = j, f i (X) = f i (X * ) with strict inequality in at least one condition, f j (X) < f j (X * ).

Due to the conflicting nature of the objective functions, the Pareto optimal solutions are usually scattered in the region S, a consequence of not being able to minimize all the objective functions simultaneously. In solving the optimization problem we obtain the Pareto set or the Pareto optimal solutions defined in the design space, and the Pareto front, an image of the objective functions, in the criterion space, calculated over the set of optimal solutions.

Necessary condition for Pareto optimality

In fact, optimizing multi-objective problems expressed by equations Eqs. (1a)-(1d) is of general character. The equations represent the problem of single-objective optimization when k = 1.

According to [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF], such as in single-objective optimization problems, the solution X * ∈ S for the Pareto optimality must satisfy the Karush-Kuhn-Tucker (KKT) condition, expressed as:

k i=1 α i ∇ f i (X * ) + m j=1 λ j ∇g j (X * ) + l i=1 μ i ∇h i (X * ) = 0 (2a) λ j g j (X * ) = 0 (2b) λ j ≥ 0 (2c) μ i ≥ 0 (2d) α i ≥ 0; k i=1 α i = 1 (2e)
where α i is the weighting factor for the gradient of the i th objective function, calculated at the point X * , ∇ f i (X * ). λ j represents the weighting factor for the gradient of the j th inequality constraint function, ∇g j (X * ), and is zero when the constraint function associated is not active, i.e., g j (X * ) > 0. μ i represents the weighting factor for the gradient of the i th equality constraint function, ∇h i (X * ).

The set of Eqs. (2a) to (2e) form the necessary conditions for X * to be a Pareto optimal. As described by [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF], it is sufficient for the complete mapping of the Pareto front if the problem is convex and the objective functions are continuously differentiable in the S space.

Otherwise, the solution will depend on additional conditions, as shown by Marler and Aurora (2004).

LITERATURE REVIEW IN METHODS FOR SOLVING MOOP

Some researchers have attempted to classify methods for solving MOOP according various considerations. Hwang & Masud (1979) and later [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF] suggested the following four classes, depending on how the decision maker (DM) articulates preferences: no-preference methods, a priori methods, a posteriori methods, and interactive methods.

In no-preference articulation methods, the preferences of the DM are not taken into consideration.

The problem can be solved by a simple method and the solution obtained is presented to the DM which will accept or reject it.

In a priori preference articulation methods, the hopes and opinions of the DM are taken into consideration before the solution process. Those methods require that the DM knows beforehand the priority of each objective transforming the multi-objective problem in a single-objective problem where the function to be optimized is a combination of objective functions.

In posteriori preference articulation methods no preferences of the DM are considered. After the Pareto set has been generated, the DM chooses a solution from this set of alternatives.

In interactive preference articulation methods the DM preferences are continuously used during the search process and are adjusted as the search continues.

Before explore the basics of some MOOP algorithms, it is opportune to define the ideal point

f o = f o 1 , f o 2 , f o 3 , . . . , f o k T ∈ R k being f o i the minimum of f i (X), X ∈ S and i = 1, 2, . . . , k.
In general, f o is unattainable, i.e., such a point in the criterion space does not map to a point in the design space.

In many cases it is advantageous to transform the original objective functions. This is especially true with scalarization methods that involve a priori articulation of preferences. As the objective functions usually have different scales, it is necessary to transform them in a non-dimensional. The relation

f s i (X) = f i (X) f o i (3)
can be used provided that f o i = 0. Such transformation results in a non-dimensional objective function with a lower limit of one and unbounded upper limit.

Pesquisa Operacional

OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO Alternatively, a most robust approach would be

f s i (X) = f i (X) -f o i f max i -f 0 i (4)
with f max i , the maximum value of f i (X) for X ∈ S. This operation is referred as normalization as the transformed function is bounded into the zero-one interval.

The weighted sum method

One of the most intuitive ways used to obtain a single unique solution for multi-objective optimization is the weighted sum method. In this approach, the MOOP are converted into a scalar preference function using a linear weighted sum function of the form, minimize:

k i=1 ω i f s i (X) (5a)
subject to: X ∈ S (5b)

ω i ≥ 0; k i=1 ω i = 1 (5c)
For this method, the weights ω i reflect, a priori, the designer's preferences. It is simple, but the proper selection of the weights may be a challenge itself.

The method can be used to find the Pareto set in many individual runs if the weights are consistently chosen for each run. Nevertheless, varying the weights continuously may not necessarily result in an even distribution of Pareto optimal points and a complete representation of the Pareto optimal set for non-convex problems as shown in Miettinen (1998).

The weighted metric methods

Other means of combining multiple objectives into a single-objective are based on the weighted distance metrics. Some DMs aim to find a feasible design that minimizes its distance from a pre-defined design as the representation of a designer's overall preferences. Assuming that the pre-defined design is the ideal design, the weighted p-norm as is expressed by

l p = k i=1 ω * i | f i -f o i | p 1 p (6)
where,

ω * i = ω i f max i -f o i is the normalizing weight. Pesquisa Operacional A NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS
The MOOP is written as minimize:

l p (7a) subject to: X ∈ S (7b)

ω i ≥ 0; k i=1 ω i = 1 (7c)
The parameter p may be chosen from 1 to infinite. With p = 2, Eq. ( 6) yields to the Euclidian distance metric.

A compromise solution, an alternative to the idea of Pareto optimality is a single point that minimizes the Euclidian distance between the potential optimal point and the ideal point.

The min-max method

The min-max solution was initially proposed by [START_REF] Lightner Mr & Director | Multiple criterion optimisation for the design of electronic circuits[END_REF]. The method try to find a feasible design that minimizes its distance from the ideal design as the representation of a designer's overall preferences. The weighted ∞-norm is used as distance metric.

In Eq. ( 6), the limit of l p as p approaches to ∞ is

l ∞ = l p→∞ yields -→ max i ω * i | f i -f o i | (8)
because the largest ω * i ( f i -f o i ) will dominate all others when taken to the infinite power. The min-max problem is then defined by minimize:

max i ω * i | f i -f o i | (9a) subject to: X ∈ S (9b)
Given different relative weights of objectives, the min-max method is capable of discovering efficient solutions of a multi-objective problem whether the problem is convex or non-convex.

The goal programming

Goal programming, originally proposed by [START_REF] Cooper Ww | Goal programming and multiple objective optimisation -Part 1[END_REF], is a technique which requires preference information before any efficient solution are generated. In fact, the method requires a designer to set goals, b i for all objectives, f i , that he wishes to achieve. It adopts the decision rule that the best compromise design should be the one which minimizes the weighted sum of deviations from the set goals, been reached. Similarly, the j th -inequality or equality constraint function are treated like the objective functions setting the goals b j = 0. Then, the MOOP is formulated as: minimize:

k i=1 ω * i (d + i + d - i ) (10a)
subject to:

f i (X) + d + i + d - i = b i , i = 1, 2, . . . , k (10b) d + i ≥ 0, d - i ≥ 0 and d + i d - i = 0, i = 1, 2, . . . , k (10c)
ω i ≥ 0 (10d) X ∈ S (10e)
The method allows the designer to assign preemptive weights to objectives and to define different achievement levels of the goals.

The ε-constraint method

Haimes et al. (1971) introduced the ε-constraint strategy that minimizes the single objective function f i (X). All other objective functions are used to form additional constraints. The MOOP is formulated as: minimize:

f i (X) (11a)
subject to:

f j (X) ≤ ε j , j = 1, 2, . . . , k; ∀ j = i (11b) X ∈ S (11c)
The definition of the limits ε j requires knowing a priori the designer's preference. A set of Pareto optimal solutions can be obtained with a systematic variation of ε j . However, improper selection of ε j ∈ R can result in a formulation with no feasible solution.

Nature inspired metaheuristic algorithms

The methods for multi-objective optimization presented thus far have involved unique formulations that are solved using standard optimization engines or single-objective optimization methods algorithms. With those methods, only one Pareto optimal solution can be expected to be found in one simulation run of a classical algorithm and not all Pareto optimal solution can be found by some algorithms in non-convex MOOP.

However, other approaches such some heuristics inspired in nature process can solve MOOP getting the Pareto set directly.

Genetic algorithms

Genetic algorithm (GA) is a metaheuristic inspired by the Darwin evolutionist theory explaining the origin of species. In nature, weak and unfit species within their environment are faced with

Pesquisa Operacional A NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS
extinction by natural selection. The strong ones have greater opportunity to pass their genes to future generations via reproduction. In the long run, species carrying the correct combination in their genes become dominant in their population. Sometimes, during the slow process of evolution, random changes may occur in genes. If these changes provide additional advantages in the challenge for survival, new species evolve from the old ones. Unsuccessful changes are eliminated by natural selection.

The concept of GA was originally proposed by Holland (1974) for applications into the control theories and it was quickly generalized to many different areas of engineering and sciences. The specific mechanics of the algorithm involve the language of microbiology and, in developing new potential solutions, mimic genetic operations. A population represents a group of potential solution points. A generation represents an algorithmic iteration. A chromosome is comparable to a design point, and a gene is comparable to a component of the design vector. Given a population of designs, three basic operations are applied: selection, crossover, and mutation. The selection operator involves selecting design vectors, called parents, in the current generation, which are combined together, by crossover, to form new chromosomes, called offspring. By iteratively applying the crossover operator, genes of good chromosomes are expected to appear more frequently in the population, eventually leading to convergence to an overall good solution. The mutation operator introduces random changes into characteristics of chromosomes. Mutation reintroduces genetic diversity back into the population and assists the search escape from local optima.

Being a population-based approach, GA is well suited to solve multi-objective optimization problems finding a set of multiple non-dominated solutions in a single run. The ability of GA to simultaneously search different regions of a solution space makes it possible to find a diverse set of solutions for difficult problems with non-convex, discontinuous, and multi-modal solutions spaces. In addition, most of multi-objective algorithms do not require the user to prioritize, scale, or weigh objectives. Consequently, GA has been the most popular heuristic approach to multi-objective design and optimization problems. Coello (2010) maintains an updated list with more than 5000 titles of publications involving different genetic algorithms.

The first multi-objective genetic algorithm, called vector evaluated genetic algorithm, or VEGA, was proposed by [START_REF] Schaffer | Multiple Objective Optimization with Vector Evaluated Genetic Algorithms[END_REF]. Afterwards, several multi-objective evolutionary algorithms were developed, as shown by [START_REF] Konak | Multi-objective optimization using genetic algorithms: A tutorial[END_REF], including the one that is used in this work to compare problems' results, the non-dominated sorted genetic algorithm, or NSGA II, proposed by Deb et al. (2000).

In the NSGA II, before selection is performed, the population is ranked on the basis of nondomination: all non-dominated individuals are classified into one category to provide an equal reproductive potential for these individuals. Since individuals in the first front have the maximum fitness value, they always get more copies than the rest of the population when the selection operator is applied. Additionally, the NSGA II estimates the density of solutions surrounding a particular solution in the population by computing the average distance of two points on either side of this point along each of the objectives of the problem. This value is called crowding distance. During selection, the NSGA II uses a crowded-comparison operator which takes into consideration both the non-domination rank of an individual in the population and its crowding distance. The non-dominated solutions are preferred over dominated solutions, but between two solutions with the same non-domination rank, the one that resides in the less crowded region is preferred. The NSGA II uses the elitist mechanism that consists of combining the best parents with the best offspring obtained. [START_REF] Goldberg | Genetic Algorithms in Search and Machine Learning[END_REF] and Deb (2001) are excellent guides for GA implementations.

Particle swarm

Particle swarm optimization algorithm (PS) was originally proposed by Kennedy & Eberhart (1995). It is a population-based search algorithm based on the simulation of the social behavior of birds within a flock. Although originally adopted for balancing weights in neural networks, PS soon became a very popular global optimizer.

There are two main distinctions between PS and GA. Genetic algorithms rely on three mechanisms in their processing. In contrast, PS only relies on two mechanisms, since PS does not adopt an explicit selection operator. The absence of a selection mechanism in PS is compensated by the use of leaders to guide the search. Such set of leaders is usually stored in a different place in an external archive, with all non-dominated solutions found so far. However, there is no notion of offspring generation in PS as with evolutionary algorithms. A second difference has to do with the way in which the individuals are manipulated. PS uses an operator that sets the velocity of a particle to a particular direction. This can be seen as a directional mutation operator in which the direction is defined by both the particle's personal best and the global best of the swarm. If the direction of the personal best is similar to the direction of the global best, the angle of potential directions will be small, whereas a larger angle will provide a larger range of exploration. In contrast, evolutionary algorithms use a mutation operator that can set an individual in any direction, although the relative probabilities for each direction may be different.

The main algorithm of PS is relatively simple, it only adopts one operator for creating new solutions and its implementation is, therefore, straightforward. Its main drawback is how to control de size of the external archive that implies in the computational time expended in the optimization process.

Simulated annealing

Simulated annealing optimization algorithm (SA) is based on an analogy of thermodynamics with the way metals cool and anneal. If a liquid metal is cooled slowly, its atoms form a pure crystal corresponding to the state of minimum energy for the metal. Serafini (1985Serafini ( , 1992)). The algorithm of the method is almost the same as the algorithm of single objective SA. The slow convergence rate for some optimization problems is the main drawback of SA, as mentioned by [START_REF] Suman | A survey of simulated annealing as a tool for singleand multiobjective optimization[END_REF].

Ant Colony optimization algorithm

The ant colony optimization algorithm (AC) is inspired by the behavior of ants and other insects that live in a colony that in spite of the simplicity of each individual, present a high level of social organization when observed together. Some examples of ant colony's capabilities found in Dorigo et al. (1999) are: division of labor and task allocation, cemetery organization and brood sorting, cooperative transport and finding the shortest path between two or more locations, often between a food source and a nest.

The first AC algorithm developed was initially applied to the traveling salesman problem, Dorigo (1992). The algorithm was based on the ant colony capability to find the shortest path between a food source and a nest. The algorithm uses artificial ants that cooperate on finding solutions to the problem through communication mediated by artificial pheromone trails.

While moving on the graph associated with the problem, artificial ants deposit pheromone on the edges traversed marking a path that may be followed by other members of the colony, which then reinforce the pheromone on that path. With this kind of communication, ants have their activities coordinated. This self-organizing behavior results in a self-reinforcing process that leads to the formation of a path marked by high pheromone concentration, while paths that are less used tend to have a diminishing pheromone level due to evaporation.

This concept can be applied to any combinatorial optimization problem for which a constructive heuristic can be defined. The process of constructing solutions can be regarded as a walk on a construction graph where each edge of the graph represent a possible step the ant can take. AC is essentially constructive, as ants generate solutions by adding solution components, corresponding to the edges chosen, to an initially empty solution until the solution is complete.

Artificial immune systems

The human immune system has as its main task the detection of the infectious foreign elements, called pathogens, that attack us, and defend us from them, i.e., its main task is to keep our organism healthy. Examples of such pathogens are bacteria and viruses. Any molecule that can be recognized by our immune system is called antigen. Such antigens provoke a specific response from our immune system. Lymphocytes of types B and T are special type of cells that play a major role in our immune system.

Upon detection of an antigen, the B cells that best recognize the antigen are cloned. Some of these cloned cells will be differentiated into plasma cells, which are the most active antibodies secretors, while others will act as memory cells. These cloned cells are subject to a high somatic mutation rate in order to increase their affinity level, i.e., their matching to the antigens. These

Pesquisa Operacional

OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO mutations experienced by the clones are proportional to their affinity to the antigen. The highest affinity cloned cells experiment the lowest mutation rates, whereas the lowest affinity cloned cells have high mutation rates. Due to the random nature of this mutation process, some clones could be dangerous to the body and are, therefore, eliminated by the immune system itself. Plasma cells are capable of secreting only one type of antibody, which is relatively specific for the antigen. Antibodies play a key role in the immune response, since they are capable of adhering to the antigens, in order to neutralize and eliminate them.

Once the antigens have been eliminated by the antibodies, the immune system must return to its normal conditions, eliminating the in-excess cells. However, some cells remain in our blood stream acting as memory cells, so that our immune system can 'remember' the antigens that have previously attacked it.

When the immune system is exposed again to the same type of antigen or a similar one, these memory cells are activated, presenting a faster and perhaps improved response, which is called secondary response.

Based on the previous explanation of the way in which human immune system works, it can be say that, from a computer science perspective, the immune system can be seen as a parallel and distributed adaptive system. Clearly, the immune system is able to learn, it has memory, and is able of tasks such as associative retrieval of information. These features make immune systems very robust, fault tolerant, dynamic and adaptive. All of these properties can be emulated in a computer.

Artificial immune systems (AIS) are composed of the following basic elements: a) a representation for the components of the system, e.g., binary strings, vectors of real numbers; b) a set of mechanisms to evaluate the interaction of individuals with their environment and with each other; such an environment is normally simulated through an affinity function, which is based on the objective functions in the case of optimization problems; c) procedures of adaptation, that indicate the way in which the behavior of the system changes over time; these procedure of adaptation consist of, for example, mutation operators.

The first direct use of the AIS for multi-objective optimization goes back to Yoo & Hajela (1999).

In their work they use a standard genetic algorithm where the immune principle of antibodyantigen affinity is employed to modify the fitness value. In a first time the population is evaluated versus the problem objectives and different scalar values are obtained by making reference to different weighting combinations. The best individual with respect to each combination is identified as antigen. The rest of the population is the pool of antibodies. Then antibodies are matched against antigens by the definition of a matching score. The best matching antibody fitness is added by this matching score, evolving the population of antibodies to cover antigens.

DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS

Usually, when solving a MOOP, the DM expects, in general, to obtain its Pareto front and its Pareto set. Although it is very interesting to know these non-dominated solutions, additional criteria are necessary to select a single solution that will be deployed. Even though useful to better understand the inter-relationships between the objectives, the Pareto front may become an obstacle since just one solution must be implemented. To make a selection, the DM will have to use an additional criterion, be it subjective or not.

Then, why not formalize and include this additional criterion into the problem and make it a control function to find the single and final solution to be deployed among those that belong to Pareto front?

Based on this idea, we propose a new methodology to formulate MOOP involving non-linear differentiable functions. The method transforms the MOOP from multi-objective into singleobjective that can be solved with the aid of any traditional single-objective optimization engine suitable for the problem in focus.

Proposal of a new methodology for decision making in MOOP

In the proposed method, the objective functions are divided into two groups: a) the control function group, or simply, control function, which contains only one function.

b) the performance functions group, which is made up of the functions that will provide the Pareto set.

Which function will be part of each group is a designer's choice and depends on his or her experience and knowledge related to the design problem. The control function can be an additional objective function or it can be elected among the problem objective functions.

The performance functions will have an important role in the process as they will provide the Pareto set as a constraint over which the control function will search the solution that optimizes it.

Thus, the multi-objective optimization problem can be written as: minimize:

f c (X), f p (X) (12a)
subject to:

g i (X) ≤ 0, i = 1, 2, . . . , m (12b) h j (X) = 0, j = 1, 2, . . . , l (12c) X inf ≤ X ≤ X sup (12d)
where f c (X) is the control function,

f p (X) = f 1 , f 2 , . . . , f p : X → R p
is the vector composed of the p objective functions in the performance functions group.

To apply the proposed methodology, the performance functions, is substituted by the KKT necessary condition, Eqs. (2a)-(2e), in such way that the problem's final solution belongs to the Pareto front of the MOOP with the performance functions only. are not known. As unknowns in the problem, they will be incorporated into the vector of design variables, defining the extended vector of unknowns:

X extended = (X, α, λ, μ) (13)
Finally, the problem is formulated as a single-objective optimization problem, with the control function, f c (X), to be minimized and constrained by the conditions for obtaining the Pareto optimal solutions considering only the performance functions, f 1 (X), f 2 (X), . . . , f p (X).

Mathematically, the optimization problem is formulated as:

find X extended that minimizes:

f c (X) (14a)
subject to:

Pareto set condition for the performance functions:

                                                     p i=1 α i ∇ f i (X) + m j=1 λ j ∇g j (X) + l i=1 μ i ∇h i (X) = 0 g i (X) ≤ 0, i = 1, 2, . . . , m h j (X) = 0, j = 1, 2, . . . , l λ j g j (X) = 0 λ j ≥ 0 μ i ≥ 0 α i ≥ 0; p i=1 α i = 1 X inf ≤ X ≤ X sup (14b) (14c) (14d) (14e) (14f) (14g) (14h) (14i)
The proposed methodology differs from those methods exposed in the items 3.1 to 3.4. Those strategies require that a complete decision-making structure of the problem is decided a priori. Although the proposed method is also based on a priori decisions, just one function shall be added or isolated as a control function. The remaining functions are incorporated into the formulation to ensure that the final solution is within the Pareto set obtained from the minimization problem of these functions, even though the Pareto set and the Pareto front is not an explicit outcome of the process.

However, it must be considered that, with the use of Eq. (2a), the methodology applies only to problems where the functions involved are continuously differentiable. Furthermore, the solution

Pesquisa Operacional A NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS
obtained is a local optimum and it can be a global optimum for convex Pareto fronts, since Eq. (2a) is a necessary condition for the existence of an extreme point and not a sufficient condition.

APPLICATIONS

To validate the proposed multi-objective optimization method, it will be used to solve three examples with increasing levels of complexity, namely the optimization problem of three quadratic functions, the design of a cantilever beam and the conceptual design of a bulk carrier.

To solve the single-objective optimization problem originated by the proposed methodology, any algorithm that works with optimization problems involving nonlinear functions and constraints can be used. Due to its accessibility, the solver fmincon was used. It is a component of the Optimization Toolbox available in the application MatLab, Matworks, version 7.1.

The NSGA II was used to compare results and to show that the results of the proposed method lay over the Pareto set from the MOOP of the performance functions group. A minor adaption was done in the NSGA II in order to interrupt the evolution if there is no significant difference between consecutive generations. Calling f i j the value of the i th objective function for the j th chromosome at n th generation of a population with pop chromosomes, the root mean square of the objective functions for the population can be defined as [START_REF] Rc | Particle swarm optimization[END_REF] with f i max and f i min being the maximum and the minimum value, respectively, of the i th objective function evaluated in the population. The alternative stop criterion was defined as

r ms n = 2 k i pop j f i j 2 f i max -f i min 2
|r ms n -r ms n-1 | ≤ tol f • r ms n-1 ( 16
)
where tol f is a small number. Sometimes, for a defined tol f value, this criterion may cause a premature interruption of the evolutionary process. To overcome this situation, the algorithm checks if it occurs in subsequent generations before discontinue the evolution. As a default, three subsequent generations were adopted with the satisfaction of Eq. ( 16) as the condition to interrupt the evolutionary process.

Minimization of three quadratic functions

Consider the minimization problem defined by: minimize:

f 1 (x 1 , x 2 ) = 3(x 1 + 5) 2 + (x 1 + 5)(x 2 -2) + (x 2 -2) 2 (17a) f 2 (x 1 , x 2 ) = (x 1 -5) 2 -(x 1 -5)(x 2 -3) + (x 2 -3) 2 (17b) f 3 (x 1 , x 2 ) = (x 1 + 2) 2 + (x 2 + 6) 2 (17c) Pesquisa Operacional
OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO Assume that the functions f 1 (x 1 , x 2 ) and f 2 (x 1 , x 2 ) integrate the performance functions group and the function f 3 (x 1 , x 2 ) is the control function. Accordingly, the problem can be formulated as: find

X extended = x 1 , x 2 , α 1 , α 2 T that minimize: f c (x 1 , x 2 ) = (x 1 + 2) 2 + (x 2 + 6) 2 (18a)
subject to:

α 1 ∇ f 1 (x 1 , x 2 ) + α 2 ∇ f 2 (x 1 , x 2 ) = 0 (18b) α i ≥ 0; 2 i=1 α i = 1 (18c)
The solution of this problem, (x

Figures (1b) and (1c)

show the Pareto set and the corresponding Pareto front for the MOOP involving the performance functions f 1 (x 1 , x 2 ) and f 2 (x 1 , x 2 ) obtained by the NSGA II. Although the points in the Pareto front seem evenly distributed, the corresponding points in the Pareto set configure a rough approximation of the solution. To get this set, a 50 chromosomes population was used with the alternative stop of evolution criteria, Eq. ( 16), tol f = 0.001. With these parameters the Pareto front stabilizes in 177 generations. The same parameter was used in the other NSGA II examples shown in Table 1. 1f), where the Pareto set is a region with infinite non-dominated points, the question that normally arises in the designer's mind is how to choose only one as a solution. In general, as more functions are aggregated to a MOOP, the more extended the Pareto set on the design space. In the example, with two functions, the Pareto set is a line. With three functions, the Pareto set is a plane surface.

Pesquisa Operacional

In addition, comparisons for the computer performance between the algorithms are shown in Table 1. As there are no other methodologies similar to the one proposed in this paper, the results were compared to two other algorithms: the weighted sum approach and the NSGA II. For both, the process involved two phases: a) the search of the non-dominated designs for the MOOP involving the performance functions group and b) the manual search in the resultant set for the alternative that minimizes the control function. The number of function calls resulted from the first part of the process. The computational performance measure adopted for all algorithms' comparison is the number of objective functions calls. Although this measure is not perfect since both algorithms -the genetic algorithm and weight sum approach -depend on the number of points to be used in the MOOP solver, it has at least a qualitative value. For sake of comparison the computational times expended to achieve the final solution in a 2 GHz dual processor computer with 3 Gb RAM are also shown.

Two other methods that return a single solution are shown in Table 1 and figure (1f), the compromise and the min-max solutions. They are very fast, but they may not reflect the best compromise preferred by the DM.

Although the example is very simple, it shows how the proposed methodology can help decision making processes. The problem solution falls within the performance functions non-dominated solutions and it is the one that optimizes the control function. In the quadratic functions example, the control function was chosen among the problem objective functions. In the next application, the MOOP is defined as minimizing the mass and tip deflection of a cantilever beam, two antagonist functions. With the proposed methodology these functions will compose the performance group. Afterward, it will be included a third function, the manufacturing cost, that will play the role of the control function.

Pesquisa Operacional A NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS

Design of a cantilever beam

As a second illustrative example, consider the design problem of a cantilever beam, adapted from Deb (2001). The beam, prismatic and with a circular section, must support the weight P at its free end. The design variables are d and L, the cross section diameter and beam length, respectively. Also consider two attributes to be minimized, the beam mass and the beam tip deflection when subjected to the weight P. Moreover, the maximum bending stress must be below of the material yield stress, σ y , and the maximum tip displacement limited to δ max . A schematic representation of the problem is shown in Figure 2. It can be demonstrated that the two selected objectives conflict since the minimization of the mass will lead to lower values of the pair (d, L) and minimization of the tip deflection will result in an increase in the cross sectional diameter with a reduction in the beam length. relationship is observed considering the cross section diameter. Fixing the beam length, the lower the diameter the higher the tip deflection and the lower the beam mass.

Similar results were obtained by the NSGA II algorithm, using a 50 chromosomes population which has evolved 58 generations, conditioned by the stop criteria, Eq. ( 16), with tol f = 0.001. If the cost function is included as the third function in the MOOP, there is a growth and dispersion of Pareto optimal solutions in the design space, as shown in Figure (5b). This spread of the nondominated alternatives in the design space makes decision making process even more difficult.

The proposed methodology can help the DM elect one alternative among all non-dominated, while preserving the concept of multi-objective optimization. The MOOP of the design of a cantilever beam can be re-formulated as:

Among all non-dominated design alternatives for the cantilever beam, where it is searched to minimize the beam tip deflection, δ(d, L), and minimize the beam mass, m(d, L), find the one that has the minimum manufacturing cost.

Another interpretation for problem formulation would be: among the that satisfy a or engineering criterion, choose the one that has the lowest manufacturing cost. Applying the proposed methodology, the MOOP can be rewritten as:

Pesquisa Operacional

find X extended = x 1 , x 2 , α 1 , α 2 , λ 1 • • • λ 6 T that minimizes: $(d, L) = 100 d -30 50 2 + d -30 50 L -400 1000 + L -400 1000 2 (21a) subject to: α 1 ∇δ + α 2 ∇m + 6 j=1 λ j ∇g j = 0 (21b) g j ≤ 0, j = 1, 2, . . . , 6 (21c)
λ j g j = 0, j = 1, 2, . . . , 6 (21d) Pesquisa Operacional
α i ≥ 0; α 1 + α 2 = 1 (21f) g 1 = σ (d, L) -σ y (21g) g 2 = δ(d, L) -δ max (21h) g 3 = 10 -d (21i) g 4 = d -50 (21j) g 5 = 200 -L (21k) g 6 = L -1000 (21l)
The fmincon was used to solve this single-objective optimization problem, which returns the following values Table 2 shows a comparison of the methodologies used to find the single solution for the beam design. Although the proposed methodology does not return the Pareto front it is very efficient to find the final solution with number of function calls significantly smaller than the other two methods used for the same purpose. To conclude, if the cost function is treated globally as the third function in MOOP, the set of Pareto optimal solutions will no longer be a line, as shown in Figure (4a); instead it is spread Pesquisa Operacional OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO over a surface, as shown in Figure (5b). In fact, the greater the number of conflicting functions in a multi-objective optimization problem, the greater the region in the space of the design variables that defines the Pareto optimal solutions. Consequently, the decision making process to elect a single solution is more difficult. How do we choose a non-dominated alternative among those present in Figure (5b)? The answer is impossible unless a fourth criterion is used to drive the choice. Rather, in the present proposal, this additional criterion is added to the MOOP as a control function to guide a convenient solution that satisfies all the constraints in a very efficient way, a solution that is non-dominated inside the group of the performance functions, and is an appropriate optimal solution obtained by optimizing the control function.

Conceptual design of a bulk carrier

The third application of the developed methodology was in the conceptual design of a bulk carrier. The conceptual design of a cargo vessel is not a trivial task. For decades, this problem has been handled in two ways, either by the adapting a known design aiming to the new requirements, or by the aid of simplified mathematical models driven by an optimization algorithm that enables obtaining the optimal solution based on technical or economic criteria previously established. This work explores the second alternative with the help of the mathematical model for conceptual design of bulk carriers developed by Pratyush & Yang (1998) and presented in Appendix A. The model is made up of functions that define the vessel attributes from which are drawn those that constrain the design space, those to be optimized and those that characterize the technical and economic performance of the vessel and allow the evaluation of each design alternative. Among them are the annual transportation cost, the annual transported cargo, the ship lightweight, the ship initial cost and other functions of the vessel design variables such as length, width, depth, draft, block coefficient and speed, respectively (L , B, D, h, C B , V k ). Pratyush & Yang (1998) proposed a general problem with explicit bounds on the block coefficient that should range between 0.63 and 0.75, and the ship speed that varies between 14 and 18 knots, and nine constraints in attributes, presented in Table A1 of Appendix A.

Parsons & Scott (2004) used the same model for comparing optimization methods, but they have altered the explicit bounds on the design variables. They limited the ship length at 274.32 m (900 ft) and the minimum deadweight at 25,000 ton, as shown in Table A2 of Appendix A. They also evaluated the problem by constraining the decision space in order to design ships suitable to cross the Panama Canal, with beam and draft limited to 32.31 m and 11.71 m, respectively, as shown in Table A3 of Appendix A.

Using Parsons and Scott's second model, Hart & Vlahopoulou (2009) constrained the depth to the minimum and the maximum of 13 m and 25 m, respectively, and extended the acceptable range for ship speed to the values of 11 knots and 20 knots, minimum and maximum respectively, as shown in Table A4 of Appendix A.

The bounds for the constraints feasibility lead to different results in optimization and the constraint bounds imposed in Pratyush & Yang's (1998) original work may be associated with the mathematical model, based on data reductions from operating ships and whose extrapolation Pesquisa Operacional may undermine not only the quality but even the validity of the results obtained by using the model. In this paper was used the constraint bounds on design variables proposed in Pratyush & Yang's (1998) original work. For practical reasons, the limits they have not mentioned were adopted wide enough to not influence the optimization results. Accordingly, the ranges for the design variables are described in Table A5 of Appendix A.

Conflicting objectives

Pratyush & Yang (1998) chose to minimize the annual transportation cost, maximize the amount of annual transported cargo and minimize vessel lightweight. In this work, we replaced the vessel lightweight by the ship initial cost. Although vessel lightweight and ship initial cost functions are related, the latter has a financial appeal like the former two objective functions. Annual transported cargo is associated with the annual income, the annual transportation cost with the annual expenses and the ship cost is related to the capital required for the ship purchase.

The design alternative that optimizes the single-objective problem can be easily obtained by taking each objective function separately. Any algorithm running with non-linear optimization problems and constraints can be used to solve this single-objective optimization problem. The fmincon was used to get the results shown in Table 3. The results indicate that the objective functions are conflicting as to maximize the annual cargo (A C ); the optimum design will be a large ship with the highest speed of 18 knots. The search reaches the upper bound for the total deadweight, 500,000 t. The opposite occurs when it minimizes the ship cost (C S ). In this case, all vessel dimensions decrease along with the ship speed, which reaches the minimum of 14 knots with a slender shape vessel with block coefficient of 0.63, assuring a low power main engine and consequently a lower main engine cost. The search for even smaller vessels is constrained by the lower bound for the total deadweight of 3,000 tons. To minimize the annual transportation cost (C T ), the solution is situated between the previous two.

Solutions for the MOOP

A clearer view of the antagonism of the objective functions in this problem can be obtained using an algorithm to find the Pareto front of multi-objective optimization problems. Figure 6 shows 10,000 feasible ship designs and highlights the non-dominated ones. To get these points, more than 12 million random design alternatives were generated and analyzed, taking 1097 seconds of computational time. Table 4 shows the design variables statistics for the feasible points. Although the problem shows a narrow surface in the criterion space, the design variables are well distributed over the design space.

With the aid of Figure 6, is possible to presume how the Pareto front of this problem should be, but, on the other hand, the random walk is not adequate to drive the search for the Pareto front, while the number of function calls becomes prohibitive.

In order to compare the results of the proposed methodology with other methods, the weighted sum approach and the genetic algorithm were used to do the bulk carrier design optimization task.

Figure 7 and Table 5 show the results of the weighted sum approach algorithm with 100 sets of random weights. As fmincon fails in find the solution for several weights sets, the process were restarted from 5 different starting points for each set single run, being chosen the best result as the resultant point. Even though points were discarded, because they have shown to be dominated solutions.

Better results are obtained from the NSGA II algorithm according to the points shown in Figure 8. These results were achieved for a population of 100 chromosomes, which has evolved 2000 generations. Table 6 shows the statistics for the non-dominated designs. The statistics for the design variables indicate a reasonable dispersion of the values. With the set of non-dominated solutions, shown in Figure 8, the question that often arises for the DM is: which design alternative should be implemented?

Pesquisa Operacional

Applying the proposed methodology

Applying the proposed methodology one attribute function can be chosen as the control function and optimized in the Pareto front generated by the optimization of the performance functions group. Consider the previous functions, the annual transportation cost, the annual cargo and the ship cost, incorporating the performance group and the voyage cost as the control function. Accordingly, the multi-objective optimization problem can be defined as: among all the non-dominated solutions that would be obtained by resolving the MOOP for the bulk carrier design, namely, minimizing the annual transportation cost (C T ), minimizing the ship cost (C S ) and maximizing the annual cargo (C A ), find the solution that minimizes the voyage cost (C V ).

The mathematical formulation of this optimization problem is shown in detail in Appendix B. Note that the problem still is multi-objective, with four objective functions and no value function involved. However, with the use of the proposed approach, it can be formulated and solved in a single passage of any algorithm dedicated to single-objective optimization with inequality and equality constraints and able to work with non-linear functions. In this manner formulated, the problem will have 30 unknowns (six design variables of the original problem, plus three weighting factors α k related to the gradients of the performance functions and more 21 weighting factors λ j related to the explicit and implicit constraints).

Table 7 shows the design alternative with some ship attributes obtained by using the fmincon algorithm. Additionally, the solution is highlighted in Figure 9 over the Pareto front of performance functions obtained thanks to the NSGA II.

Table 8 compares the results obtained by the proposed methodology with those obtained by the weighted sum approach and the NSGA II. The proposed methodology is very efficient in finding the optimized design as the number the function calls is significantly smaller.

It can be observed that the use of the proposed methodology leads to a solution that satisfies all the constraints and provides a good technical performance and its result is optimal according to the control function.

CONCLUSION

Most engineering design problems are multi-objective and the cases where the objective functions do not conflict are rare. To solve these kinds of problems, many researchers developed methods to search for the solution of multi-objective optimization without simplifying the problem to single-objective and having to decide a priori how to group the objective functions into a single scale. The evolutionary methods are frequently used to locate the set of solutions of multi-objective problems. These algorithms provide a discrete picture of the Pareto front in the criterion space. It was observed that the greater the number of objective functions, the more scattered the set of non-dominated solutions is in the decision space and the harder it is for the DM to choose an alternative to be deployed. This paper proposes a new methodology to solve multi-objective optimization problems in which one objective function is proposed or isolated and is treated as a control function that will drive the decision-making process. The other objective functions will form the performance functions group. Through this strategy, a single-objective optimization problem is formulated, in which the control function will be optimized over the Pareto set that would result from the optimization problem established by the performance functions, if this problem was previously solved.

Then the resultant single-objective optimization problem can be solved by any classical standard optimization engines, the limitations will be those present in the method used, such as the nature of the functions, whether they are linear or nonlinear, or continuous, regarding the presence or not of inequality and equality constraints.

With the proposed method a minimum knowledge is needed a priori. There is no need to know a value function relation to the objectives before starting to solve the problem as other a priori articulation of preferences methods.

Pesquisa Operacional

OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO As drawbacks, the method doesn't provide a set of non-dominated solutions and as in many other single-objective optimization algorithms, the solution converges to a local optimum since the conditions of Karush-Kuhn-Tucker, chief pillar of the proposed methodology, are the necessary conditions but not sufficient to ensure that the solution is in the "global" Pareto front. If the problem is convex and the functions involved are regular, the local solution also will be a global solution, ensuring a unique result for the problem. Another disadvantage, the functions involved in the problem shall be continuously differentiable.

Although these conditions are limiting, the proposed methodology is very efficient in solving engineering design problems, as demonstrated by the examples solved with its use.

APPENDIX B: BULK CARRIERS DESIGN MODEL WITH THE PROPOSED METHODOLOGY

To apply the proposed methodology is convenient to define the following entities:

a) The control function:

A

  NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS (a) Solution by the proposed methodology when f 3 is the control function. (b) Pareto set for the performance functions f 1 and f 2 obtained by the NSGA II algorithm. (c) Pareto front for performance functions f 1 and f 2 obtained by the NSGA II algorithm. (d) Pareto set for the performance functions f 1 and f 2 obtained by the weighted sum approach. (e) Solutions applying the metodology with different choices of the control function. (f) Pareto set for the three quadratic functions MOOP obtained by the weighted sum approach.

Figure 1 -

 1 Figure 1 -Results for the quadratic functions MOOP.

Figure

  Figure (1d) shows the Pareto set for the same problem obtained by using the weighted sum approach, with 50 weight vectors chosen at random. This approach returns a well-matched Pareto set. Although the NSGA II is a powerful algorithm when solving MOOPs, it focuses on the criterion space and consequently can generate rough solutions in the design space as show the results in Figures (1b) and (1d).

Figure

  Figure (1e) shows the results when another function is chosen to play the role of the control function. The resultant points are distinct but they are over the corresponding Pareto set resulting from the MOOP of the corresponding performance function group. Moreover, they belong to the Pareto set that results from the MOOP composed of the three quadratic functions, shown in Figure (1f). This Pareto set was obtained by the weighted sum approach with 1,200 sets of weights chosen at random. Observing Figure (1f), where the Pareto set is a region with infinite non-dominated points, the question that normally arises in the designer's mind is how to choose only one as a solution. In general, as more functions are aggregated to a MOOP, the more extended the Pareto set on the design space. In the example, with two functions, the Pareto set is a line. With three functions, the Pareto set is a plane surface.

Figure 2 -

 2 Figure 2 -Cantilever beam.

  Figure (4a) shows that all alternatives for the beam dimensions have length of 200 mm and cross section diameter ranging from 19 mm to 50 mm, approximately. These results are expected because, for a fixed diameter, to reduce the mass and the tip deflection the length decreases and the lower bound for this design variable is 200 mm, as shown in Figure (3a), but an inverse

Figure 3 -

 3 Figure 3 -Cantilever beam attributes.

  For comparison, Figures (4e) and (4f) show the results found by the weighted sum approach. Fifty random sets of weights were used to reduce the normalized mass and deflection beam functions to a single scale. Now suppose that the problem must include the beam manufacturing cost and that it should be minimized. Additionally, suppose the cost function being defined by the function: $(d, L) = 100 d -The contours of the cost function are shown in Figure (5a).

A

  NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS (a) Feasible and non-dominated designs for the cantilever beam in the design space. (b) Feasible and non-dominated designs for the cantilever beam in the criterion space. (c) Pareto set of the MOOP involving the beam mass and tip deflection obtained by the NSGA II algorithm. (d) Pareto front of the MOOP involving the beam mass and tip deflection obtained by the NSGA II algorithm. (e) Pareto set of the MOOP involving the beam mass and tip deflection obtained by the weighted sum approach. (f) Pareto front of the MOOP involving the beam mass and tip deflection obtained by the weighted sum approach.

Figure 4 -

 4 Figure 4 -Cantilever beam Pareto set and Pareto front found by distinct algorithms.

Figure 5 -

 5 Figure 5 -Adding the cost function to the MOOP involving the design of a cantilever beam.

A

  NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS λ j ≥ 0 (21e)

T

  (d, L) = 35.0077, 200.0000 α = 0.8144 0.1856 T λ = 0.0000 0.0000 0.0000 0.0000 0.0035 0.0000 Note that the design variables (d, L) and the weighting factors (only λ 5 is not zero, revealing that only the constraint equation g 5 , Eq. (21k), is active) converged to a point on the Pareto optimal set for the MOOP involving the performance functions only, as shown in Figure(4a). For the sake of illustration, the values of the manufacturing cost calculated over the Pareto set obtained for the bi-objective optimization problem are shown in Figure(5d). The minimum cost is located in the diameter axis at the neighborhood of 35 mm.

A

  NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS (a) Feasible and non-dominated designs for the bulk carrier obtained by a random walk in the design space. (b) x-y view of the feasible and non-dominated designs. (c) x-z view of the feasible and non-dominated designs. (d) y-z view of the feasible and non-dominated designs.

Figure 6 -

 6 Figure 6 -Feasible and non-dominated designs for the MOOP with the annual transportation cost (C T ) minimization, the ship cost (C S ) minimization and the annual cargo (C A ) maximization, obtained by a random walk over the design variables.

  OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO (a) Pareto front obtained by the weighted sum approach. (b) x-y view of the Pareto front. (c) x-z view of the Pareto front. (d) y-z view of the Pareto front.

Figure 7 -

 7 Figure 7 -Pareto front for the MOOP with the annual transportation cost (C T ) minimization, the ship cost (C S ) minimization and the annual cargo (C A ) maximization obtained by the weighted sum approach with 100 random weight sets.

A

  NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS (a) Pareto front for the bulk carrier design obtained by the NSGA II algorithm. (b) x-y view of the Pareto front. (c) x-z view of the Pareto front. (d) y-z view of the Pareto front.

Figure 8 -

 8 Figure 8 -Pareto front for the MOOP with the annual transportation cost (C T ) minimization, the ship cost (C S ) minimization and the annual cargo (C A ) maximization obtained by the NSGA II algorithm.

A

  NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS (a) Solution for the bulk carrier design obtained by the proposed methodology. (b) x-y view of Pareto front of the bulk carrier design obtained by the NSGA II. (c) Pareto front of the bulk carrier design obtained by the NSGA II. (d) Pareto front the bulk carrier design obtained the NSGA II.

Figure 9 -

 9 Figure 9 -Pareto front by the NSGA II for the MOOP involving the performance group {annual transportation cost (C T ), ship (C S ), and annual cargo (C A )} and solution obtained by the proposed methodology with voyage cost (C V ) as the control function.

9 T g 1 = 6 - 32 g 5 = 15 T g 10 = 21 TT

 916325151021 voyage cost: C v b) The performance functions: f p = C S , C T , A C T ship cost: C S transportation cost: C T annual cargo: A C c) The implicit constraints:g c = g 1 , g 2 , . . . , g L/B g 6 = 0.07B -K B -B M T -K G g 2 = L/D -15 g 7 = dwt -500000 g 3 = L/T -19 g 8 = 3000 -dwt g 4 = T -0.45dwt 0.31 g 9 = F n -0.T -0.7D + 0.7 d)The explicit lower bound constraints:g inf = X inf -X =g 10 , g 11 , . . . , g 60 -L g 13 = 3 -T g 11 = 10 -B g 14 = 0.63 -C B g 12 = 4 -D g 15 = 14 -V k e) The explicit upper bound constraints: g sup = X -X sup = g 16 , g 17 , . . . , g 21 T g 16 = L -600 g 19 = T -30 g 17 = B -100 g 20 = C B -0.75g 18 = D -40 g 21 = V k -18f) The (unknown) design variables:X = L , B, D, T, C B , V k Tg) The (unknown) performance functions gradient multipliers:α = α 1 , α 2 , α 3 T h)The (unknown) implicit constraints gradient multipliers:λ c = λ 1 , λ 2 , . . . , λ 9 Ti) The (unknown) explicit lower bound constraints gradient multipliers:λ inf = λ 10 ,λ 11 , . . . , λ 15 T Pesquisa Operacional OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO j) The (unknown) explicit upper bound constraints gradient multipliers:λ sup = λ 16 , λ 17 , . . . , λLet the Karush-Kuhn-Tucker condition be defined by:F(X, α, λ c , λ inf , λ sup ) = ∇f T k α + ∇g T c λ c + ∇g T inf λ inf + ∇g T sup λ supThen, with the proposed methodology, the multi-objective optimization problem involving the bulk carrier design can be rewritten as the following single-objective optimization problem: find X extended = (X, α, λ c , λ inf , λ sup that minimizes the voyage cost:g j ≤ 0, j = 1, 2, . . . , 21(B5)λ j ≥ 0 (B6)α i ≥ 0, i = 1, . . . , 3(B7)This single-objective optimization problem was solved by fmincon which returns the following values: X = 98.78, 12.76, 6.59, 5.22, 0.63, 14.00 T α = 0.00, 0.00, 1.00 T λ c = 0.40, 0.02, 0.03, 0.13, 0.81, 0.01, 0.00, 0.01, 0.97 T λ inf = 0.00, 0.00, 0.00, 0.02, 1.00, 0.46 T λ sup = 0.00, 0.00, 0.00, 0.00, 0.00, 0.01T Pesquisa Operacional

  k i=1 ω * i |d i |, where d i is the deviation from the goal b i for the i th objective. To model the absolute values, d i is split into positive and negative parts such that d i = d +

	OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO
	i + d -i , with d + i ≥ 0, d -i ≥ 0 and d + i d -i = 0. Deviations d -i and d + i denote underachievement and overachievement, respectively, where achievement implies that a goal has
	Pesquisa Operacional

  It should be noted that the weighting factors of the Eq. (2a), i.e., α = α 1 , α 2 , . . . , α p T , λ = λ 1 , λ 2 , . . . , λ m

	OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO
	Pesquisa Operacional

T and μ = μ 1 , μ 2 , . . . , μ l T

Table 1 -

 1 Results for the minimization of quadratic functions.

	control function	method	x 1	x 2	α 1	α 2	function computational calls time (s)
		proposed	0.2955	-4.4558 0.2323 0.7677	40	0.041
	f 1	weighted NSGA II	0.2942 0.2000	-4.4603 0.2317 0.7683 -3.7026 na na	880 6921 *	0.878 5.591
		proposed -2.5707 -3.5394 0.3626 0.6374	47	0.051
	f 2	weighted -2.5466 -3.5769 0.6423 0.3577 NSGA II -2.9838 -1.3697 na na	1001 2739 * *	0.947 2.210
		proposed -1.1142 -0.0226 0.3019 0.6981	35	0.035
	f 3	weighted -1.1242 -0.0245 0.3029 0.6971 NSGA II -0.1769 -0.3563 na na	898 8434 * * *	0.951 6.744
	compromise solution	-1.9484 -1.5833	na	na	22	0.012
	min-max solution	-1.5845 -0.5042	na	na	213	0.036

1 , x 2 , α 1 , α 2 ) = (-1.1142, -0.0226, 0.3019, 0.6981), is shown in Table 1 and in the Figure (1a). The point that minimizes the control function falls into the Pareto set resulting from the MOOP formed by the performance functions. The Pareto sets shown in Figures (1a)-(1e) are not obtained by the methodology and they were computed for comparison only by a suitable algorithm. * 141 generations completed; * * 57 generations completed; * * * 177 generations completed.

Table 2 -

 2 Results for the design optimization of a cantilever beam.

	control function	method	d	L	α 1	α 2	function computational calls time (s)
	manufacturing cost	proposed 35.0 200.0 0.8142 0.1858 weighted 34.1 200.0 0.7883 0.2117 NSGA II 34.3 200.0 na na	297 2774 2791 *	1.014 1.778 4.075
	compromise solution	29.9 200.0	na	na	75	0.0552
	min-max solution	26.5 223.6	na	na	277	0.1176

* 58 generations completed.

Table 3 -

 3 Results of the single-objective optimization of each objective function.

	Design variables and attributes	Units	Lower	Bounds Upper	Minimize	Objective functions Minimize	Maximize
	Design Variables	Length Beam Depth Draft Block Coefficient Speed	(m) (m) (m) (m) (knots)	60 10 4 3 0.63 14	600 100 40 30 0.75 18	221.9 37.0 19.8 14.6 0.72 14.0	83.9 14.0 6.7 5.4 0.63 14.0	470.6 78.4 34.4 24.8 0.63 18.0
	Ship Attributes	Annual Cost (TC) Vessel Cost (CS) Annual Cargo (CA) Deadweight Main Power Metacentric Height Froud	(10£/t) (10 8 £) (10 6 t/year) (t) (kW) (m)	3000	500000	0.797 0.186 0.750 74513 7968 4.13 0.15	2.259 0.029 0.058 3000 1336 1.35 0.25	1.519 0.795 1.271 500000 53666 10.18 0.14
		Round Trips				10.2	21.1	2.6

Table 4 -

 4 Statistics of the non-dominated solutions obtained by the random walk.

	Design variables	L (m)	B (m)	D (m)	T (m)	C B	V k (knots)
	xi σ i	352.5 54.1 30.3 19.6 0.70 75.4 11.6 6.6 4.1 0.03	16.0 1.2
	x min	88.4	13.6	6.4	5.1	0.63	14.0
	x max	493.5 77.1 40.0 26.3 0.75	18.0

Table 5 -

 5 Statistics for the Pareto set obtained by the weighted sum approach.

	Design variables	L (m)	B (m)	D (m)	T (m)	C B	V k (knots)
	xi σ i	267.5 44.0 23.8 17.4 0.74 82.2 13.7 7.4 5.1 0.03	14.4 0.8
	x min	94.2	15.7	7.6	6.0	0.63	14.0
	x max	425.5 68.9 36.6 26.3 0.75	18.0

Table 6 -

 6 Statistics for the Pareto set obtained by the NSGA II.

	Design variables	L (m)	B (m)	D (m)	T (m)	C B	V k (knots)
	xi σ i	238.0 38.3 20.3 14.6 0.71 113.5 18.7 10.1 7.0 0.05	14.8 1.1
	x min	84.2	14.0	6.4	5.0	0.63	14.0
	x max	443.1 70.2 38.3 26.3 0.75	18.0

Table 7 -

 7 Results of the design optimization with the proposed methodology.

Table 8 -

 8 Results for the design optimization of a bulk carrier.

	method	L , B, D, h (m)	C B	V (knots) function (£) control	function computational calls time (s)
	proposed	98.8, 12.8, 6.6, 5.2	0.63	14.0	12855	1031	7.872
	weighted	94.2, 15.7, 7.6, 6.0	0.63	14.0	16466	69666	73.236
	NSGA II	86.3, 13.5, 6.7, 5.4	0.64	14.0	13660	75783 *	272.125
	compromise 270.3, 45.1, 22.4, 16.4 0.64	14.0	na	239	0.206
	min-max	258.8, 43.1, 23.5, 17.2 0.75	14.0	na	519	0.474

* 797 generations completed.

Table A1 -

 A1 Port days (D P ) 2(dwt C /H R + 0.5) Round trips per year (R T PY ) 350/(D S + D P ) Pratyush S., Yang J.B., 1998, (13 constraints).

	Handling rate (H R )	8000	(t/day)
	• Voyage costs (C V ) • Ship cost (C S )	(C F + C P )R T PY 1.3(2000w 0.85 s + 3500w o + 2400P 0.8 )	(£) (£)
	Capital cost (C C )	0.2C S	(£)
	Running cost (C R )	40000dwt 0.3	(£)
	Annual cost (C A )	C C + C R + C V	(£)
	• Annual cargo (A C )	dwt c R T PY	(t/year)
	• Transportation cost (C T )	C A /A C	(£/t)
	Vertical center of buoyance (K B)	0.53T	(m)
	Metacentric radius (B M T )	(0.085C B -0.002)B 2 /(T C B )	(m)
	Vertical center of gravity (K G)	1.0 + 0.52D	(m)
	Constraints
	Explicit (4)	Implicit (9)
	0.63 ≤ C B ≤ 0.75	L/B ≥ 6
	14 ≤ V k ≤ 18	L/D ≤ 15
		L/T ≤ 19 T ≤ 0.45dwt 0.31
		T ≤ 0.7D + 0.7
		K B + B M T -K G ≥ 0.07B
		3000 ≤ dwt ≤ 500000
		F n ≤ 0.32

Table A2 -

 A2 Parsons M.G., Scott R.L., 2004 (model a), (14 constraints).

		OSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO
	Explicit (5)	Implicit (9)
	L ≤ 274.32	L/B ≥ 6
	0.63 ≤ C B ≤ 0.75	L/D ≤ 15
	14 ≤ V k ≤ 18	L/T ≤ 19 T ≤ 0.45dwt 0.31
		T ≤ 0.7D + 0.7
		Pesquisa Operacional

K B + B M T -K G ≥ 0.07B 25000 ≤ dwt ≤ 500000 F n ≤ 0.32

Table A3 -

 A3 Parsons M.G., Scott R.L., 2004 (model b), (16 constraints).

	Explicit (7)	Implicit (9)
	L ≤ 274.32	L/B ≥ 6
	T ≤ 11.71	L/D ≤ 15
	B ≤ 32.31 0.63 ≤ C B ≤ 0.75	L/T ≤ 19 T ≤ 0.45dwt 0.31
	14 ≤ V k ≤ 18	T ≤ 0.7D + 0.7
		K B + B M T -K G ≥ 0.07B
		25000 ≤ dwt ≤ 500000
		F n ≤ 0.32

Table A4 -

 A4 Hart C.G., Vlahopoulos N., 2009, (21 constraints).

	Explicit (12)	Implicit (9)
	150 ≤ L ≤ 274.32	L/B ≥ 6
	20 ≤ B ≤ 32.31	L/D ≤ 15
	13 ≤ D ≤ 25 10 ≤ T ≤ 11.71	L/T ≤ 19 T ≤ 0.45dwt 0.31
	0.63 ≤ C B ≤ 0.75	T ≤ 0.7D + 0.7
	14 ≤ V k ≤ 20	K B + B M T -K G ≥ 0.07B
		25000 ≤ dwt ≤ 500000
		F n ≤ 0.32

Table A5 -

 A5 Adopted in the present work (21 constraints).

	Explicit (12)	Implicit (9)
	60 ≤ L ≤ 600	L/B ≥ 6
	10 ≤ B ≤ 100	L/D ≤ 15
	4 ≤ D ≤ 40 3 ≤ T ≤ 30	L/T ≤ 19 T ≤ 0.45dwt 0.31
	0.63 ≤ C B ≤ 0.75	T ≤ 0.7D + 0.7
	14 ≤ V k ≤ 18	K B + B M T -K G ≥ 0.07B
		3000 ≤ dwt ≤ 500000
		F n ≤ 0.32
	Pesquisa Operacional

A NEW METHOD FOR DECISION MAKING IN MULTI-OBJECTIVE OPTIMIZATION PROBLEMS

Pesquisa OperacionalOSCAR BRITO AUGUSTO, FOUAD BENNIS and STEPHANE CARO

APPENDIX A: CONCEPTUAL MODEL FOR BULK CARRIERS

Design Variables [START_REF] Deb K | A Fast Elitist Non-Dominated Sorting Genetic Algorithm for Multi-Objective Optimization[END_REF] Length (L) (m)

Ship attributes

Steel weight


Luc Jaulin
email: luc.jaulin@ensta-bretagne.fr

Fabrice Le

Fabrice Le Bars
email: fabrice.le_bars@ensta-bretagne.fr

Sailboat

Jaulin Luc

Ihsev

Sailboat as a windmill

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Sailboat robots (see e.g. [START_REF] Neumann | Feasibility of basic visual navigation for small sailboats[END_REF] [16] [START_REF] Erckens | Navigation Strategy and Trajectory Following Controller for an Autonomous Sailing Vessel[END_REF] [2] [START_REF] Cruz | Ocean sampling and surveillance using autonomous sailboats[END_REF]) need energy for the actuators, for the sensors [START_REF] Sauze | An autonomous sailing robot for ocean observation[END_REF], for the embedded computer and for communication [START_REF] Stelzer | A study on potential energy savings by the use of a balanced rig on a robotic sailing boat[END_REF] [START_REF] Elkaim | An Energy Scavenging Autonomous Surface Vehicle for Littoral Surveillance[END_REF]. Sonar panels cannot be considered as sufficient in many situations (during the night, or in cloudy areas) and we would like to consider other sources of energy that do not depend on the sun. A wind turbine or water turbine have sometimes been used, but the energy brought cannot be considered as significant [START_REF] Rynne | Unmanned autonomous sailing: Current status and future role in sustained ocean observations[END_REF]. In this paper, we propose to use the sailboat itself as a huge wind turbine, or equivalently to reconstruct a mobile windmill. The windmill behavior of the robot assumes the boat is in a station keeping mode. Such a mode can be chosen in case where the robot has to wait for a rendezvous, or when the robot has its batteries almost empty. We assume here that the robot has only two actuators: the rudder and a blocker for tuning the sail. The corresponding controller is illustrated by Figure 1, where u 1 , u 2 correspond to the inputs (i.e., the rudder angle and the tuning of the sail) and m, θ , ψ are the outputs Fig. 1 The controller makes the robot rotating on itself as a windmill in order to produce energy.

(i.e., the GPS, the compass and the weather vane). If we consider that the blocker does not consume any energy, the only energy used for control is the rudder which consume less than 0.1W, if it is well balanced [START_REF] Stelzer | A study on potential energy savings by the use of a balanced rig on a robotic sailing boat[END_REF]. When the locker is open and the sail is opening pushed by the wind, the positive power delivered by the wind through the sail can be collected by a generator and stored inside batteries. A spring makes it possible to maintain the mainsheet tight, i.e., when the sail is in a flag mode the spring will rewind the mainsheet and close the sail. The power can be collected either at the mainsheet level via a winch or at the mastfoot. The purpose of this paper is to demonstrate the feasibility of the approach and to evaluate the amount of energy we could expect to collect with this technique.

The paper is decomposed as follows. Section 2 presents a model for the sailboat taking into account the energy and the blocker. Section 3 proposes a control strategy giving the robot windmill like behavior to produce energy. Some simulated experiments detailed on Section 4 show that it is possible to solve the station keeping problem [START_REF] Elkaim | Station Keeping and Segmented Trajectory Control of a Wind-Propelled Autonomous Catamaran[END_REF] while collecting an average of 100W for the batteries.

State space model

Different types of models exists for sailboats [START_REF] Gale | Development of a sailing dinghy simulator[END_REF], [START_REF] Guillou | Architecture multi-agents pour le pilotage automatique des voiliers de compétition et extensions algébriques des réseaux de petri[END_REF] [10]. To our knowledge, the most accurate one has been provided by Xiao and Jouffroy [START_REF] Xiao | Modeling and nonlinear heading control of sailing yachts[END_REF]. Here, to describe the dynamic of the sailboat robot, we propose a model that is sufficiently accurate to illustrate the behavior of our controller and able to give an approximation of the energy that could be collected. Classically, a state space model for a robot has the form ẋ = f (x, u)

where x is the state vector and u is the input vector. Sometimes, it is more convenient to write this state equation under the form

ẋ = g (x, z, u) z = h (x, u) where g (x, z, u) = g (x, h (x, u) , u) = f (x, u) .
The vector z contains link variables which are intermediate variables that are used to shorten the equation. Link variables can correspond to forces, angles, . . . and are often needed for the simulation to draw the robot and also to control that some feasibility state constraints are satisfied.

Model. The model is given by the following state space equations (see Figure 2).

                     (i) ẋ = v cos θ + p 1 a cos ψ (ii) ẏ = v sin θ + p 1 a sin ψ (iii) θ = ω (iv) v = f s sin δ s -f r sin u 1 -p 2 v 2 p 9 (v) ω = f s (p 6 -p 7 cos δ s )-p 8 f r cos u 1 -p 3 ωv p 10 (vi) l = u 2 if γ > 0 (vii) Ė = p 6 | f s | u 2 (1)
where the link variables are given by

                               (viii) w ap = a cos (ψ -θ ) -v a sin (ψ -θ ) (ix) ψ ap = atan2(w ap ) (x) a ap = w ap (xi) γ = cos ψ ap + cos ℓ (xii) ℓ = |δ s | if γ ≤ 0 (xiii) δ s = -tan -1 tan ψ ap if γ ≤ 0 -ℓ sign sin ψ ap otherwise (xiv) f s = p 4 a ap sin δ s -ψ ap (xv) f r = p 5 v sin u 1
This model is close to the models developed in [START_REF] Jaulin | Représentation d'état pour la modélisation et la commande des systèmes[END_REF], except that here, (a) we added the direction of the wind ψ and its amplitude a as parameters, (b) the control is not anymore the sail angle, but the length of the mainsheet, which is more realistic, (c) the speed of the robot is not considered as small compared to the true wind (the notion of apparent wind has thus to be introduced), (d) the angular friction now depends on the speed, which is more consistent with actual sailboats and (e) the length ℓ of the mainsheet and the energy of the batteries E are introduced as state variables. All quantities are expressed using the international unit system. For simplicity, the lenght ℓ of the mainsheet is expressed in radian (rad). ℓ corresponds to the absolute value of the maximal angle δ s that could reach the sail when the mainsheet is tight. Let us now describe more deeply all variables involved in this model.

Inputs. The sailboat has two inputs. The first input u 1 = δ r is the angle between the rudder and the sailboat. The second input u 2 corresponds to the blocker. When u 2 = 1, the locker is unblocked and the length of the mainsheet ℓ may increase (if the direction of the wind allows it). Otherwise, u 2 = 0 and the blocker is active.

State variables. The state variable occurring in our model ( 1) are x, y, θ , v, ω, ℓ, E where (x, y) are coordinates of the robot, θ is its heading, v is its speed along the main axis, ω is its rotational speed. The energy of the batteries E will increase with time. The length of the mainsheet ℓ corresponds of the maximal angle of the sail. In the particular case where the mainsheet is tight, it can be computed from the state variables θ , u 2 , ψ, v and thus it cannot be considered as a state variable anymore. Therefore, the dimension of the state vector (either 6 or 7) changes with time. The sailboat thus corresponds to an hybrid system.

Parameters. In our model, p 1 is the drift coefficient, p 2 is the tangential friction, p 3 is the angular friction, p 4 is the sail lift, p 5 is the rudder lift, p 9 is the mass of the boat and p 10 is the mass moment of inertia. The distances p 6 , p 7 , p 8 are represented in Figure 2. All parameters p i are assumed to be known exactly. Two other quantities should also be considered as parameters: the speed a of the wind and its direction ψ.

Link variables. These variables are used to shorten the expression of the state equations. (viii) The vector w ap corresponds to the apparent wind expressed in the robot frame. The amplitude (ix) and the angle (x) of w ap (in the robot frame) are denoted by a ap and ψ ap . (xi) The coefficient γ is positive if the mainsheet is tight. (xii) In this case, ℓ is a state variable and its evolution obeys to the differential equation l = u 2 . Otherwise, the mainsheet is tight, ℓ is a link variable and its value is equal to |δ s |. This change of status of ℓ is typical of what happen for hybrid State equations. The two first equations (i),(ii) of (1) express that the boat follows its heading, but always looses with respect to the wind. Equations (iv) and (v) are obtained using the Newton laws. Equation (vi) tells us that the length ℓ of the mainsheet can only increases when the sail is inflated and when the blocker is off (i.e., u 2 = 1). Equation (vii) provides the power delivered to the batteries: when u 2 = 1, and f s = 0, the sail opens with an angular velocity of 1rad. sec -1 and the power collected is p 6 | f s |. Note that this model for the sailboat could be made more realistic by adapting the modeling tools described by Fossen in the context of marine vessel [START_REF] Fossen | Guidance and Control of Ocean Vehicles[END_REF] to sailboats. But to our knowledge, realistic state equation for sailboats do not exist yet.

Controller

A classical approach to build controllers is to take a realistic model of the system to be controlled (such as [START_REF] Gale | Development of a sailing dinghy simulator[END_REF] for the control of sailboats) and then to use classical control methods to get the controller. Here, we follow a pragmatic approach influenced by the potential field strategy proposed by [START_REF] Petres | Reactive path planning for autonomous sailboat[END_REF] for sailboat robots. Our sailboat robot is assumed to have three sensors and two actuators. The controller will have some parameters which are easy to tune, two state variables q ∈ {1, 2, 3, 4} , and

t 0 ∈ R + , two outputs u 1 ∈ [-π 4 , π 4 ], u 2 ∈ {0, 1} and three inputs m ∈ R 2 , θ ∈ R, ψ ∈ R.
Let us now describe all these variables.

Sensors (which also correspond to the input of the controller). The heading θ of the robot is measured by a compass. The angle of the wind ψ is returned by a weather vane (even if this sensor can sometimes be omitted as shown in [START_REF] Xiao | A wind-independent control strategy for autonomous sailboats based on voronoi diagram[END_REF]). The position m is given by a GPS. These sensors are the inputs of our controller.

Actuators (which also correspond to the output of the controller). The inputs of the robot are the angle

u 1 ∈ [-π 4 , π 4
] of the rudder and the blocker u 2 ∈ {0, 1} which makes it possible to tune indirectly the length of the mainsheet.

Parameters. δ max r is the maximal angle of the rudder is taken as

δ max r = π 4 .
ζ is the close hauled angle. For the simulation, we will choose ζ = 1rad.

State variable. The discrete variable q ∈ {1, 2, 3, 4}, gives three modes: the traction (q = 1), the rewind (q = 2) and the positioning (q = 3, q = 4). The start time t 0 corresponds to the time at which the timer is started when the controller is in the positioning mode.

The basic idea of the controller is to decompose the plane into three cones, the intersection of which corresponds to the origin, as illustrated by Figure 4. In the mill cone (painted gray), the robot rotates as a windmill to produce energy. In the mill cone, the robot slowly moves downwind. The points m that are inside the mill cone satisfy the inequality cos(ψ-arg(m)) >cos ζ . In the hatched cone, the controller will carry favor to the close-hauled heading π + ψζ . In the white cone, it will prefer the heading π + ψ + ζ .

We now give the details of the controller which is clearly influenced by the line following controller proposed in [START_REF] Jaulin | Suivi de route pour un robot voilier[END_REF] [START_REF] Jaulin | A simple controller for line following of sailboats[END_REF] and already experimented [START_REF] Bars | An experimental validation of a robust controller with the VAIMOS autonomous sailboat[END_REF] on the sailboat robot VAIMOS. Function in: m, θ , ψ; out: u 1 , u 2 ; inout: q,t 0 1 if (q = 1 and ψ ∼ θ ) then q = 2; 2 if (q = 2 and ψ ∼ θ + π) then 3 if cos(ψ-arg(m)) >cos ζ 4 t 0 = t 5 if (sin(ψarg(m)) > 0) then q = 3; else q = 4 6 else q = 1 7 if (q ∈ {3, 4} and tt 0 > 30) then q = 1;

8 if q = 1 then θ = ψ 9 if q = {2, 4} then θ = π + ψ + ζ 10 if q = 3 then θ = π + ψ -ζ 11 if ( cos θ -θ ≤ 0 or q ≤ 2) 12 then u 1 = π 4 .sign sin θ -θ 13 else u 1 = π 4 . sin θ -θ 14 l = π 2 . cos(ψ-θ )+1 2 . 15 if l > ℓ then u 2 = 1 else u 2 = 0.
Steps 1 to 7 correspond to the discrete event part of our hybrid controller. It is illustrated by the Petri net of Figure 5. The gray places represent actual states and white places represent fake states (the token leaves a fake place as soon as it enters Fig. 5 Petri net associated with our windmill sailboat in it). Bold arrows have a higher priority and are necessary to make the Petri net deterministic. Let us now describe the different discrete states for q.

• Traction (q = 1). The controller opens the sail and maneuvers to go downwind (see Step 8) as fast as possible (see Steps 11,[START_REF] Jaulin | A simple controller for line following of sailboats[END_REF]. The controller escapes the state q = 1 at Step 1 as soon as ψ ∼ θ (i.e., ψ = θ ± 2kπ). When q = 1, the mainsheet pulls the generator and energy is produced. • Rewind (q = 2). The controller makes the boat rotating toward the wind, in order to close the sail. When the robot is upwind (ψ ∼ θ + π), then the rewind mode terminates (see Step 2). If the robot is inside the mill cone, the controller goes to the state q = 1 at Step 6. Otherwise, depending on which cone the robot is, the controller chooses the states q = 3 or q = 4 at Step 5. • Positioning (q ∈ {3, 4}). The controller chooses a close-hauled heading for 30 second, in order to bring closer to the mill cone.

Steps 8 to 10 provide the desired heading θ to follow, depending of the value of q. When q = 1, the controller asks to go downwind (Step 8), Otherwise, it ask to go to a close hauled mode (Steps 9,10). Steps 11,12,13 tune the rudder (see [START_REF] Jaulin | Suivi de route pour un robot voilier[END_REF] for more explanations). If the boat goes toward the wrong direction (cos θθ ≤ 0) or if q ∈ {1, 2}, the rudder at its maximum, i.e., ± π 4 . otherwise, a proportional control is proposed (Step 13). For tuning of the sail, we propose to take a Cardioid model [START_REF] Jaulin | An interval approach for stability analysis; Application to sailboat robotics[END_REF] at Step 14 to compute the right angle for the sail. At Step 15, the controller suggests to open the sail (u 2 = 1) in order to reach the desired length l, by opening the blocker. This will mainly happen when q = 1.

Test-case

In order to illustrate the principle of the controller, we now propose a simulation of the controlled sailboat robot. For the parameters, we have chosen

p 1 = 0.05, p 2 = 0.2 kg • s -1 , p 3 = 6000 Kg • m, p 4 = 1000 kg • s -1 , p 5 = 2000 kg • s -1 , p 6 = 1m, p 7 = 1m, p 8 = 2m, p 9 = 300 kg, p 10 = 10000 Kg • m 2 .
Except for p 6 , these values correspond approximately to the coefficients of the sailboat robot VAIMOS [START_REF] Jaulin | Suivi de route pour un robot voilier[END_REF]. The value for p 6 is almost zero for VAIMOS, due to the balestron rig (or balanced rig). For our application, it is important to have an important p 6 for the energy production. For the speed a of the wind and its direction ψ, we took.

a = 4 m • s -1 , ψ = π.
For the parameter of the controller, we have chosen ζ = 1 rad. The resulting trajectory is illustrated by Figure 6 where the robot is initialized at m = (-400, 200) (small black disk). The trajectory corresponds to 30 minutes and the average of the collected power is around 93W. On the picture, we clearly see that on the mill cone, the boat rotates and move downwind. As soon as it goes outside the cone, it comes back to the mill cone choosing the right close hauled heading. The executable program and the C++ source code of the simulator that has generated Figure 6 can be found at http://www.ensta-bretagne.fr/jaulin/mill.html

Remark. Betz's law [START_REF] Betz | Introduction to the Theory of Flow Machines[END_REF] claims that the maximum power that can be extracted from the wind, independent of the design of a wind turbine, is given by

P Betz = 8 27 ρSa 3 ,
where S is the surface of the turbine, ρ is the fluid density and a is the speed of the wind. From this formula, we can deduce that to get the same power than that collected by the batteries in our test-case, a wind turbine with a diameter of 2.4 meters would be needed. Of course, such a wind turbine would change significantly the dynamic performances of the sailboat.

Conclusion

In this paper we have presented a new concept that allows a sailboat to take advantage of a station keeping mode in order to charge its batteries. The basic idea is to transform the sailing boat into a windmill using a hybrid controller. When the wind Fig. 6 Trajectory of the sailboat robot which tries to remain inside the circle and also to collect energy from the wind inflates the sail, the mainsheet pull a generator that produces energy for the batteries. A test-case has shown that a mean power of 93W could be collected for a wind speed equal to 4 m•s -1 . All computations made by the controller can be performed using any cheap and low-powered microcontroller.

Fig. 2

 2 Fig. 2 Sailboat to be used as a windmill

Fig. 3

 3 Fig. 3 Differential graph of our sailboat robot

  Figure 3 represents the differential graph of the state equations. The state variables are represented by grey nodes and the inputs by square nodes. The integral relations are represented by bold arrows and the link relation by dotted arrows. The two bold dotted arrows illustrate that the differential dependency between u 2 , γ and ℓ are valid for some conditions only.

Fig. 4

 4 Fig.4Principle of the controller to maintain the robot around the origin still trying to spin as a mill in order to charge the batteries.

Acknowledgements The robot VAIMOS is the result of a collaboration between LPO (Laboratoire de Physique des Océans), RDT (Recherches et Développements Technologiques) of Ifremer (Institut Franc ¸ais de Recherche pour l'Exploitation de la Mer) and ENSTA Bretagne (Ecole Nationale Supérieure de Techniques Avancées). The authors render thanks to all people involved in the project: Y. Auffret, S. Barbot, L. Dussud, B. Forest, E. Menut, S. Prigent, L. Quemeneur, P. Rousseaux (RDT, Ifremer); F. Gaillard, T. Gorgues, O. Ménage, J. Moranges, T. Terre (LPO) and B. Clément, Y. Gallou, O. Reynet, J. Sliwka and B. Zerr (ENSTA Bretagne).