A multiobjective evolutionary approach for linearly constrained project selection under uncertainty

European Journal of Operational Research 179 (2007) 869–894 www.elsevier.com/locate/ejor

A multiobjective evolutionary approach for linearly constrained project selection under uncertainty Andre´s L. Medaglia

a,*

, Samuel B. Graves b, Jeffrey L. Ringuest

b

a

b

Departamento de Ingenierı´a Industrial & Centro para la Optimizacio´n y Probabilidad Aplicada, Universidad de los Andes, Bogota, Colombia Operations, Information, and Strategic Management Department, Boston College, Chestnut Hill, MA, USA Received 13 May 2004; accepted 7 March 2005 Available online 6 December 2005

Abstract In the project selection problem a decision maker is required to allocate limited resources among an available set of competing projects. These projects could arise, although not exclusively, in an R&D, information technology or capital budgeting context. We propose an evolutionary method for project selection problems with partially funded projects, multiple (stochastic) objectives, project interdependencies (in the objectives), and a linear structure for resource constraints. The method is based on posterior articulation of preferences and is able to approximate the efficient frontier composed of stochastically nondominated solutions. We compared the method with the stochastic parameter space investigation method (PSI) and illustrate it by means of an R&D portfolio problem under uncertainty based on Monte Carlo simulation.  2005 Elsevier B.V. All rights reserved. Keywords: Multi-objective optimization; Project selection; Multi-objective evolutionary algorithms

1. Introduction In the project selection problem a decision maker allocates limited resources to a set of competing projects. In addition, while selecting projects and allocating resources, the decision maker must take into consideration one or more corporate goals or objectives.

*

Corresponding author. Tel.: +57 1 3394949; fax: +57 1 332432. E-mail addresses: [email protected] (A.L. Medaglia), [email protected] (S.B. Graves), [email protected] (J.L. Ringuest). URLs: http://wwwprof.uniandes.edu.co/~amedagli (A.L. Medaglia), http://www2.bc.edu/~graves/ (S.B. Graves), http:// www2.bc.edu/~ringuest/ (J.L. Ringuest). 0377-2217/$ - see front matter  2005 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2005.03.068

870

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

There are two major taxonomies of the project selection literature. The first taxonomy is based on the application context. Within this classification, the project selection problem arises primarily in the context of research and development (R&D) [1,15,19,33,34], information technology (IT) [25,28,29,37–40,43], and capital budgeting [17,31,48,49]. The second taxonomy is based on the nature of the tools used for solving the project selection problem. This taxonomy classifies project selection models into one of two streams: the management science stream and the financial optimization stream. For more information about the project selection problem refer to [18]. In this paper, we propose a method for solving any project selection problem with the following characteristics: 1. Partially funded projects. Projects can be partially funded, or conversely, fully or not funded at all (i.e., zero-one). The proposed method allows the decision maker to engage in partially funded projects. For a sample of zero-one project selection refer to [22,39]. 2. Multiple objectives. The decision maker may wish to simultaneously satisfy several corporate goals (e.g., minimize time-to-market and maximize economic return). Often, these goals will be in conflict, resulting in a more complicated, yet realistic, decision making process. The proposed model takes into consideration multiple (possibly conflicting) objectives. 3. Posterior articulation of preferences. The methods for multiobjective project selection can be classified based on how the decision maker expresses his/her preferences in the decision making process. Some methods have been proposed based on prior [23] and progressive [36] articulation of preferences, but few based on posterior articulation of preferences. The proposed method, based on posterior articulation of preferences, requires the least amount of information from the decision maker and could be completely automated. 4. Uncertainty in the objectives. Some elements of the project selection process may be uncertain. For instance, some uncertain elements that could affect a projects evaluation are interest rates, opportunity cost, risk, revenues, inflation, and cash flows, among others. During project evaluation, the proposed method allows uncertain elements in the objectives to be included as random coefficients in complex expressions (even nonlinear). Moreover, the method also allows the treatment of the whole project as a black-box. In this case, the method simulates the project and evaluates its performance measures (objectives) subject to random events. 5. Project interdependencies. In the project selection literature, there are mainly three types of interdependencies: resource, benefit, and technical interdependencies [38]. Interdependencies preclude the linearity of the models. For example, suppose there are two projects, A and B, with interdependent benefits of $500 and $1000, respectively. The benefits of implementing both projects could exceed $1500 (=$500 + $1000) due to the additional benefits accrued from a pooling effect. This benefit interdependence could be modeled using non-linear terms as in [15,38]. The proposed model allows project interdependencies to be modelled in the objective functions. 6. Linearly-constrained resources. Project selection problems are constrained by their use of scarce resources (e.g., capital, labor, space requirements, and personnel, among others). Linear optimization models [12] have been popular due in part to the fact that many problem conditions are suitable to be expressed in linear form. The proposed method uses common, yet expressive, linear relations to model the restrictions on the various resources. The stochastic parameter space investigation (PSI) method [35] is able to solve the project selection problem with the above characteristics. The stochastic PSI proceeds as follows: (1) project allocations are generated based on box constraints on the variables representing them (i.e., lower and upper bounds); (2) if a project allocation does not satisfy the constraints on the resources, that configuration is discarded; (3) if a project allocation is feasible, the objectives are simulated and evaluated; (4) statistics of central tendency

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

871

and deviation are computed based on the simulation results (i.e., mean and kth percentile); and (5) if the project allocation is nondominated (based on Definition 1 in Section 2), it is preserved for future presentation to the decision maker at the end of the procedure (otherwise, it is discarded). The stochastic PSI method has several advantages. Among them, it can be applied to a very large class of project selection problems, it is very easy to implement, makes no assumption about the behavior of the decision maker, and accounts for uncertainty in the objectives. However, all the flexibility of the stochastic PSI method does not come for free. From the work on the deterministic PSI by Steuer and Sun [47] and our own experience, the limitations of the stochastic PSI are: the method can only handle project selection problems of limited size (less than ten variables), its computational efficiency is dependent on the relative size of the feasible region to the size of the hypercube that encloses the feasible region (given by the bounds on the variables), and it requires explicit and tight bounds on the values of the decision variables. In this paper we present a new multiobjective evolutionary algorithm that serves as an alternative to the stochastic PSI method. The proposed method is applicable to project selection problems with partially funded projects, multiple (stochastic) objectives, project interdependencies (in the objectives), and linearly-constrained resources. Furthermore, being a method based on posterior articulation of preferences, it is able to approximate the efficient frontier composed of stochastically nondominated solutions. This paper complements previous work done in [26]. This paper is organized as follows. In Section 2 we present the proposed multiobjective evolutionary algorithm. In Section 3 we present our computational experience with the proposed algorithm. In Section 3.1 we tune the parameters; in Section 3.2 we present a comparison between the stochastic PSI and the proposed evolutionary approach in problems found in the literature; and in Section 3.3 we show an R&D portfolio optimization example that illustrates the flexibility of our approach. This example uses Monte Carlo simulation for evaluating a given set of project allocations. Finally, conclusions and current research directions are given in Section 4.

2. Proposed method 2.1. Project selection as a stochastic multiobjective linearly-constrained optimization problem Consider the project selection problem seen as the following stochastic multiobjective linearly-constrained optimization problem: max ðf1 ðx; xÞ; . . . ; fk ðx; xÞ; . . . ; fK ðx; xÞÞ subject to x 2 X;

ð1Þ ð2Þ

where x is a solution (vector of n decision variables) representing a nth dimensional array of project allocations; fk(x, x) is the objective function evaluating the kth criterion and x represents the stochastic effect in the objectives; and X ¼ fx 2 Rn : gi ðxÞ 6 bi ; i ¼ 1; . . . ; m; x P 0g is the polyhedral space of feasible project allocations defined by the intersection of the linear constraints gi(x) 6 bi (i = 1, . . . , m). Furthermore, without loss of generality we assume X 5 ; and bounded. Definition 1 (Stochastic Pareto solution). Let x, y 2 X be a pair of feasible solutions for the project selection problem described in (1) and (2). We say that y stochastically dominates x (i.e., y  x) if and only if the following two conditions hold: (i) E[fk(y, x)] P E[fk(x, x)] and P{fk(y, x) P Tk} P P{fk(x, x) P Tk}, for all k; and (ii) there exists a k such that E[fk(y, x)] > E[fk(x, x)] and/or P{fk(y, x) P Tk} > P{fk(x, x) P Tk},

872

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

where E[fk(x, x)] is the expected value of the kth objective and P{fk(x, x) P Tk} is the probability that fk(x, x) exceeds (or at least reaches) the target value Tk specified by the decision maker. A feasible solution which is not dominated by any other solution is called a Pareto solution or simply nondominated. Definition 2 (Pareto optimal set, P ). The Pareto optimal set P is composed of the feasible solutions which are not dominated by any other solution. Therefore, P ¼ fx 2 X : there is no x0 2 X; x0  xg. Definition 3 (Efficient (or Pareto) Frontier, PF ). The efficient frontier is the image of the Pareto optimal set in the objective space, this is, PF ¼ ffðxÞ ¼ ðf1 ðxÞ; . . . ; fk ðxÞ; . . . ; fK ðxÞÞ : x 2 P g. 2.2. Multiobjective evolutionary algorithms An evolutionary algorithm is a stochastic search procedure inspired by the evolution process in nature. Individuals (i.e., solutions) evolve and the fitter ones have a better chance of reproduction and survival in the selection process. The reproduction mechanisms favor the characteristics of the stronger parents and hopefully produce better children guaranteeing the presence of those characteristics in future generations. Several lines of research on algorithms sharing the same evolutionary principle have matured: evolutionary programming [13], evolutionary strategies [32,44], and genetic algorithms [20,16]. Evolutionary algorithms have been applied to solve complex problems where traditional optimization algorithms fail. For further information on evolutionary algorithms we refer the reader to [5–7,27]. Evolutionary algorithms are well suited for solving multiobjective optimization problems. According to Zitzler and Thiele [50], the use of a population of solutions favors the exploration of a vast set of alternatives in parallel. In addition, Coello [8] highlights the fact that evolutionary algorithms are able to find members of the Pareto optimal set in a single run, instead of multiple runs often used in mathematical programming-based algorithms. Multiobjective evolutionary algorithms (MOEA) have evolved since the pioneering work of Schaffers vector evaluation genetic algorithm (VEGA) [42]. In later years, the most active MOEA line of research has become the design of algorithms based on posterior articulation of preferences. Basically, this class of MOEA algorithms evolve the population of solutions towards the efficient frontier. The most popular MOEAs based on posterior articulation of preferences are the multiobjective genetic algorithm (MOGA) [14], the Niched Pareto genetic algorithm (NPGA) [21], the nondominated sorting genetic algorithm (NSGA [45] and NSGA II [10]), the Strength Pareto evolutionary algorithm (SPEA) [50], and the Pareto archived evolution strategy (PAES) [24]. For more information on multiobjective evolutionary algorithms the reader is referred to Coello et al. [9]. 2.3. Project selection problem: A MOEA approach The proposed algorithm described in this section borrows elements from NSGA-II [10], an efficient multiobjective evolutionary framework based on posterior articulation of preferences, and from linear optimization theory [12]. The key elements implemented in this NSGA-based algorithm are: 1. Elitist strategy. Elitism promotes the preservation of the best individuals. In a multiobjective context it means that (stochastically) nondominated individuals are passed from generation to generation. 2. Parameter-less diversity. Often, evolutionary algorithms depend on tuning and setting a large number of parameters for preserving diversity in their population. The drawback of this approach is that, although it makes the algorithm flexible, a great deal of experience is needed to know what knob to tweak. The proposed algorithm uses a parameter-less approach based on crowding distance that ensures diversity in the population.

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

873

3. Fast stochastic dominance mechanism. Algorithms based on posterior articulation of preference often impose a burden on keeping nondominated solutions from generation to generation. The proposed algorithm uses a fast mechanism proposed by Deb et al. [10] that alleviates this computational bottleneck. 4. Efficient constraint-handling mechanism. The constraint handling mechanism has been designed to take advantage of the linearly constrained solution space. The algorithm uses fundamental concepts such as the resolution theorem and the simplex and dual simplex methods [12].

2.3.1. The multiobjective evolutionary algorithm Let PðtÞ ¼ fx1 ðtÞ; . . . ; xP ðtÞjxp ðtÞ 2 X; p ¼ 1; . . . ; P g be the population at generation t, where xp(t) is the pth chromosome, P is the maximum number of individuals in the population, and X is the feasible set of project allocations as described in (2). Let Cc ðtÞ be the children population after recombining PðtÞ; Cm ðtÞ be the population of fresh individuals injected at generation t; and CðtÞ the whole children population composed of Cc ðtÞ [ Cm ðtÞ. Let EðtÞ be the expanded population composed of the population and their children (i.e., PðtÞ [ CðtÞ). Let FðtÞ ¼ F1 ðtÞ [    [ Fl ðtÞ [    [ FLðtÞ ðtÞ be the set of efficient frontiers organized by levels of nondomination. The efficient frontier F1 ðtÞ is formed by the nondominated individuals of EðtÞ; F2 ðtÞ is the frontier with the nondominated individuals of EðtÞ, once the individuals from F1 ðtÞ are excluded; F3 ðtÞ is the frontier with the nondominated individuals of EðtÞ, once the individuals from F1 ðtÞ [ F2 ðtÞ are excluded; and so on. The last frontier SLðtÞ1is FLðtÞ ðtÞ, it contains the nondominated individuals from EðtÞ, after excluding the individuals from l¼1 Fl ðtÞ. Furthermore, the individuals from FLðtÞ ðtÞ do not dominate any individual from EðtÞ. Note that the total number of efficient frontiers L(t) is not fixed, depends on the individuals found in EðtÞ, and varies from generation to generation. Fig. 1 illustrates a population of 10 individuals classified in three efficient frontiers, namely, F1 , F2 , and F3 . This example assumes, without loss of generality, a given generation and two deterministic functions f1(Æ) and f2(Æ), that are to be maximized. Algorithm 1 illustrates step by step the proposed evolutionary mechanism, where T is the maximum number of generations. At the end, PðT Þ contains an approximation of the efficient frontier composed of stochastically nondominated solutions. Algorithm 1 (NSGA-based multiobjective evolutionary algorithm for project selection) 1: t 1 2: Initialize Pð1Þ {see Section 2.3.3} 1

2

Fig. 1. Example of a ten-individual population classified by levels of nondomination.

874

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

3: Evaluate Pð1Þ {see Section 2.3.6} 4: while t 6 T do 5: Recombine PðtÞ to generate Cc ðtÞ {see Section 2.3.4} 6: Inject fresh individuals Cm ðtÞ {see Section 2.3.5} 7: CðtÞ ¼ Cc ðtÞ [ Cm ðtÞ 8: Evaluate CðtÞ {see Section 2.3.6} 9: EðtÞ ¼ PðtÞ [ CðtÞ 10: Sort EðtÞ by stochastic nondomination to generate FðtÞ {see Section 2.3.7} 11: Select Pðt þ 1Þ from FðtÞ {see Section 2.3.7} 12: t t+1 13: end while The following sections contain further detail on several aspects of the algorithm. 2.3.2. Solution encoding In evolutionary algorithms, solutions are conveniently encoded into chromosomes or individuals. An arbitrary feasible solution x 2 X (see (1) and (2)) is represented in the proposed evolutionary algorithm by the pth individual xp(t) of PðtÞ using floating-point encoding [27]. In this type of encoding (or representation), xp ðtÞ ¼ ½xp1 ðtÞ; . . . ; xpj ðtÞ; . . . ; xpn ðtÞ 2 X, where xpj ðtÞ 2 R for j = 1, . . . , n. 2.3.3. Generating new individuals Using a mechanism based on fundamental concepts of linear programming we guarantee that every generated individual is feasible. Algorithm 2 illustrates the mechanism for generating the initial population. Using Algorithm 3 we generate an individual that maps to an extreme point of polytope X. If the individual xp(1) is already in the initial population Pð1Þ, we generate a sequence of linear convex combinations between this individual and randomly picked individuals xq ð1Þ 2 Pð1Þ, until a new individual xp(1), not already in Pð1Þ, is generated. Algorithm 2 (Algorithm for generating the initial population Pð1Þ) 1: p 1 2: while p 6 P do 3: xp(1) extreme(X) {See Algorithm 3} 4: while xp ð1Þ 2 Pð1Þ do 5: Pick randomly xq ð1Þ 2 Pð1Þ 6: Generate randomly k from a uniform distribution between 0 and 1. 7: xp(1) kxp(1) + (1  k)xq(1) 8: end while 9: Pð1Þ Pð1Þ [ fxp ð1Þg 10: p p+1 11: end while One of the key components of Algorithm 2 is line 3, where an extreme point of the polytope X is generated by means of Algorithm 3. Without loss of generality, assume that a slack si P 0 is added to each 0 constraint gi(x) 6 bi of X (i = 1, . . . , m). After adding these slacks, we obtain a new space  X ¼  fx 2 x n m 0 0 R ; s 2 R : gi ðxÞ þ si ¼ bi ; i ¼ 1; . . . ; m; x; s P 0g. In matrix form, we write X ¼ x ¼ 2 s

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894



875

Rnþm : Gx0 ¼ b; x0 P 0 . The columns of matrix G are rearranged in basic and nonbasic columns, namely, G = [GBjGN]. Associated with these columns, we have basic and nonbasic variables of x 0 , x0B and x0N (=0), respectively. In lines 1–4 of Algorithm 3, a potential basis GB is generated by randomly choosing m out of the m + n columns of G. If GB is full-row rank (i.e., rank(GB) = m), a basis is obtained; if not, we try another set of randomly picked columns. However, having a basis does not guarantee that a basic feasible solution is achieved. In lines 5–11 of Algorithm 3, a basic feasible solution of X 0 is obtained starting from basis GB. If GB does not generate a feasible solution (i.e., x0B ¼ G1 B b P 0), we perform dual-simplex pivoting [12], until one is achieved. Once we obtain a feasible solution x 0 2 X 0 , we map it back to the original variable space, namely, xp(t) 2 X, which in turn becomes an individual of Pð1Þ. By design, the initial population has a strong bias toward extreme points. The idea behind this concept, is that because of the resolution theorem of linear programming [12], extreme points are sufficient to generate any feasible solution of the polytope X. Algorithm 3 (Algorithm extreme for generating feasible extreme points of X) Input: Polytope definition X Output: Feasible (extreme-point) individual xp(t) 2 X 1: Initialize GB = 0 2: while rank(GB) < m do 3: Generate a potential basis GB by randomly picking m out of m + n columns of G 4: end while 5: x0B ¼ G1 B b 6: while x0B j0 do 7: Pick a column of GB which is leaving the basis {follows dual simplex pivoting scheme} 8: Pick a column from GN which is entering the basis {follows dual simplex pivoting scheme} 9: Update the basis GB with the leaving and entering columns 10: x0B ¼ G1 B b 11: end while 12: Project x 0 2 X 0 to xp(t) 2 X

2.3.4. Recombining individuals via crossover The purpose of the crossover operation is to generate new individuals (children) from parents in PðtÞ. This operation has been designed so that it generates a single child with traits of both parents. Let XðtÞ be the crossover pool of parents and pc the chance of any individual in PðtÞ of being selected as a parent. Algorithm 4 shows how the crossover operation is performed. Again, the idea behind the operation is the convex combination between feasible solutions of X. Algorithm 4 (Algorithm for the crossover operation) Input: Population of parents PðtÞ and probability of crossover pc Output: Children from the crossover operation Cc ðtÞ 1: XðtÞ ; 2: Cc ðtÞ ; 3: p 1 4: while p 6 P do 5: Generate randomly u from a uniform distribution between 0 and 1.

876

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

6: if u < pc then 7: XðtÞ XðtÞ [ fxp ðtÞg 8: end if 9: p p+1 10: end while 11: while jXðtÞj > 1 do 12: Pick randomly two distinct individuals from XðtÞ, namely xp(t) and xq(t) 13: Generate randomly k from a uniform distribution between 0 and 1. 14: x(t) = kxp(t) + (1  k)xq(t) 15: Cc ðtÞ Cc ðtÞ [ xðtÞ 16: XðtÞ XðtÞ n fxp ðtÞ; xq ðtÞg 17: end while 2.3.5. Injecting fresh individuals A successful evolutionary strategy balances exploitation and exploration. On one hand, exploitation allows an evolutionary algorithm to get the most out of a promising vicinity of the search space. The crossover operation described in Section 2.3.4 was designed for this purpose. On the other hand, exploration allows an evolutionary algorithm to jump to unexplored regions of the solution space. Being extreme points the most valuable individuals (due to the reasons explained in Section 2.3.3), the mechanism described in Algorithm 5 guarantees that individuals with information of (or closely related to) the extreme points of X are constantly injected into the population. The fraction of the population that will be injected with information closely related to the extreme points of X is denoted by pm (henceforth we call it injection rate). Algorithm 5 (Algorithm for injecting fresh individuals 5) Input: Population of parents PðtÞ and injection rate pm Output: Population of fresh individuals Cm ðtÞ 1: Cm ðtÞ ; 2: p 1 3: while p 6 bpm Æ Pc do 4: xp(t) extreme(X) {See Algorithm 3} 5: while xp ðtÞ 2 PðtÞ do 6: Pick randomly xq ðtÞ 2 PðtÞ 7: Generate randomly k from a uniform distribution between 0 and 1 8: xp(t) kxp(t) + (1  k)xq(t) 9: end while 10: Cm ðtÞ Cm ðtÞ [ fxp ðtÞg 11: p p+1 12: end while 2.3.6. Evaluating populations A given individual xp(t) is evaluated by the vector of objectives f(xp(t)) = (f1(xp(t), x), . . . , fk(xp(t), x), . . . , fK(xp(t), x)). However, this vector of objectives is a random vector due to the stochastic effect x present in the project evaluation process (e.g., risk, inflation, and variation on the cash flows, among others). Therefore, we have chosen to quantify the merit of each project allocation using a vector of expected values (E[f1(xp(t))], . . . , E[fk(xp(t))], . . . , E[fK(xp(t))]) and as a proxy of its variation, the vector of percentiles (T1, . . . , Tk, . . . , TK), where P{fk(xp(t), x) P Tk} = Fk and Fk is the probability level specified by the decision maker.

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

877

It is often the case that either it is difficult or impossible to find closed-form solutions for (E[f1(xp(t))], . . . , E[fk(xp(t))], . . . , E[fK(xp(t))]) and (T1, . . . , Tk, . . . , TK). In this situation, it is reasonable to resort to a black-box approach as illustrated in Fig. 2. The input is a given configuration of project allocations xp(t). This set of allocations is evaluated by means of a Monte Carlo simulator. Statistics are gathered until a significant amount of replications are performed. The output is an estimate of the vector of expected values (E[f1(xp(t))], . . . , E[fk(xp(t))], . . . , E[fK(xp(t))]) and the percentiles (T1, . . . , Tk, . . . , TK). 2.3.7. Sorting and selecting the population for the next generation In Section 2.3.1, we describe how the expanded population EðtÞ , PðtÞ [ CðtÞ is sorted in levels of stochastic nondomination FðtÞ , F1 ðtÞ [    [ Fl ðtÞ [    [ FLðtÞ ðtÞ (refer to line 10 of Algorithm 1). Using the nondomination sorting algorithm presented in Deb et al. [10], the complexity of this process is 2 OðKjEðtÞj Þ. We select the population for the next generation P(t + 1) from an enlarged sampling space, namely, the expanded population sorted in fronts FðtÞ. The selection mechanism used in the algorithm is known as (l + k) selection [4]. In (l + k) selection, l parents (i.e., PðtÞ) and k children (i.e., CðtÞ) compete for survival, being the best P individuals from (l + k) (i.e., EðtÞ) the chosen ones. The best individuals are selected from EðtÞ once the expanded population has been sorted in levels of stochastic nondomination FðtÞ. First, all the individuals in the first frontier F1 ðtÞ are selected; then, the individuals in the second frontier F2 ðtÞ are selected; and so on. The last frontier that does not fit completely in Pðt þ 1Þ (limit of P individuals) is sorted based on the crowding distance proposed by Deb et al. [10]. After sorting this last frontier, the best individuals (based on crowding distance) are picked to complete Pðt þ 1Þ. The (l + k) selection mechanism is illustrated in Algorithm 6. Fig. 3 illustrates what is done in the crowding distance calculation of Algorithm 6. In this example, the first frontier that does not fit in Pðt þ 1Þ has six solutions. Both x1 and x6, the extremes of the frontier, are assigned crowding distances of b(x1) = 1 and b(x6) = 1 (so they will almost certainly be selected).

Fig. 2. Black box approach to project evaluation.

[

E [ f1 ( x )

1

2 3

4 5

6

Fig. 3. Example of crowding distance calculation.

[

E [ f2 ( x)

878

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

Because individual x5 is in a less crowded region of the front than individual x3, b(x5) > b(x3), thus the selection process will favor x5 over x3. Note also that x5 is bounded by a larger cuboid than x3. The rationale is that individual x5 provides more unique information about the front than individual x3. For instance, individuals x2, x3, and x4, contain similar information, so any of those could represent the others. Algorithm 6 ((l + k) selection with crowding distance) Input: Expanded population sorted by stochastic nondomination F ¼ fF1 ; . . . ; Fl ; . . . ; FLðtÞ g Output: Population for the next generation Pðt þ 1Þ 1: l 1 2: while l 6 L(t) and Pðt þ 1Þ < P do 3: if jPðt þ 1Þj þ jFl j 6 P then 4: Pðt þ 1Þ Pðt þ 1Þ [ Fl 5: else E½fk ðxp ðtÞÞ  k ðxp ðtÞÞ , qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ; 8xp ðtÞ 2 Fl ; k ¼ 1; . . . ; K 6: E½f P 0 E½fk ðxp ðtÞÞ2 0 xp ðtÞ2Fl

7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21:

0, 8xp ðtÞ 2 Fl {initialize crowding distance} b(xp(t)) k 1 while k 6 K do Sort Fl by the kth objective (increasing order). The result is Fl ¼ fxð1Þ ðtÞ; . . . ; ðpÞ  k ðxðpþ1Þ ðtÞÞ P E½f  k ðxðpÞ ðtÞÞ. x ðtÞ; . . . ; xðjFl jÞ ðtÞg, where E½f (1) b(x (t)) 1 bðxðjFl jÞ ðtÞÞ 1    k ðxðpþ1Þ ðtÞÞ  E½f  k ðxðp1Þ ðtÞÞ ; p ¼ 2; . . . ; jFl j  1 bðxðpÞ ðtÞÞ bðxðpÞ ðtÞÞ þ E½f k k+1 end while Sort Fl by b(xp(t)) (decreasing order). The result is Fl ¼ fxð1Þ ðtÞ; . . . ; xðpÞ ðtÞ; . . . ; xðjFl jÞ ðtÞg, where b(x(p  1)(t)) P b(x(p)(t)) P ¼ P  jPðt þ 1Þj  Pðt þ 1Þ Pðt þ 1Þ [ fxð1Þ ðtÞ; . . . ; xðP Þ ðtÞg end if l l+1 end while

3. Computational experiments In this section we report our computational experience with the proposed evolutionary algorithm for project selection. In Section 3.1 we show how the algorithms parameters affect time and quality of the solution. Section 3.1.2 illustrates our experience with a variation of the crossover operator that favors the best parents. In Section 3.2 we compare our approach to the stochastic parameter space investigation method [35], a general purpose method for project selection. For comparison, we use a set of test problems also found in [35]. In Section 3.3 we illustrate by means of a variation of a problem found in the literature [23], how our proposed approach can be hooked to a Monte Carlo simulator to solve a realistic project selection example. The evolutionary algorithm was coded in the Java language [11] and executed using the Java HotSpot Virtual Machine version 1.4.2 from Sun Microsystems. All the experiments were conducted in a Dell machine with one processor Intel Pentium III running at 1GHz with 256MB of RAM.

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

879

3.1. Tuning the algorithms parameters 3.1.1. Injection rate and crossover probability We conducted an experiment to test how the algorithms parameters injection rate (pm) and probability of crossover (pc) affected the performance (computing time) and quality of the solutions. We designed the experiment by varying pm and pc from 0.1 to 0.9 in steps of 0.1. For each of the 81 experiments, we conducted 5 independent runs, fixing P = 20 and T = 40. Fig. 4 shows graphically that both injection rate and the crossover probability affect computing time. Furthermore, we adjusted a linear regression to the 81 samples with injection and crossover rates as the independent variables and computing time as the dependent variable. The model had a R2 = 0.98 and we rejected the hypothesis that both coefficients for the independent variables were 0. This fact provides us evidence to believe that by increasing any of the parameters, we have larger computing times. Both graphically and statistically, injection proved to be more expensive than crossover in terms of computing time. To measure the quality of the solution, we compared the stochastically nondominated front obtained with the proposed evolutionary algorithm (with a given combination of pm and pc) against the one obtained with stochastic PSI. We report the average fraction of solutions provided by the evolutionary algorithm in the aggregated front obtained by combining the solutions of both algorithms. In the remaining part of this section, we call this fraction absolute quality. Fig. 5 shows graphically that both injection rate and the crossover probability affect the absolute quality of the evolutionary algorithm. We adjusted a linear

Fig. 4. Average computational performance (in ms) varying injection rate (pm) and crossover probability (pc).

880

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

Fig. 5. Average absolute quality varying injection rate (pm) and crossover probability (pc).

regression to the 81 samples with injection and crossover rates as the independent variables and absolute quality as the dependent variable. We rejected the hypothesis that both coefficients for the independent variables were 0. This fact provides us evidence to assert that by increasing any of them, we have better absolute quality. Note that from Fig. 5 we see that if both injection and crossover are low (say 20%) the absolute quality is low (about 50%). By increasing injection or crossover (or both) we can see that there is an increasing trend on the absolute quality (reaching about 85% with high values of injection and crossover). Due to the fact that both parameters affect quality, we decided to set pm and pc to 0.5 in subsequent runs. This combination also provides a balanced tradeoff with respect to computing time (about 80,000 milliseconds in average). 3.1.2. Testing a modified crossover operator We implemented a variation of the crossover operator shown in Section 2.3.4 that favors the selection of the best parents. The parents are chosen by a tournament, based on stochastic nondomination, between two randomly chosen individuals from PðtÞ. In case there is no winner, we break ties in the tournament by picking one of the competing individuals randomly. The crossover between the two (winning) parents follows the mechanism outlined in Section 2.3.4. Henceforth we call this type of crossover SNDX (for stochastic nondomination), whereas we call RndX (for random) the traditional crossover that does not favor the best parents.

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

881

To compare SNDX against RndX, we picked the hardest problem from the set proposed in Ringuest and Graves [35], namely the one labeled ST03 in Section 3.2. Based on the results from Section 3.1.1 we set both pc and pm to 0.5. We also set P = 20 and T = 40. We executed 20 independent runs for each algorithm implementing SNDX and RndX. We first checked if there was a significant difference with regard to performance in terms of computing time. Table 1 summarizes the performance of SNDX and RndX measured in CPU time (in ms) for the 20 runs. We computed a t-test and could not reject the hypothesis that the average time for SNDX and RndX were equal. Furthermore, the 95% confidence interval for (lSNDX  lRndX) is [1331, 5806.2] (it contains 0). We further investigated if there was a significant difference between SNDX and RndX in terms of the quality of the solutions. We compared the Pareto fronts obtained with the algorithms implementing SNDX and RndX by calculating an aggregated front. We count how many of the stochastically nondominated solutions came from each of the procedures and calculated the fraction of the aggregated nondominated front provided by each of the algorithms (see Table 2). By using a paired t-test we found that it seems that there is no significant difference at a 95% confidence level. Furthermore, the 95% confidence interval for (lSNDX  lRndX) is [0.049, 0.2346] (it contains 0). Even though we did not find a significant difference between the crossover operators in terms of performance nor in quality, we found that the average ratio between tournaments with winners and the total number of crossovers in SNDX was only 58:1000. The reason being is that stochastic nondominance is not very discriminatory between individuals. It looks for dominance not only on the expected values of the criteria, but also on a given percentile (spread). Therefore, many of the tournaments in SNDX end up with ties, ultimately resembling RndX. We conducted a test with a larger population of 200 in a single run of ST03. For this single run the ratio increased to 111:1000. Furthermore, the fraction of nondominated vectors in the aggregated front was 53.2% and 46.8% for SNDX and RndX, respectively. On the other hand, the computing time was 764,640 and 784,007 millisecond for SNDF and RndX, respectively. This simple experiment shows a promising behavior for SNDX over RndX in runs with larger populations. 3.2. Comparison to the stochastic PSI In this section we summarize the results of the comparison between the proposed algorithm and the stochastic parameter space investigation (PSI) presented in Ringuest and Graves [35]. The five test problems used in this comparison are the ones presented in [35] based on ADBASE instance generator [46]. The problems objective function coefficients follow a triangular distribution with minimum of 0, maximum of 10, and random mode (within 0 and 10). Henceforth, we call these problems ST01, ST02, ST03, ST04, and

Table 1 Comparing SNDX vs. RndX: computing time (in milliseconds) Crossover type

Mean

Minimum

Maximum

Standard deviation

SNDX RndX

85,172.00 82,934.20

73,986.00 70,801.00

94,736.00 93,845.00

5018.69 6079.04

Table 2 Comparing SNDX vs. RndX: quality as a fraction of the aggregated stochastically nondominated front Crossover type

Mean (%)

Minimum (%)

Maximum (%)

Standard deviation (%)

SNDX RndX

45.4 54.6

4.8 31.8

68.2 95.2

15.2 15.2

882

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

ST05. Results for ST01 were shown in [26] with a preliminary version of the proposed evolutionary algorithm. For the evolutionary algorithm we executed 20 independent runs (different random seeds). The parameters were set as follows: population size P = 20, maximum number of generations T = 40, crossover rate pc = 0.50, and the injection rate pm = 0.50. For every feasible solution found by the algorithm, 10,000 replications of the objective functions were executed. For the stochastic PSI we executed 20 independent runs (different random seeds). Several trials of candidate solution vectors were generated until 100 points fell within the constraint set. For every feasible solution found by the algorithm, 10,000 replications of the objective functions were executed. PSI was coded in the SAS language [41]. We measure the computational performance in elapsed time. In total, 200 runs (=20 runs · 5 instances · 2 algorithms) were performed. Fig. 6 summarizes the comparison for problem ST01. The stochastic PSI took in average 61.6% longer to solve (see Fig. 6(a)). The evolutionary approach consistently generated 20 nondominated solutions in each run (except by one), while the stochastic PSI produced less and it was highly variable. On average, the sto21

150,000

19.9

140,000

Minimum Mean

130,000

Maximum Q3

19

123.320

120,000

110,000

100,000

90,000

Nondominated Solutions Generated

Milliseconds

Q1

80,000

(a)

7 5

4.95

Evolutionary

(b)

Q1 Minimum Mean Maximum Q3

60.00% 50.00% 40.00% 30.00% 20.00% 10.00% 0.00%

0.00%

Evolutionary

100.00%

90.00%

Q1 Minimum

80.00%

Mean

70.00%

Maximum Q3

60.00% 50.00% 40.00% 30.00% 20.00% 10.00% 0.00%

0.00%

Stochastic PSI Algorithm

Stochastic PSI Algorithm

100.00%

Fraction of solutions from combined pool of nondominated solutions

Fraction of nondominated solutions from total of generated solutions by algorithm

9

100.00%

100.00%

(c)

11

Stochastic PSI Algorithm

70.00%

Maximum Q3

1

Evolutionary

80.00%

13

3

76,305

70,000

90.00%

15

Q1 Minimum Mean

17

Evolutionary

(d)

Stochastic PSI Algorithm

Fig. 6. Comparison of the proposed method vs. stochastic PSI on instance ST01. (a) Computational performance, (b) number of solutions, (c) relative quality of the nondominated solutions and (d) absolute quality of the nondominated solutions.

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

883

chastic PSI generated five nondominated solutions, but in a given run it only generated 2 (and at most 10 in another) (see Fig. 6(b)). Among the nondominated solutions generated by the evolutionary approach, all were truly nondominated (relative quality of 100%); whereas, all the solutions generated by the stochastic PSI were dominated by those of the evolutionary approach. Fig. 7 summarizes the comparison for problem ST02. On average, the stochastic PSI took roughly twice the time needed by the evolutionary approach (see Fig. 7(a)). The evolutionary approach consistently generated 20 nondominated solutions in each run, while the stochastic PSI generated more. Nonetheless, by design the number of nondominated solutions generated by the evolutionary algorithms is bounded above by the population size of 20. As for ST01, the stochastic PSI was highly variable, it generated from 25 to 41 nondominated solutions in a given run, averaging 33 (see Fig. 7(b)). However, when comparing the quality of the solutions, 98.3% of those generated by the evolutionary algorithm were nondominated, while only 44.0% of those generated by the stochastic PSI were nondominated (see Fig. 7(c)). Among all the nondominated solutions, 57.8% and 42.2% where generated by the evolutionary and stochastic PSI, respectively (see Fig. 7(d)). Fig. 8 summarizes the comparison for problem ST03. Problem ST03 was particularly difficult for the stochastic PSI. The stochastic PSI took roughly 15 times longer than the evolutionary approach, averaging 42

210,000

170,000

40

Minimum Mean

38

Maximum Q3

Milliseconds

160,792

150,000

130,000

110,000

Nondominated Solutions Generated

190,000

Q1

90,000

Q1

36

Minimum Mean Maximum

34

Q3 33.4

32 30 28 26 24 22 20

20 73.364

70,000

18

Evolutionary

(a) 100.00%

70.00% 60.00%

Q1 Minimum Mean Maximum Q3

50.00% 43.96%

40.00% 30.00% 20.00% 10.00% 0.00%

90.00% 80.00% 70.00%

Q1 Minimum Mean Maximum Q3

60.00%

57.76%

50.00% 42.24%

40.00% 30.00% 20.00% 10.00% 0.00%

Evolutionary

(c)

Stochastic PSI Algorithm

100.00%

98.25%

90.00% 80.00%

Evolutionary

(b)

Fraction of solutions from combined pool of nondominated solutions

Fraction of nondominated solutions from total of generated solutions by algorithm

Stochastic PSI Algorithm

Stochastic PSI Algorithm

Evolutionary

(d)

Stochastic PSI Algorithm

Fig. 7. Comparison of the proposed method vs. stochastic PSI on instance ST02. (a) Computational performance, (b) number of solutions, (c) relative quality of the nondominated solutions and (d) absolute quality of the nondominated solutions.

884

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894 30 1,270,000

1,247,308

Minimum

Minimum Mean

Nondominated Solutions Generated

1,070,000

Milliseconds

Maximum 870,000

Q1

28

Q1

Q3

670,000

470,000

270,000

Mean

26

Maximum 24

Q3

22 20 18.75

18

17.4

16 14 12

81,117

70,000

10

Evolutionary

(a)

Stochastic PSI

Evolutionary

(b)

Algorithm

100.00% 92.38%

90.00% 80.00%

Q1 Minimum

70.00%

Mean Maximum Q3

60.00% 50.00% 40.00%

38.72%

30.00% 20.00% 10.00%

Fraction of solutions from combined pool of nondominated solutions

Fraction of nondominated solutions from total of generated solutions by algorithm

100.00%

0.00%

90.00% 80.00% 73.92%

70.00% 60.00% 50.00%

Q1 Minimum

40.00%

Mean

30.00%

Maximum Q3

26.08%

20.00% 10.00% 0.00%

Evolutionary

(c)

Stochastic PSI Algorithm

Stochastic PSI Algorithm

Evolutionary

(d)

Stochastic PSI Algorithm

Fig. 8. Comparison of the proposed method vs. stochastic PSI on instance ST03. (a) Computational performance, (b) number of solutions, (c) relative quality of the nondominated solutions and (d) absolute quality of the nondominated solutions.

20.8 minutes per run (see Fig. 8(a)). On average, the evolutionary approach generated 18.8 nondominated solutions per run, while the stochastic PSI generated about 17 (see Fig. 8(b)). However, with regard to quality of the solutions, 92.4% of those generated by the evolutionary algorithm were nondominated, while only 38.7% for the stochastic PSI (see Fig. 8(c)). Among all the nondominated solutions, 73.9% and 26.1% were produced by the evolutionary and stochastic PSI, respectively (see Fig. 8(d)). On one hand, we confirmed what was found by Ringuest and Graves [35], problem ST03 was particularly difficult for the stochastic PSI due to the fact that the feasible region filled an extremely small portion of the hypercube formed by the variable bounds. On the other hand, the evolutionary algorithm proved to be robust and solved it in roughly the same time as any other instance. Fig. 9 summarizes the comparison for problem ST04. The stochastic PSI took roughly 2.8 times longer than the evolutionary approach (see Fig. 9(a)). On average, the evolutionary approach generated twice as many nondominated solutions (see Fig. 9(b)). With regard to quality, 99.8% and only 10.4% of the nondominated solutions generated by the evolutionary algorithm and the stochastic PSI were nondominated, respectively (see Fig. 9(c)). Among all the nondominated solutions, 96.5% and just 3.5% were produced by the evolutionary and stochastic PSI, respectively (see Fig. 9(d)).

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894 270,000

20

20

230,000

Minimum Mean

226,446

Maximum Milliseconds

210,000

Q3

190,000 170,000 150,000 130,000 110,000

Nondominated Solutions Generated

Q1 250,000

885

Q1 Minimum

18

Mean 16

Maximum Q3

14 12 10

9.65

8 6

90,000 81,117

70,000

4

Evolutionary

(a)

80.00%

Q1 Minimum Mean

70.00% 60.00%

Maximum Q3

50.00% 40.00% 30.00% 20.00% 10.40%

10.00%

96.54%

90.00% 80.00% 70.00% 60.00%

Q1 Minimum Mean Maximum Q3

50.00% 40.00% 30.00% 20.00% 10.00% 3.46%

0.00%

0.00%

Evolutionary

(c)

Stochastic PSI Algorithm

100.00%

99.75%

Fraction of solutions from combined pool of nondominated solutions

90.00%

Evolutionary

(b)

Algorithm 100.00%

Fraction of nondominated solutions from total of generated solutions by algorithm

Stochastic PSI

Stochastic PSI Algorithm

Evolutionary

(d)

Stochastic PSI Algorithm

Fig. 9. Comparison of the proposed method vs. stochastic PSI on instance ST04. (a) Computational performance, (b) number of solutions, (c) relative quality of the nondominated solutions and (d) absolute quality of the nondominated solutions.

Fig. 10 summarizes the comparison for problem ST05. The stochastic PSI took roughly twice as much time as required by the evolutionary approach (see Fig. 10(a)). The evolutionary approach consistently generated 20 nondominated solutions (except once), while the stochastic PSI averaged only 14, with a high variance (see Fig. 10(b)). With regard to quality, 99.8% and just 7.3% of the nondominated solutions generated by the evolutionary algorithm and the stochastic PSI were nondominated, respectively (see Fig. 10(c)). Among all the nondominated solutions, 96.3% were contributed by the evolutionary algorithm (see Fig. 10(d)). 3.3. An R&D portfolio optimization example using Monte Carlo simulation Khorramshahgol and Gousty [23] proposed a goal programming model in which Delphic inquiry was used for R&D project optimization. Khorramshahgol and Gousty presented a model for a four-project R&D selection problem with two objectives: profit and market share. They use a five-year planning horizon, a discount rate of 10%, and expected cash flows. Projects can be partially funded, but just the equivalent of one full project is to be funded. Investment, returns and market share are proportional to the level of funding. Using a Delphic technique, they established a market share goal of 30% and a profit goal of $60 million.

886

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894 22

Milliseconds

170,000

Q1 Minimum Mean

20

Maximum Q3

165,828

150,000

130,000

110,000

90,000

Nondominated Solutions Generated

190,000

19.95

Q1 Minimum

18

Mean Maximum Q3

16 14

13.9

12 10 8 6

78,579

70,000

4

Evolutionary

Stochastic PSI

(a)

Stochastic PSI Algorithm

100.00%

99.75%

96.27%

90.00% 80.00% 70.00%

Q1 Minimum Mean Maximum Q3

60.00% 50.00% 40.00% 30.00% 20.00% 10.00%

7.28%

Fraction of solutions from combined pool of nondominated solutions

Fraction of nondominated solutions from total of generated solutions by algorithm

100.00%

Evolutionary

(b)

Algorithm

90.00% 80.00% 70.00% 60.00%

Q1 Minimum Mean Maximum Q3

50.00% 40.00% 30.00% 20.00% 10.00% 3.73%

0.00%

0.00%

Evolutionary

(c)

Evolutionary

Stochastic PSI

Stochastic PSI

(d)

Algorithm

Algorithm

Fig. 10. Comparison of the proposed method vs. stochastic PSI on instance ST05. (a) Computational performance, (b) number of solutions, (c) relative quality of the nondominated solutions and (d) absolute quality of the nondominated solutions. Table 3 Cash flows for R&D project selection example under uncertainty Year

Cash flows for project

Budget

A

B

C

D

0 1 2 3 4

30 Tria.(10, 20, 30) Expo(10)   Erlang 15 ; 2 0

30 Tria.(0, 10, 30) Expo(10) 1  Erlang 10 ;2 Erlang(1, 5)

0 30 Expo(30)   Erlang 15 ; 2 Erlang(1, 5)

0 30 Expo(10) 2  Erlang 35 ;2 Erlang(1,5)

30 25 0 0 0

Market share (in %)

Tria.(10, 20, 30)

Tria.(0, 30, 40)

Tria.(0, 25, 35)

Tria.(25, 35, 45)

–

In Khorramshahgol and Goustys original model cash flows are assumed to be expected values. To illustrate how our proposed method can deal with uncertainty, we remove this restrictive assumption and model cash flows and market shares as random variables. Table 3 shows the cash flows (in millions), budget, and market share for each project. The cash flows are random variables following triangular, exponential, and

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

887

Erlang distributions. In the table, a triangular random variable with minimum value a, maximum value c, and mode b is noted by Tria.(a, b, c); an exponential random variable with mean 1k is noted by Expoð1kÞ; and an Erlang random variable with scale parameter a and shape parameter k (with expected value ak) is noted by Erlang(a, k). This example was carefully designed so that the central moments of the random parameters share the same values of their deterministic counterparts in Khorramshahgol and Goustys original model. Let xi be the fraction of funding in project i (i 2 I ¼ fA; B; C; Dg). Let r be the discount rate. Let f~ ij be P

f~ ij ~ i be the (random) net the (random) cash flow of project i in year j (j = 0, . . . , 4). Let ~vi ¼ 4j¼0 ð1þrÞ and m j present value and (random) market share if project i is fully funded, respectively. Let bj be the budget avail~ i , and bj are shown in Table 3. The project selection model under able for year j. The parameters for f~ ij , m uncertainty follows:   ~ B ~ max V~ ; M; ð3Þ subject to X b1 ~vi  xi ; V~ ¼ b0 þ þ ð4Þ 1 þ r i2I X ~ ¼ ~ i  xi ; M ð5Þ m i2I

~ ¼ b1 þ B

X

f~ i1  xi ;

ð6Þ

i2I

b0 þ f~ A0 xA þ f~ A0 xB P 0; X xi ¼ 1;

ð7Þ ð8Þ

i2I

0 6 xi 6 1;

ð9Þ

where the first, second, and third objectives in (3) maximize net present value, market share, and the budget left at the end of year 1, respectively. The constraint (7) says that the initial budget minus the cash flows for projects A and B, should be nonnegative. Constraint (8) says that we can only fund at most the equivalent of one full project. Finally, constraints (9) allow partial funding for each project. We embedded a Monte Carlo simulator (written in Java) into the proposed evolutionary algorithm for project evaluation. The parameters for the evolutionary algorithm were set as follows: population size P = 20, maximum number of generations Nmax = 40, crossover rate pc = 0.35, and the injection rate pm = 0.20. For every feasible solution, 1000 replications of the simulator were executed. We used a discount rate of 10% and data from Table 3. We executed 20 independent runs of the experiment. Fig. 11 shows the expected values of the nondominated stochastic feasible solutions from all runs. From the nondominated stochastic solutions generated by the algorithm, in average, 86.8% were feasible (with a standard deviation of 2.9%). Although this fraction is high (and not variable), it is not possible to guarantee feasibility of all the Pareto solutions found by the algorithm. This is due to the fact that the budget left at the end of year 1 is a constraint with stochastic cash flows, so it could not be modelled as a (hard) linear constraint, but as a third objective to be maximized. For the original problem, Khorramshahgol and Gousty solve a goal program and obtain xA = 0.19, xB = 0.44, xC = 0, and xD = 0.37, with a profit V = 60 and M = 30. Among all stochastically nondominated feasible solutions shown in Fig. 11, we found many solutions that are comparable to that of Khorramshahgol and Gousty (see Table 4). Note that all of these solutions (in objective space) meet the ~ P 30, E½B ~ P 0, and T B~ P 0 (where P ðB ~ > T B~ Þ ¼ 0:90). In other words, following criteria: E½V~  P 60, E½M all the stochastically nondominated feasible solutions found with the evolutionary algorithm provide

888

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

70 60 50

VPN

40 30 20 10 0 35 30 M

50

ts

ke

ar

25 re

ha

40 20

30 20 15

10 0

left Budget

(end of

year 1)

Fig. 11. Projected nondominated feasible solutions for R&D project selection example under uncertainty.

enough evidence to show that they are able to meet the budget constraints as well as the target profit and market share of $60 million and 30%, respectively. Furthermore, note that as with Khorramshahgol and Goustys original solution, all the Pareto solutions shown in Table 5 do not find attractive to fund project C. The column ID in Tables 4 and 5 is an identifier composed of the run number (1–20) and the sequence within the Pareto optimal front (1 to P = 20). In Table 6, we show for each Pareto optimal front obtained with the algorithm, the average allocation for each project (in solution space). Similarly, in Table 7 we calculated the average nondominated criteria vectors for each Pareto optimal front obtained with the algorithm. These averages are given to provide an indication of the magnitude of the allocations and objectives and the methods robustness (low variation). However, caution should be exercised with these results because our method is based on posterior articulation of preferences and its main objective is to unveil the (unknown) efficient frontier and not to provide aggregated measures.

4. Concluding remarks This paper presents a new evolutionary method for solving project selection problems with partially funded projects, multiple objectives, project interdependencies (in the objectives), and a linear structure for resource constraints. The algorithm allows uncertain elements in the objectives to be modelled as random coefficients in complex (nonlinear) expressions. The algorithm also allows project evaluation being treated as a black box. In the latter case, the method simulates the project and evaluates its performance

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

889

Table 4 Nondominated feasible solutions for R&D project selection example under uncertainty ID

E½V~  (NPV)

~ E½M (Market share)

~ E½B (Budget)

T V~ ðP ðV~ > T V~ Þ ¼ 0:90Þ

~ > T M~ Þ ¼ 0:90Þ T M~ ðP ðM

~ > T B~ Þ ¼ 0:90Þ T B~ ðP ðB

6–3 11–2 13–14 3–9 7–20 18–2 5–18 17–17 14–2 16–3 5–19 11–19 12–2 20–14 6–14 19–13 19–8 15–14 2–2 3–6 19–14 7–7 18–19 9–1 7–16 20–13 20–17 4–11 18–8 1–16 10–15 4–8 11–12 13–19 17–7 16–1 3–7 2–1 8–18 10–6 12–12 4–6 13–18 6–15 18–13 1–6

64.27 63.91 63.60 63.46 63.33 63.19 63.11 62.92 62.83 62.78 62.76 62.73 62.72 62.68 62.65 62.62 62.59 62.59 62.55 62.41 62.39 62.35 62.30 62.29 62.22 62.21 62.06 62.04 62.03 62.02 61.99 61.96 61.93 61.89 61.80 61.61 61.47 61.47 61.39 61.35 61.32 61.32 61.28 61.19 60.82 60.35

32.97 32.52 32.29 32.86 32.85 32.16 30.46 31.37 31.32 32.24 33.21 31.54 31.99 31.59 32.69 32.18 32.94 32.47 32.50 30.46 30.40 30.86 32.77 30.86 30.10 30.66 32.30 32.35 31.00 31.30 32.89 30.77 30.55 31.32 30.38 30.90 31.66 30.86 30.11 30.10 30.51 31.12 30.41 30.65 30.09 30.84

1.92 3.59 3.63 2.87 3.97 5.27 10.95 5.62 8.87 4.04 1.30 6.72 5.04 6.27 2.86 5.50 1.79 4.01 2.99 11.44 11.62 9.17 2.78 10.28 12.02 10.32 3.88 3.58 8.61 7.17 2.36 9.13 9.80 7.52 10.63 8.42 6.64 9.75 11.88 11.12 9.59 7.77 10.64 9.34 11.48 9.52

45.58 43.71 44.64 44.58 45.09 44.91 46.75 46.42 46.62 44.62 44.10 45.06 44.84 44.76 44.98 45.70 43.18 44.61 44.44 47.83 47.63 45.85 44.72 46.99 46.10 46.42 44.43 44.15 45.10 45.53 43.54 45.79 46.67 46.14 47.15 44.26 44.91 45.92 46.96 46.21 45.06 44.85 46.43 45.46 45.65 45.39

28.25 27.94 27.62 27.86 28.04 27.67 26.53 26.72 26.54 27.55 28.43 27.57 27.60 27.43 27.87 27.60 28.23 27.78 27.81 26.29 25.99 26.85 28.22 25.85 26.03 26.62 27.89 27.83 26.95 26.87 28.34 26.84 26.54 27.09 26.43 26.70 27.34 27.01 25.99 26.22 26.45 26.71 26.33 26.75 26.25 26.88

1.23 2.75 2.79 1.43 2.32 4.20 9.43 3.68 6.61 3.11 0.73 5.45 4.10 5.19 2.12 4.43 1.12 2.98 2.19 9.55 9.83 7.68 2.04 7.78 10.28 8.74 3.00 2.74 7.08 5.87 1.66 7.73 8.46 6.36 9.20 6.89 5.54 8.36 10.16 9.47 7.95 6.37 8.96 7.86 9.70 8.11

measures (objectives) subject to random elements. The method is based on posterior articulation of preferences and is able to approximate the efficient frontier composed of stochastically nondominated solutions.

890

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

Table 5 Pareto feasible solutions for R&D project selection example under uncertainty ID

xA

xB

xC

xD

6–3 11–2 13–14 3–9 7–20 18–2 5–18 17–17 14–2 16–3 5–19 11–19 12–2 20–14 6–14 19–13 19–8 15–14 2–2 3–6 19–14 7–7 18–19 9–1 7–16 20–13 20–17 4–11 18–8 1–16 10–15 4–8 11–12 13–19 17–7 16–1 3–7 2–1 8–18 10–6 12–12 4–6 13–18 6–15 18–13 1–6

0.08 0.14 0.14 0.00 0.00 0.11 0.24 0.01 0.03 0.18 0.10 0.22 0.17 0.20 0.10 0.12 0.11 0.07 0.14 0.13 0.14 0.26 0.14 0.04 0.32 0.16 0.15 0.14 0.27 0.22 0.11 0.24 0.20 0.15 0.22 0.27 0.19 0.23 0.19 0.31 0.28 0.25 0.16 0.24 0.32 0.24

0.07 0.04 0.04 0.18 0.21 0.11 0.09 0.24 0.28 0.00 0.03 0.02 0.04 0.03 0.06 0.10 0.03 0.13 0.02 0.22 0.22 0.03 0.02 0.31 0.02 0.17 0.03 0.03 0.00 0.03 0.04 0.05 0.11 0.12 0.11 0.00 0.05 0.08 0.16 0.02 0.01 0.00 0.17 0.05 0.01 0.06

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.04 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.85 0.82 0.82 0.82 0.79 0.78 0.67 0.71 0.69 0.82 0.87 0.76 0.79 0.76 0.83 0.78 0.86 0.80 0.84 0.64 0.64 0.71 0.84 0.65 0.66 0.66 0.81 0.82 0.73 0.75 0.85 0.71 0.69 0.73 0.66 0.73 0.76 0.69 0.64 0.67 0.71 0.74 0.66 0.70 0.67 0.70

We compared the proposed method with the stochastic parameter space investigation method (PSI) proposed by Ringuest and Graves [35]. Compared to the stochastic PSI, the results showed that the method is faster and more robust, provides higher quality of the nondominated solutions, and is also able to consistently generate a controlled number of solutions that approximate the efficient frontier. In problem

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

891

Table 6 Average allocation for R&D project selection example under uncertainty Run

xA

xB

xC

xD

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

0.45 0.47 0.45 0.46 0.52 0.43 0.46 0.41 0.44 0.53 0.38 0.51 0.39 0.37 0.40 0.46 0.46 0.42 0.44 0.47

0.13 0.14 0.14 0.11 0.13 0.19 0.16 0.23 0.24 0.12 0.17 0.10 0.24 0.27 0.22 0.12 0.20 0.14 0.18 0.13

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.01 0.00 0.00 0.00

0.42 0.39 0.41 0.44 0.34 0.38 0.37 0.36 0.33 0.34 0.44 0.39 0.37 0.35 0.38 0.41 0.33 0.44 0.38 0.40

Average Std. dev.

0.45 0.04

0.17 0.05

0.00 0.00

0.38 0.04

Table 7 Average nondominated criteria vectors for R&D project selection example under uncertainty Run

E½V~  (NPV)

~ E½M (Market share)

~ E½B (Budget)

T V~ ðP ðV~ > T V~ Þ ¼ 0:90Þ

~ > T M~ Þ ¼ 0:90Þ T M~ ðP ðM

~ > T B~ Þ ¼ 0:90Þ T B~ ðP ðB

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

21.44 21.98 22.09 22.93 21.38 20.51 21.45 20.75 18.46 20.75 22.19 21.68 21.64 19.58 22.09 23.10 21.47 25.05 19.34 17.64

21.16 21.61 20.04 17.80 22.67 22.43 20.86 22.27 26.52 24.99 18.56 21.76 18.90 21.32 20.91 18.78 21.63 11.37 28.90 33.35

0.47 0.47 0.44 0.42 0.54 0.44 0.45 0.38 0.52 0.57 0.40 0.54 0.36 0.22 0.35 0.46 0.43 0.30 0.55 0.70

0.14 0.14 0.14 0.11 0.09 0.20 0.16 0.27 0.21 0.11 0.15 0.06 0.22 0.45 0.26 0.09 0.19 0.11 0.20 0.08

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00

0.39 0.39 0.42 0.47 0.37 0.36 0.39 0.35 0.26 0.32 0.45 0.39 0.42 0.33 0.39 0.45 0.36 0.59 0.26 0.22

Average Std. dev.

21.28 1.66

21.79 4.45

0.45 0.11

0.17 0.09

0.00 0.00

0.38 0.08

892

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

instances where the feasible region is relatively small compared to the enclosing hypercube formed by the bounds on the variables, the stochastic PSI takes a very long time to solve the project selection problem. In contrast, the evolutionary algorithm takes roughly the same time to solve the problem regardless of the size of the feasible region. Nonetheless, the stochastic PSI method has several advantages. Among them, it can be applied to a larger class of project selection problems and it is very easy to implement. We illustrated the proposed method with a variation of the Khorramshahgol and Gousty [23] model for R&D project optimization. We relaxed deterministic assumptions and added stochastic components to the model, namely, cash flows and market shares were modelled by random variables. The method was able to approximate the efficient frontier with stochastically nondominated solutions. The evolutionary algorithm has been designed to solve larger problems than the ones handled by the stochastic PSI. The language in which the algorithm was developed (Java) allows for portability and extensibility. The format of the input files allows for solving larger problems expressing the linear constraints in a compact MPS-like format [30]. An example of the input format is given in [26]. Furthermore, it is worth mentioning that the proposed algorithm outlined in this paper is able to tackle a class of stochastic multiobjective optimization problems with linear constraints, where the project selection problem described herein is a particular case. Aristeguieta et al. [2,3] have used the proposed evolutionary algorithm in the context of project optimization in the oil and gas industry. In their work, they use the proposed evolutionary algorithm on an example of an exploration and production company. This company is presented with a set of eight projects and its goal is to maximize the expected net present value and the cumulative production subject to a budget constraint. They assume the technical and economic teams have previously evaluated the projects using Monte Carlo simulation. They compare the algorithm with classical approaches such as the rank-andcut method and a mean–variance approach. The authors conclude that the evolutionary approach shown in this paper is a very promising technique to optimize oil and gas portfolios due to the fact that unveils solutions hidden to the competing methods. We are currently conducting joint research to include project dependencies and scheduling in this context. We are currently conducting research on an extension of the algorithm for solving go-no-go projects. Note that binary decisions with regard to funding, destroy the convexity of the solution space. Another promising line of research is the design and implementation of alternative heuristics based on simulated annealing and tabu search for the project selection problem presented herein and its outlined variants.

Acknowledgement This study was partially funded by the grant CIFI-2057 of the Research Fund for Assistant Professors at Universidad de los Andes.

References [1] D.A. Aaker, T.T. Tyebjee, A model for the selection of interdependent R&D projects, IEEE Transactions on Engineering Management 25 (2) (1978) 30–36. [2] O. Aristeguieta, S. Begg, R. Bratvold, Improving oil and gas business performance through project portfolio optimisation, In: 6th International Conference on Optimization: Techniques and Applications (ICOTA), Ballarat, Australia, 2004, Oral presentation. [3] O. Aristeguieta, S. Begg, R. Bratvold, Improving oil and gas business performance through project portfolio optimisation, Working Paper, University of Adelaide, Australian School of Petroleum, 2004. [4] T. Back, F. Hoffmeister, Extended selection mechanisms in genetic algorithms, in: R. Belew, L. Booker (Eds.), Proceedings of the Fourth International Conference on Genetic Algorithms, Morgan Kaufmann Publishers, San Mateo, CA, 1991, pp. 92– 99.

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

893

[5] D. Beasley, D.R. Bull, R.R. Martin, An overview of genetic algorithms: Part 1. Fundamental, University Computing 15 (2) (1993) 58–69. [6] D. Beasley, D.R. Bull, R.R. Martin, An overview of genetic algorithms: Part 2. Research topics, University Computing 15 (4) (1993) 170–181. [7] J.E. Beasley, Population heuristics, in: P.M. Pardalos, M.G.C. Resende (Eds.), Handbook of Applied Optimization, Oxford University Press, 2000. [8] C.A. Coello, A short tutorial on evolutionary multiobjective optimization, in: E. Zitzler, K. Deb, L. Thiele, C.A. Coello, D. Corne (Eds.), Evolutionary Multi-Criterion Optimization, Springer, Washington, DC, 2001, pp. 21–40. [9] C.A. Coello, G.A. Lamont, D.A. Van Veldhuizen, Evolutionary Algorithms for Solving Multi-Objective Problems, Kluwer, New York, 2002. [10] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multi-objective genetic algorithm: NSGA-II, IEEE Transactions on Evolutionary Computation 6 (2) (2002) 182–197. [11] B. Eckel, Thinking in Java, third ed., Prentice-Hall, 2002. [12] S.C. Fang, S. Puthenpura, Linear Programming and Extensions, Prentice-Hall, Englewood Cliffs, NJ, 1993. [13] L.J. Fogel, A.J. Owens, M.J. Walsh, Artificial Intelligence Through Simulated Evolution, John Wiley, New York, 1966. [14] C.M. Fonseca, P.J. Fleming, Genetic algorithms for multi-objective optimization: Formulation, discussion and generalization, in: S. Forrest (Ed.), Proceedings of the Fifth International Conference on Genetic Algorithms, Morgan Kaufmann, San Mateo, CA, 1993, pp. 416–423. [15] G.E. Fox, N.R. Baker, J.L. Bryant, Economic models for R and D project selection in the presence of project interactions, Management Science 30 (7) (1984) 890–902. [16] D.E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley, Reading, MA, 1989. [17] J. Goldwerger, J. Paroush, Capital budgeting of interdependent projects: Activity analysis and approach, Management Science 23 (11) (1977) 1242–1246. [18] S.B. Graves, J.L. Ringuest, Models & Methods for Project Selection: Concepts from Management Science, Finance and Information Technology, Kluwer Academic Publishers, 2003. [19] K. Heidenberger, C. Stummer, Research and development project selection and resource allocation: A review of quantitative modelling approaches, International Journal of Management Reviews 1 (2) (1999) 197–224. [20] J.H. Holland, Adaptation in Natural and Artificial Systems, University of Michigan Press, Ann Arbor, 1975. [21] J. Horn, N. Nafploitis, D.E. Goldberg, A niched Pareto genetic algorithm for multi-objective optimization, in: Z. Michalewicz (Ed.), Proceedings of the First IEEE Conference on Evolutionary Computation, IEEE Service Center, Piscataway, NJ, 1994, pp. 82–87. [22] A.J. Keown, B.W. Taylor, C.P. Duncan, Allocation of research and development funds: A zero-one goal programming approach, Omega 7 (4) (1979) 345–351. [23] R. Khorramshahgol, Y. Gousty, Delphic goal programming (DGP): A multiobjective cost/benefit approach to R&D portfolio analysis, IEEE Transactions on Engineering Management 33 (3) (1986) 172–175. [24] J.D. Knowles, D.W. Corne, Approximating the nondominated front using the Pareto archived evolution strategy, Evolutionary Computation 8 (2) (2000) 149–172. [25] H.C. Lucas, J.R. Moore Jr., A multiple-criterion scoring approach to information system project selection, INFOR 14 (1) (1976) 1–12. [26] A.L. Medaglia, Models and Methods for Project Selection: Concepts from Management Science, Finance, and Information Technology, Kluwer Academic Publishers, Boston, 2003, Chapter: An evolutionary algorithm for project selection problems based on stochastic multiobjective linearly constrained optimization, pp. 163–189. [27] Z. Michalewicz, Genetic Algorithms + Data Structures = Evolution Programs, third ed., Springer, 1996. [28] K. Muralidhar, R. Santhanam, M. Schniederjans, An optimization model for information system project selection, Journal of Management Science and Policy Analysis 6 (1) (1988) 53–62. [29] K. Muralidhar, R. Santhanam, R.L. Wilson, Using the analytic hierarchy process for information system project selection, Information and Management 18 (1) (1990) 87–95. [30] B.A. Murtagh, Advanced Linear Programming: Computation and Practice, McGraw-Hill International Book Co., 1981. [31] G.L. Nemhauser, Z. Ullmann, Discrete dynamic programming and capital allocation, Management Science 15 (9) (1969) 494–505. [32] I. Rechenberg, Evolutionsstrategie—Optimierung technischer Systeme nach Prinzipien der biologischen Evolution, FrommannHolzboog Verlag, Stuttgart, Germany, 1973. [33] J.L. Ringuest, S.B. Graves, The linear multi-objective R&D project selection problem, IEEE Transactions on Engineering Management 36 (1) (1989) 54–57. [34] J.L. Ringuest, S.B. Graves, The linear R&D project selection problem: An alternative to net present value, IEEE Transactions on Engineering Management 37 (2) (1990) 143–146. [35] J.L. Ringuest, S.B. Graves, A sampling-based method for generating nondominated solutions in stochastic MOMP problems, European Journal of Operational Research 126 (2000) 651–661.

894

A.L. Medaglia et al. / European Journal of Operational Research 179 (2007) 869–894

[36] J.L. Ringuest, T.R. Gulledge Jr., Interactive multiobjective complex search, European Journal of Operational Research 19 (1985) 362–371. [37] R. Santhanam, G.J. Kyparisis, A multiple criteria decision model for information system project selection, Computers and Operations Research 22 (8) (1995) 807–818. [38] R. Santhanam, G.J. Kyparisis, A decision model for interdependent information system project selection, European Journal of Operational Research 89 (1996) 380–399. [39] R. Santhanam, K. Muralidhar, M. Schniederjans, A zero-one goal programming approach for information system project selection, Omega—International Journal of Management Science 17 (6) (1989) 583–593. [40] R. Santhanam, M. Schniederjans, A model formulation system for information system project selection, Computers and Operations Research 20 (7) (1995) 757–767. [41] SAS Institute Inc., SAS Language: Reference, Version 6, first ed., SAS Institute Inc., Cary, NC, 1990. [42] J.D. Schaffer, Multiple objective optimization with vector evaluated genetic algorithms, in: J.J. Grefenstette (Ed.), Proceedings of the 1st International Conference on Genetic Algorithms, Erlbaum, Hillsdale, NJ, 1985, pp. 93–100. [43] M. Schniederjans, R. Santhanam, A multi-objective constrained resource information system project selection method, European Journal of Operational Research 70 (2) (1993) 244–253. [44] H.P. Schwefel, Numerical Optimization of Computer Models, Wiley, Chichester, 1981. [45] N. Srinivas, K. Deb, Multiobjective optimization using nondominated sorting in genetic algorithms, Evolutionary Computation 2 (3) (1994) 221–248. [46] R.E. Steuer, ADBASE: A multiple objective linear solver, Technical report, Terry College of Business, University of Georgia, Athens, Georgia, 2000. [47] R.E. Steuer, M. Sun, The parameter space investigation method for multiple objective nonlinear programming: A computational investigation, Operations Research 43 (4) (1995) 641–648. [48] H.M. Weingartner, Mathematical Programming and the Analysis of Capital Budgeting Problems, Prentice-Hall, Englewood Cliffs, NJ, 1963. [49] H.M. Weingartner, Capital budgeting of interrelated projects: Survey and synthesis, Management Science 12 (7) (1966) 485–516. [50] E. Zitzler, L. Thiele, Multiobjective evolutionary algorithms: A comparative case study and the strength Pareto approach, IEEE Transactions on Evolutionary Computation 3 (4) (1999) 257–271.