Supply Chain RedesignâMultimodal Optimization Using a Hybrid

11094

Ind. Eng. Chem. Res. 2009, 48, 11094–11107

Supply Chain RedesignsMultimodal Optimization Using a Hybrid Evolutionary Algorithm P. K. Naraharisetti,† I. A. Karimi,*,‡ and Rajagopalan Srinivasan*,†,‡ Institute of Chemical and Engineering Sciences, ASTAR (Agency for Science, Technology & Research), 1 Pesek Road, Jurong Island, Singapore 627833, and Department of Chemical & Biomolecular Engineering, National UniVersity of Singapore, 4 Engineering DriVe 4, Singapore 117576

Supply chain redesign (SCR) involves decisions regarding the timings, amounts, and locations of the investment and disinvestment in facilities, production, material purchase, product sales, contracts, capital-raising loans and bonds, etc. such that the profit is maximized. SCR is a heavily constrained problem; hence as the problem size increases, the MILP formulations (Naraharisetti, P. K.; Karimi, I. A.; Srinivasan, R. Supply Chain Redesign through Optimal Asset Management and Capital Budgeting. Comput. Chem. Eng. 2008, 32, 3153-3169) become increasingly difficult to solve. In addition, MILP solvers typically give only one solution, while multiple optimal solutions may be desirable in practice. Hence, an alternative optimization technique is warranted. In this work, we propose a hybrid MILP-evolutionary algorithm strategy for supply chain redesign and present progress on three fronts: (a) a novel reformulation of the MILP in which most decision variables are unconstrained and the rest can be easily repaired to satisfy constraints, (b) a single-objective hybrid optimization algorithm that uses an evolutionary search and reaches 97% of the objective value attained by CPlex 9.0 on a small example, while outperforming CPlex 9.0 on a large SCR problem, and (c) a multimodal algorithm that identifies multiple supply chain networks with 90-95% of the objective value obtained by CPlex 9.0. Finally, we analyze the effect of uncertainty on each supply chain network identified by our multimodal algorithm. 1. Introduction With the advent of globalization, new markets are opening up in various parts of the world and organizations are venturing into these markets to exploit new opportunities and maximize shareholder value. This requires efficient management of their supply chains and assets. Typical supply chain assets include various production and inventory facilities, raw material and product inventories, technological expertise, and financial assets such as capital from loans and bonds and contracts for supply and sales. To stay profitable, the organizations continually invest in high-profit assets and disinvest the low-profit ones. While making these supply chain management decisions, they must also consider various issues such as logistics costs,2,3 regulatory factors,4,5 and risk6 due to uncertainties. Hence, a methodology that considers these issues and gives a plan to maximize shareholder value is of utmost importance. 1.1. Chemical Supply Chain Redesign. Considerable literature exists on the specific issues of new facility location and/ or technology upgrade of an existing facility. A comprehensive review on the design and operation of supply chains with emphasis on financial strategies and considering several issues such as international tax planning, financing of investments, transfer pricing and uncertainties, etc. was presented by McDonald and Reklaitis.7 It also describes various modeling and solution strategies such as simulation-optimization, real options, and multistage models, among others. Several capacity expansion models exist in the literature. Bjorkqvist and Roslof8 * To whom correspondence should be addressed. (I.A.K.) E-mail: [email protected]. Tel.: +65 6516 6359. Fax: +65 6779 1936. (R.S.) E-mail: [email protected]. Tel.: +65 6516-8041. Fax: +65 6779-1936. † ASTAR. ‡ National University of Singapore.

assumed capacity increments in multiples of some discrete amount (defined a priori). Bok et al.9 addressed capacity expansion for the Korean petrochemical industry, and others10-12 addressed capacity expansion and production-distribution in the pharmaceutical industry. In addition, Oh and Karimi4 have incorporated the effect of regulatory factors in a strategic capacity planning model, which is relevant in this era of globalization. A multiobjective supply chain design model considering capacity additions and the design of chemical supply chains under demand uncertainty were presented by Guillen and others.13,14 Several researchers have considered the problem of production allocation and distribution. To date, Arntzen and others15 have probably presented the most comprehensive model to address various issues in the supply chain of the electronics industry. Production planning and distribution by grouping products into families and including transportation times in the formulation was considered by McDonald and Karimi.16 Wilkinson and others17 presented an integrated production and distribution model and solved it by using aggregated periods. A production allocation and distribution model, while considering economies of scale in material purchase/sales, was presented by Tsiakis and others.18 Gjendrum and others19 presented a nonlinear model to consider transfer pricing in the supply chain. Models for a global logistics system and, considering transfer pricing20 a model for production-distribution planning of multiproduct batch plants and the use of model predictive control strategy on a rolling horizon,21 an enterprise-wide production and distribution model considering various operational policies22 were also presented, while Oh and Karimi5 considered a production-distribution planning model considering duty

10.1021/ie9002574 CCC: $40.75  2009 American Chemical Society Published on Web 11/02/2009

Ind. Eng. Chem. Res., Vol. 48, No. 24, 2009

drawbacks, which are becoming increasingly important due to globalization and competition between organizations. On the basis of our literature review, a formulation that considers the issues of disinvestments (thus relocation) of facilities, technology upgrade, strategic contracts related to the sales and purchase of materials, and loans and bonds for raising capital is not present in the literature. In addition, the following simplifying assumptions are common (of course, exceptions do exist) in existing literature and need further attention: (1) Discrete amounts (multiples of a fixed amount) of investments are assumed instead of continuous amounts within some limits. (2) Lower limits on the capacity of a new construction and that of the expansion of an existing facility are assumed identical. (3) Existing facilities are not allowed to shutdown, so a maintenance shutdown cannot occur. (4) Lower limits on the capacity utilizations of production facilities are assumed to be zero (rather than nonzero) in order to avoid infeasibilities in the absence of raw materials. We previously presented a supply chain redesign model that addresses the above issues.1 Briefly, we considered issues such as investment, disinvestment, relocation, regulatory factors, transportation cost, contracts for strategic material supply, and loans and bonds for raising capital, among others. Our MILP model for a small academic example was implemented in GAMS 21.7 and solved using Cplex 9.0 to a gap of 9% in 24 h. Clearly, the computational effort is large even for a small example. Considering a typical problem in the industry, the MILP model may not even give a feasible solution in reasonable time. Hence, we need alternative solution techniques. 1.2. Multimodal Optimization. A nonlinear optimization problem with an objective f(x,y) ) g(x,y) + Rxβyγ (where f and g are some functions of variables x and y, and R, β, and γ are parameters) may exhibit multimodality. Similarly, a linear optimization problem with an objective f(x,y) ) g(x,y) + Rx + βy may also exhibit multimodality. For example, in a traveling salesman problem, multiple optimal routes may exist for visiting all locations. MILP solvers based on the branch and bound algorithm generally give a single solution, which is guaranteed to be optimal given sufficient solution time. Often, several solutions/plans are required instead of a single solution. Having a set of solutions gives the decision maker an opportunity to analyze the performance of each plan with respect to the uncertain parameters and hence a multimodal optimization algorithm that can identify multiple alternate optimal solutions is necessary. As the MILP size increases, the chances of not getting a good solution in the given time also increase. This is because a rigorous branch and bound algorithm attempts to guarantee optimality using a systematic search. It normally maintains a single best solution and discards others. Evolutionary algorithms (EAs) on the other hand generate random solutions and maintain a set of solutions. In addition, in contrast to an MILP solver, they also give the flexibility to initialize the search with a rough solution derived by the planner. The supply chain redesign (SCR) MILP model has thousands of variables. An MILP solver solves a series of linear programming (LP) problems by fixing some binary variables. However, the main feature is that it maintains solution feasibility throughout. In contrast, EAs are best suited for unconstrained problems. The various solutions encoded as chromosomes are generally created randomly and may easily violate constraints. Because of the large numbers

11095

of variables and constraints in a large MILP, it is challenging to generate chromosomes that always remain in the feasible solution space. However, it is important to note that the chromosome infeasibilities in the SCR problem are largely due to the continuous variables and the binary variables can easily be remodeled such that they create few or no infeasibilities. In this study, we reformulate the constraints involving some binary variables to obtain a formulation that is partially unconstrained and use repair to keep the rest of the binary variables in the feasible region. Once a feasible set of binary variables is obtained, we solve the reduced MILP to get the objective. A large variety of evolutionary algorithms exists in the literature, including genetic algorithm,23 differential evolution,24 tabu search,25,26 scatter search,27 ant colony optimization,28 and simulated annealing.29,30 The choice of the algorithm depends on the characteristics of the problem and there is no “single best algorithm”. Hence, we have developed a hybrid evolutionary algorithm customized to the SCR problem. It is our belief that our proposed hybrid algorithm would perform well with other problems as well. For implementing the hybrid algorithm, a GAMS-Matlab link31 is used. In addition to allowing the use of algorithms (either built in or user developed) from Matlab to be combined with GAMS, the link also combines the optimization capabilities of GAMS with the visualization capability of Matlab. A decomposition technique for a large scale stochastic programming problem was developed where the GAMS-Matlab link was used.14 In their study, the supply chain design decisions are made by the genetic algorithm in Matlab and a smaller stochastic programming problem for operational decisions is solved in GAMS. In our work, the evolutionary algorithm is developed in Matlab and the MILP is solved in GAMS. Finally, it should be noted that several other methodologies such as Lagrangean, Benders, and bi-level decomposition do exist for solving large MILPs, and would be worth applying to this problem. However, at present, we have chosen evolutionary algorithms for this study. 1.3. This Work. In this work, we present progress on three fronts: (a) a novel reformulation of the SCR MILP in which most decisions are unconstrained and the rest are easily repaired to satisfy constraints; (b) a single-objective hybrid optimization algorithm, which uses an evolutionary search and reaches 97% of the objective value attained by CPlex 9.0 for a small example, while outperforming CPlex 9.0 for a large SCR problem; (c) a multimodal algorithm that reaches 90-95% of the objective value obtained by CPlex 9.0, while identifying multiple supply chain networks. Finally, we analyze the effect of uncertainty on each supply-chain network identified by our multimodal algorithm. The overview of this work is presented in Figure 1. Constraint handling and reformulation are presented in section 2. The development of the evolutionary algorithm and the strategy for multimodal optimization are presented in section 3. Two examples, one on identifying multiple optimal solutions and another on a large supply chain, are presented in section 4. 2. Constraint Handling Considering the large number of constraints in our model, the crossover and mutation operations often generate infeasible strings. There are several methods of handling these

11096


Figure 1. (Top) Reformulation of the supply chain redesign problem into small MILP (sMILP) and chromosome for use in evolutionary algorithm; (bottom) the use of reformulation and multimodal hybrid evolutionary for supply chain redesign.

infeasibilities and some of them are to (a) reject infeasible strings, (b) penalize objective when infeasible strings are generated, (c) repair the infeasible strings, and (d) generate strings such that they are always feasible. The strategies of rejecting the infeasible strings or penalizing the objective

function perform poorly because the SCR problem is heavily constrained. Our strategy is to generate most of the decision variables in a string in the feasible region and repairing the infeasibilities in the others. This is partially similar to that presented elsewhere.32 In the literature,32 researchers have


handled the constraints by generating the decision variables in the feasible region. This is achieved by modifying the algorithm from case to case. In our strategy, the evaluation function varies from case to case and the algorithm remains the same. In this paper, we will elaborate on this remodeling and highlight the relative advantages and disadvantages of solving an MILP model directly in GAMS vis-a`-vis the hybrid algorithm. In addition, we present an algorithm for multimodal optimization. 2.1. MILP Formulation. First, consider the constraints involving the binary variables in the SCR problem, which are reformulated in this work. Detailed formulation can be found elsewhere,1 where the MILP uses a discrete-time representation with H periods of equal length in the horizon. Capacity changes in the supply chain can be modeled by including all the facilities in the networksboth existing and potential facilities. We use binary variables EYet, EXet, TUet and 0-1 continuous variables EEYet, ESYet to model the supply chain redesign. The variables are written to denote the existence of a facility at time t (EEYet), capacity addition at time t (EYet), disinvestment at time t (EXet), existence in the past followed by a disinvestment of a facility (ESYet) and finally, technology upgrade (TUet). Hence, 3H binary and 2H (0-1) continuous variables are required for each facility to model SCR. Our primary assumptions are that a capacity expansion takes time and only one expansion can be planned at any given time. Hence, a second expansion can start only when the first is completed. A facility that is disinvested cannot be purchased back. The number of times a facility may expand in the planning horizon is limited and is known a priori. Although no capacity may be added, a technology upgrade is also considered as an addition with reference to this limit on expansions as work needs to be done. Because of the time required for construction of a facility, that is construction lead time (CLTe), a facility cannot expand (EYet ) 1) more than once in any interval of length CLTe. EYet e 1

t > CLTe

(1)

t-CLTe+1

If a facility e is not under construction at the beginning of the horizon, then no capacity addition is possible for t e CLTe. Hence, the above equations are written for t > CLTe only. In practice, a few large capacity additions would be preferred over several small ones, because a capacity addition involves a lot of construction work and it is convenient to have fewer additions. Therefore, it may be desirable to limit the number of capacity additions and technology upgrades for a facility by using,

∑ EY

et

t

disinvested facility, because a disinvested facility cannot be purchased backsESYet is used to represent this) at time t + 1 (EYet+1 ) 1), but a disinvestment is not possible because the facility does not exist in the first place. A facility e cannot change its status (exists or does not) for t e min[CLTe,TDe] and this is given by equation FX1. All equations with labels starting with “FX” are not written as explicit constraints, but are meant for fixing variables. EEYet ) EEYe,1

t e min[CLTe, TDe]

(FX1)

If a facility starts production due to a capacity addition at time t, but did not exist at time t - 1, then it should be taken as existing at time t. In other words, EYet ) 1 and EEYe(t-1) ) 0 wEEYet ) 1. EEYet g EYet

t > CLTe

EEYet e EEYe(t-1) + EYet' - EXet'' t > min[CLTe, TDe], t' > CLTe,

(3) t'' > TDe (4)

Equation 4 also ensures that a facility cannot be disinvested at t, if it does not exist at t - 1, and it can exist at t, only if it either existed at t - 1 or it adds capacity at t. Similarly, it also ensures that a facility cannot add capacity and disinvest at the same time. Also, a facility existing at time t - 1 must either exist or be sold off at time t (eq 5). EEYe(t-1) e EEYet + EXet'

t' > TDe, t > min[TDe, CLTe]

(5)

Once a facility is disinvested, it cannot be purchased. Equations 6-8 ensure this. ESYet g EXet ESYe(t+1) g ESYet

t > TDe

(6)

H > t > TDe

(7)

A facility cannot exist and be sold off at the same time.

t

∑

11097

+ TUEt' e NEe

t > CLTe,

(2)

As mentioned, technology upgrade is also considered in this limit on the number of expansions. UP is the set of facilities that can undergo a technology upgrade. Now, if a facility e exists (EEYet ) 1) at time t, a capacity addition (EYet ) 1) does not change its status, but a disinvestment is possible at time t + 1 (EXet+1 ) 1). If a facility e does not exist at time t (EEYet ) 0), then a capacity addition (Greenfield addition) is possible (if this facility is not a

t > TDe

(8)

Once a facility undergoes a technology upgrade (TUet ) 1wUFet ) 1), it is taken as an upgraded facility (UFe(t-1) ) 1wUFet ) 1). If a facility is not undergoing an upgrade at the beginning of the horizon, then it cannot be upgraded before its construction lead time for upgrade (TUPe). However, if an upgrade is in progress at the beginning of the horizon, then the variables must be initialized appropriately. UFet ) UFe(t-1) + TUet

t e TUPe,

e ∈ UP (9)

UFet ) 0

E ∈ UP, t' > TUPe

ESYet + EEYet e 1

t > TUPe,

e ∈ UP

(FX2)

2.2. Reformulation of Strategic Binary Variables. For reformulation, we define two variables for facility existence (EEYse, EEYde), one variable for technology upgrade (TUse), and three variables for three capacity expansions (rEYie). All of these are real numbers between 0 and 1. In the reformulation, EEYet, ESYet, and EXet can be replaced by two variables: (a) the time at which the facility begins to exist EEYse, and (b) the duration for which a facility exists, EEYde. Our hybrid algorithm generates strings with variables EEYse, EEYde, TUse, and rEYie and reformulates these variables to obtain the binary variables in the MILP as described below:

11098


Initialization: EEYet ) 0; ESYet ) 0; EYet ) 0; EXet ) 0; TUet ) 0

∀e, t

and warehouse facilities that process the same material either expand or are disinvested at the same time. This is done to eliminate unnecessary evaluations of the objective when a production facility expands and the corresponding warehouses do not or when a production facility is disinvested, but the corresponding warehouses are not. Once the algorithm searches for the near-optimum using the first stage evaluation for a prespecified time, it would use the second evaluation stage, which does not employ any heuristic. We use 60% of the computational time on the first stage and the rest on the second stage. For this two stage approach, it is important to properly identify those facilities that are interlinked and for the SCR problem it is not difficult to perform this analysis for a person working in this domain. Hence user knowledge can be put to use and a solution could be obtained in less time, when compared to the single stage approach. In the reformulation, we use “n” variables to define the number of expansions allowed at a given facility, which also define whether the nth expansion occurs. Hence, we have only n instead of H variables to model capacity expansions at a given facility. It can thus be seen that 3H binary variables (EYet, EXet, and TUet) and 2H continuous variables (EEYet and ESYet) are represented by three (EEYse, EEYde, and TUse) + n(rEYie) variables in the reformulation. Since the number of expansions is usually limited to three or four in a planning horizon, the total number of variables in the reformulation would be only about six or seven for a given facility. This is far less than the 120 (H ) 40) binary and 80 continuous (0-1) variables required in an MILP. The six variables can then be decoded to a feasible set of binary and continuous variables using the decoding procedure described earlier. The binary and continuous variables thus obtained are then passed to the MILP solver, where these are fixed, and the smaller MILP (tactical decisions are still in MILP) is solved to a gap of 1% and the objective is passed to the hybrid algorithm. The 0-1 continuous and binary variables are constrained due to the construction lead time, the assumption that a facility once disinvested cannot be purchased back, and the limits on the numbers of expansions and technology upgrades. It can thus be seen that some of the 0-1 continuous and binary variables that are constrained in the MILP become unconstrained in the reformulation. 3. Evolutionary Algorithm

In the reformulation, the existence of the facilities is remodeled such that no infeasibilities are generated in the MILP and the number of decision variables is reduced. However, the binary variables for capacity expansion (EYet) may create infeasibilities. This is due to (a) the three expansions are generated independently by the hybrid algorithm and they may overlap each other due to the construction lead time, and (b) the investment in a new facility would create an additional EYet ) 1 and this may cause four additions, which is infeasible and has to be repaired. The expansions that may occur prior to the existence of a facility are repaired in the reformulation by the statement that no expansion is allowed if a facility does not exist. Hence, the instance when the facility begins to exist (for a facility that does not exist at the beginning of the horizon) can be considered as the first expansion. If any expansion overlaps with a previous one due to its construction lead time, then this expansion is considered not to occur and the corresponding binary variable is set to zero. The overview of this work is presented in Figure 1. 2.3. Expansion Heuristic. To facilitate the search, the objective is evaluated in two stages. In the first stage, all the facilities that are interlinked, that is, all production facilities

The many algorithms such as genetic algorithm, differential evolution, simulated annealing etc., are convenient, but no single algorithm works well for all problems. The specific choice of one versus another, in our opinion, is mostly a matter of user preference rather than performance. Therefore, we have synthesized our algorithm by combining features of many different existing algorithms. Hence, we prefer to call it as a generic “evolutionary algorithm”, rather than assigning a specific name. Since we combine MILP and evolutionary algorithm (EA), we call it a “hybrid” evolutionary algorithm. Population-based EAs select one or more strings from the population and create new strings by modifying them. These modifications are termed mutations and crossovers. The new strings are compared with the old and a selection procedure selects the best of them for the next generation. In general, all EAs need to be fine-tuned for a given problem by adjusting its parameters. We select our parameters either from a uniform distribution or as the average of the bounds, as we believe that arbitrary and random trials to decide the best parameters are not practical for new and large problems.


3.1. Initialization. We generate two strings, one at each of the bounds (string of zeros and string of ones) and calculate the objective for each. Next, we create a temporary population by augmenting the child strings created from a generation to the initial population. In this temporary population, the top part consists of the latter, and the bottom part is an empty/random population that is updated, as mutations and crossovers are performed. For reproduction, parents are selected only from the top part. 3.2. Reproduction. For describing the mutations and crossovers for reproduction, we employ the following nomenclature. A string is a set of variables. r1(i) is a random string of zeros and ones and r2(i) is a random string of real numbers between zero and one, where “i” indicates the ith variable in the string and ones(i) is a string of ones. The length of r1, r2, and ones(i) is equal to the length of the parent string. 3.2.1. Mutations. A multinon-uniform mutation is performed on a single parent (p) to obtain a single child (c). If p(i) is not at the upper bound, then c(i) ) p(i) + (1 - p(i)) r1(i) (1 - m)3 If p(i) is not at the lower bound, then c(i) ) p(i) - p(i) r1(i) (1 - m)3 where, m ) time when current generation started/total run time. If p(i) is not at the bounds, c(i) can take either of the above two values with equal probability. A nonuniform mutation is similar to the multinon-uniform mutation, but is performed on only one randomly selected variable in a string. A uniform mutation is performed on a single parent to obtain a single child. This is performed only on a single variable in the parent string. A random variable is selected first. If it is above 0.5, then a child string is obtained by replacing it by a random number below 0.5 and vice versa. A boundary mutation is similar to a uniform mutation except that the randomly selected variable is pushed to the bounds. If it is already at one bound, then it is flipped to the other bound. Four strings (p1, p2, p3, and p4) are selected from the original population and one child is formed in differential eVolution as follows: c(i) ) p1(i) r1(i) + [p2(i) ( r2(i){p3(i) - p4(i)}]{1 - r1(i)}

In the above child, addition is initially selected as a thumb rule in place of “(”. This does not bias the mutation toward any particular string/direction, because p3 and p4 are randomly selected, and their difference is unrestricted. If a variable goes beyond a bound because of the addition operation, then that particular value is discarded and subtraction is done. Thus, all variables are kept within bounds. While it is common practice to reset the values beyond the bounds to be at the bounds, and this may produce the same string in some cases, we maintain diversity by replacing addition by subtraction. 3.2.2. Crossovers. Two strings (p1 and p2) are selected from the original population and two child strings are formed in a random crossoVer as follows: c1(i) ) p1(i) r1(i) + p2(i)[1 - r1(i)] c2(i) ) p1(i)[1 - r1(i)] + p2(i) r1(i) An arithmetic crossoVer is performed on two parents to obtain two child strings. The two child strings are produced as follows:

11099

c1(i) ) p1(i) r2(i) + p2(i)[1 - r2(i)] c2(i) ) p1(i)[1 - r2(i)] + p2(i) r2(i) An arithmetic crossover can be done on all variables or only on some randomly selected variables, while retaining the others from the parent string. Note that the difference between the random and arithmetic crossovers is the random number, which is binary in for the former and 0-1 continuous for the latter. A heuristic crossoVer is an arithmetic crossover where one parent is selected randomly to produce two child strings. The first (second) child is created by the parent undergoing an arithmetic crossover with a dummy parent with all variables at their upper (lower) bounds as follows. c1(i) ) p1(i) r2(i) + [1 - r2(i)] c2(i) ) p1(i)[1 - r2(i)] Clearly, several types of crossovers and mutations are possible, and one has to decide the sequence in which and how often they should be used before a new generation is completed. Since there is no obvious reason why one crossover or mutation should be selected before others, we decide a priori the sequence and frequency of each as follows: 6 random crossovers, 3 differential evolutions, 12 multi-nonuniform mutations, 3 uniform mutations, 3 arithmetic crossovers, 3 nonuniform mutations, and 6 boundary mutations. While the crossovers produce two children each, the mutations produce one child each. Overall, the operations produce 51 child chromosomes in each generation. 3.3. Selection. Our optimization algorithm can be used either to identify one good solution or for identifying multiple different solutions. The difference between the two is performance. When the algorithm is used to identify a single solution, it generally gives a relatively better objective value. We present the two cases separately. 3.3.1. Procedure for Single Solution. When a child is produced, its objective is compared with that of its parent. When there are two parents and two children, the first child is compared with the first parent and the second with the second parent. Ideally, the two child chromosomes must be compared with both the parents; however, we eliminate this step, as the selection procedure described next obviates its need. If the parent is in the best 60% of the top part (initial population) of the temporary population, and the child’s objective is better, then the child replaces the parent and the parent is moved to the bottom part (child population) of the temporary population. If not, the parent retains its position, and the child is placed in the bottom part of the temporary population. If the parent is in the worst 40% of the initial population (the top part of the temporary population), then the child always replaces its parent, and the parent is placed in the bottom part of the temporary population. Hence, a better child may participate in the reproduction for the next mutation/crossover in the current generation. In this way, the bottom part of the temporary population is filled, when all mutations and crossovers end. The bottom part of the temporary population has the child strings produced by all mutations and crossovers. When all the mutations and crossovers end, the temporary population consists of all members in the initial population and all the child strings. From this population, the best strings can be selected to form the initial population for the next generation. Our updating procedure is an advanced version of the updating procedure in modified differential evolution (MDE).33 In MDE, either the parent or the child string is discarded, if it is found to be inferior. This may sometimes lead to premature

11100


convergence, because the discarded string may be better than the other strings in the population. Angira and Babu33 did not consider this. One modification to their procedure is to consider both the initial population of the current generation and the current population for updating the population at the end of the generation, so that some of the lost parents could be considered. However, if a child that replaced a parent participates in mutations/crossovers and gets replaced by its child, then it gets lost permanently. A second modification is to replace the child with the worst string in the population every time. However, this also has a drawback of finding the worst string each time a mutation/crossover is performed. This is computationally expensive for problems that require a large population. Our method is an alternative to the second method described above and is computationally less expensive, because the selection of best strings is done only once every generation. The best available strings on hand may not always participate in reproduction, but one has to trade-off between computation time and choosing the best strings for reproduction. Also, by storing the parent/child in the bottom part of the temporary population, no strings are lost. For the SCR problem, our advanced updating procedure gave better results than that of Angira and Babu.32 In addition to the above changes, we check if a string is evaluated previously, first by checking the temporary population and later by checking the set of strings that were evaluated in the past. Thus, we reduce redundant evaluations and improve computational time. One notable difference between our algorithm and other algorithms is the population size. It is generally believed that the population size should be about 5-10 times the number of variables in a string. Considering the large number of variables in the SCR problem, such a large population set would only increase the computation time, because the algorithm would be spending more time on strings that may not be very useful. In our hybrid algorithm, we have used a population size of only 50 for a problem with over 100 variables without compromising solution quality. Only for multimodal optimization do we use a larger population size. 3.3.2. Procedure for Multimodal Optimization. One advantage of using a population-based optimization algorithm is the ability to get multiple solutions. While distinct multiple solutions (called multimodal optimization) would be preferred, the algorithm may get stuck in one region of the search space, if proper diversifying criteria are not used. To identify multiple solutions, we use a niche count34 (a diversity preserving factor). First, a secondary string is defined as the string containing the binary variables modeling the existence or upgrade of facilities. For a network of e facilities and horizon length H, the length of this secondary string would be 2eH, eH for facility existence and eH for technology upgrade. We define this, because it would be sufficient to compare two supply chains by comparing variables at the same locations in the strings and checking for similarity (both facilities exist or do not exist) or dissimilarity (one facility exists and other does not). Now, we define dij as the number of variables that are dissimilar between two secondary strings i and j. Clearly, dij cannot exceed the length of the secondary string. Next, σshare is defined as the extent to which two strings are similar. This is the fraction of the length of the secondary string used to measure similarity/dissimilarity threshold. Then, a sharing function Shij of strings i and j in the population is defined as

(

Shij ) 1 -

)

dij R if dij e σshare and Shij ) 0 otherwise σshare

Table 1. Material Balance for Various Production Facilities in Example 1 facility (nation)

mass balance for production facility

52 (1) 53 (1) 54 (1) 55 (2) before upgrade 55 (2) after upgrade 56 (3) 57 (4) 58 (3) before upgrade 58 (3) after upgrade 59 (3) 60 (2) 61 (4)

m6 + 0.5m7 ) 0.7m20 + 0.8m16 m8 + 0.7m9 ) 0.9m17 + 0.8m21 m16 + 0.6m17 ) 0.8m19 + 0.8m22 0.9m10 + 0.7m19 ) 0.8m23 + 0.8m24 0.9m10 + 0.7m19 ) 0.8m23 + 0.8m31 0.7m11 + 0.9m12 ) 0.9m18 + 0.7m25 m13 + 0.4m18 ) 0.6m26 + 0.8m27 0.8m14 + 0.9m15 ) m28 + 0.7m29 0.8m14 + 0.9m15 ) m28 + 0.7m30 m6 + 0.5m7 ) 0.7m20 + 0.8m16 m16 + 0.6m17 ) 0.8m19 + 0.8m22 0.7m11 + 0.9m12 ) 0.9m18 + 0.7m25

Shij for a given string with itself is one and that with a different string is zero. Now, niche count nci is defined as nci ) ∑jShij. The niche count thus calculated is used in general multimodal optimization to guide the selection of chromosomes for reproduction. It can be seen from the definition that a high niche count for a given string in a population means that several strings in a population are similar to the string under consideration and a niche count of one means that the string under consideration is totally dissimilar to the rest. Taking this into account, a secondary fitness function (say f 2) is defined as f 1/nci, where f 1 is the actual fitness/objective. A high value of f 2 means that a string is under-represented in the population, and vice versa. f 2 is used as an objective for applying the selection operator so that an under-represented string can participate more in reproduction and thus maintain population diversity.34 In this work, we use σshare ) 0.25 L (L ) length of the string) and R ) 1. Unlike in the literature,34 we do not use f 2 for applying the selection operator, but use random selection as defined previously in section 2.3. Once we have the final population at the end of a generation, that is, the temporary population that is filled, we define f 2 as above and f 3 as nci-1. The temporary population is sorted three times on the basis of the three objectives (f 1, f 2, and f 3) and we select the equivalent of one-third of the initial population from each of the sorted lists to obtain the full starting population for the next generation. By doing so, the strings that are good in all the lists would have more copies, and elitism and diversity are preserved. For sorted lists based on f 2 and f 3, a string is sent to the starting population only when its corresponding f 1 is greater than 75% of the best f 1. If the f 2 or f 3 sorted lists do not have sufficient members which are above 75% to fill one-third of the starting population, then additional members from f 1 sorted list (those that are not chosen previously) are chosen to fill the starting population. 3.4. Analysis of Final Population. 3.4.1. Single Solution. When the selection procedure for a single solution is used, we obtain the final population where the objective of each solution is usually different. However, a closer analysis reveals that the supply chain configurations of solutions are largely the same and the differences in objective are due to some minor changes in operational decisions such as the amounts of materials stored or transported. 3.4.2. Multiple Solutions. Using the selection procedure for multiple solutions, we obtain a final population that is diverse. Once multiple optimal supply chain networks are generated by using the two additional objectives, multiple networks must be identified automatically from the list. This is done by using simple clustering based on the value of nci of the final population and selecting the best network from each of the clusters. We use K-means iteratively and select the solution that has the


11101

Table 2. Location of Inventory Facilities and the Material They Hold in Example 1 facility (nation)

inventory facility (material it holds)

32-37, 62-67, 82-83, 94-99 (1) 38-39, 48-49, 68-69, 78-79, 84-86, 92, 100-108, 124-126 (2) 40-41, 44-47, 70-71, 74-77, 87-89, 109-117 (3) 42-43, 50-51, 72-73, 80-81, 90-93, 118, 123, 127-129 (4)

32 (6), 33 (7), 34 (8), 35 (9), 36 (16), 37 (17), 62 (20), 63 (16), 64 (17), 65 (21), 66 (19), 67 (22), 82-83 (20-31), 94-99 (20-31) 38 (10), 39 (19), 48 (16), 49 (17), 68 (23), 69 (24, 31), 78 (19), 79 (22), 84-86 (20-31), 92 (20-31), 100-108 (20-31), 124-126 (20-31) 40 (11), 41 (12), 44 (14), 45 (15), 46 (6), 47 (7), 70 (18), 71 (25), 74 (28), 75 (29, 30), 76 (20), 77 (16), 87-89 (20-31), 109-117 (20-31) 42 (13), 43 (18), 50 (11), 51 (12), 72 (26), 73 (27), 80 (18), 81 (25), 90-93 (20-31), 118-123 (20-31), 123 (20-31), 127-129 (20-31)

Table 3. Computational Statistics for Full MILP Model in GAMS and the Hybrid Evolutionary Algorithm for Example 1 MILP

hybrid

binary variables 0-1 continuous variables constraints variables nonzeros CPU time (days)

15 888 10 720 189 834 18 134 1 765 314 5

NPV ($ bn) best possible NPV ($ bn)

-2.43 10.17

8800 in sMILP 7088 from GA chromosome length ) 417 + 1 7920 continuous 0-1 variables are fixed from GA 189834 18134 103302 6 (400s-600s for each sMILP, including model loading time, to reach 1% gap) 1.45 10.17

Table 4. Capacity Profiles of the Production Facilities for Example 1 facility (nation)

capacity (kton/qtr) (using MILP)

capacity (kton/qtr) (using hybrid algorithm)

52 (1) 53 (1) 54 (1) 55 (2) 56 (3) 57 (4) 58 (3) 59 (3) 60 (2) 61 (4)

6000 for periods 1-7, 0 for periods 8-40 6000 for periods 1 to 4, 0 for periods 5 to 40 6000 for periods 1 to 8, 0 for periods 9 to 40 6000 for periods 1 to 40 6000 for periods to 8, 0 for periods 9 to 40 6000 for periods 1 to 10, 0 for periods 11 to 40 6000 for periods 1 to 40 does not exist in the planning horizon does not exist in the planning horizon does not exist in the planning horizon

6000 for periods 1-40 6000 for periods 1 to 10, 0 for periods 11 to 40 6000 for periods 1 to 40 6000 for periods 1 to 40 (no technology upgrade) 6000 for periods 1 to 40 6000 for periods 1 to 40 6000 for periods 1 to 28, 0 for periods 29 to 40 (no technology upgrade) does not exist in the planning horizon does not exist in the planning horizon does not exist in the planning horizon

Table 5. Computational Statistics for Full MILP Model in GAMS and the Hybrid Evolutionary Algorithm for Example 2

binary variables 0-1 continuous variables constraints variables nonzeros CPU time (hrs) NPV ($ bn) best possible NPV ($ bn)

MILP

hybrid

2220 2160 26 390 18 134 103 302 24 15.08 16.59

400 in sMILP 1,820 from GA chromosome length ) 108 + 1 2040 continuous 0-1 variables are fixed from GA 26 390 18 134 103 302 120 (30-50s for each reduced MILP) 14.6 16.59

minimum variance. As a general rule in the iterations, the maximum number of clusters is set arbitrarily between 10 and 15. Alternatively, we plot nci versus i and determine (visually) the number of clusters as the number of bands. The multiple supply-chain networks thus identified are used for further analysis to study performance under different scenarios or uncertainty.

3.6 GHz processor. The hybrid algorithm was first tested on another example for which a single global optimal solution exists. While the data is not shown here, the example is of a

4. Examples The SCR problem is heavily constrained and hence it is unlikely that the proposed hybrid algorithm can give a better solution than that possible with GAMS for small-size problems. However, our hybrid algorithm is suitable for multimodal optimization and has the potential to fare better when the problem size is too large for an MILP solver. To demonstrate these two features of our algorithm, we solve a large SCR problem in example 1 and a multimodal optimization problem in example 2. The algorithm is implemented in Matlab R2007b and the MILP model is solved using CPlex 9.0 in GAMS 21.7 on a Windows XP-based HP workstation with an Intel Xeon (Dual)

Figure 2. Capacity profiles for example 2 from the full MILP. The technology upgrade of facility 1 occurs at period 22, and the profit is $15.08 bn.

11102


Figure 3. Part 1 of 3

size similar to that of example 2 in this paper. GAMS/CPlex takes about 1 day of CPU time to get a solution. Our algorithm attains >97% of the same solution in about 5 days of CPU time. This shows that the algorithm is able to find a reasonably good solution, even if it takes a longer time. To illustrate the performance on a much larger example, we use example 1. All examples are solved with the objective of maximizing the shareholder value.

4.1. Example 1sLarge SCR Problem. The supply-chain network (SCN) consists of 5 raw material suppliers supplying 10 raw materials and 4 intermediate materials, 10 production facilities (7 existing and 3 potential facilities) producing the 4 intermediates and 12 products. There are 75 existing facilities and 23 potential ones. Two inventory facilities exist at the input to each product facility and two at the output. The products and intermediates are sent to 12 distribution centers (10 existing


11103

Figure 3. Part 2 of 3

and 2 potential facilities) and these serve 36 customer centers (30 existing and 6 potential facilities). Our model includes inflation, depreciation, regulatory factors, material supply contracts, loans and bonds for raising capital for investments, maintenance shutdowns, and stockouts among others. To solve this problem, we consider a planning horizon of 10 years divided into 40 periods of equal length. The details of the supply chain showing the material flow and storage are given in Tables 1 and 2. For the sake of focus, we give only the basic network

structure to compare the two solutions and do not give the detailed cost parameters. The MILP model for this large SCR problem has 15888 binary variables. Of these, 7088 binary variables are related to EYet, EXet, TUet and the remaining are related to the contracts for the supply and production of materials. Since, the latter (8800) are relatively less constrained, they are easy to solve and we do not need to reformulate them even for a large SCR problem. The computational complexity and difficulty of solving

11104


Figure 3. Part 3 of 3 (a) Capacity profiles for network 1 in example 2 obtained from the hybrid algorithm in the multimodal optimization mode. Technology upgrades do not occur and the profit is $14.52 bn. (b) Capacity profiles for network 2 in example 2 obtained from the hybrid algorithm in the multimodal optimization mode. Technology upgrades do not occur and the profit is $14.30 bn. (c) Capacity profiles for network 3 in example 2 obtained from the hybrid algorithm in the multimodal optimization mode. Technology upgrades of facility 1 occurs at period 8 and the profit is $14.09 bn. (d) Capacity profiles for network 4 in example 2 obtained from the hybrid algorithm in the multimodal optimization mode. Technology upgrades of facility 1 occurs at period 33 and the profit is $14.46 bn. (e) Capacity profiles for network 5 in example 2 obtained from the hybrid algorithm in the multimodal optimization mode. Technology upgrades of facility 1 occurs at period 5 and the profit is $13.62 bn. (f) Capacity profiles for network 6 in example 2 obtained from the hybrid algorithm in the multimodal optimization mode. Technology upgrades of facility 1 occurs at period 8 and the profit is $14.16 bn. (g) Capacity profiles for network 7 in example 2 obtained from the hybrid algorithm in the multimodal optimization mode. Technology upgrades of facility 1 occurs at period 20 and the profit is $14.08 bn. (h) Capacity profiles for network 8 in example 2 obtained from the hybrid algorithm in the multimodal optimization mode. Technology upgrades of facility 1 occurs at period 20 and the profit is $14.12 bn. Table 6. Profit, Standard Deviation, and Skewness of Eight Different Networks Identified by the Hybrid Algorithm network

profit ($) 1.0 × 1010

standard deviation ($) 1.0 × 108

skewness

1 2 3 4 5 6 7 8

14.52 14.30 14.09 14.46 13.62 14.16 14.08 14.12

4.87 5.00 4.67 4.58 4.52 4.68 4.64 4.61

-0.4676 -0.4060 -0.3921 -0.4719 -0.3800 -0.3971 -0.3831 -0.4163

the MILP are due to the 7088 binary variables, and hence only these are reformulated. According to the reformulation, these would require 75 + 2 × 23 ) 121 variables for their existence (only one variable defining the duration for which a facility exists is required for facilities existing at the start of the horizon), 3 × 98 variables for expansions, and 2 variables for technology upgrade, as only 2 facilities are eligible for technology upgrade in this example. Hence, the reformulation has 417 variables and

the length of a chromosome string is 418 with the objective value. In addition to the binary variables, 7920 out of the 10720 continuous (0-1) variables are also obtained from the reformulation and fixed in the MILP. We use a population size of 50 for this large problem. The computational statistics are given in Table 3. It can be seen that GAMS/CPlex gives a solution with a profit of -$2.43 bn after five days of CPU time, but the hybrid algorithm gives a far better profit of $1.45 bn. However, the profit from the hybrid algorithm ($1.45 bn) is far less than the best profit ($10.17 bn, corresponding to upper bound of the relaxed MILP, when the optimization is stopped) possible from GAMS/CPlex. Further fine-tuning of the reproduction operators and selection procedures is necessary to improve the quality of the solution. Table 4 shows the capacity profiles of the production facilities. The solution from the hybrid algorithm has fewer disinvestments and no capacity additions or technology upgrades, which possibly increase the profit. In this work, we used the boundary conditions of the chromosomes for initializing the hybrid


11105

Figure 4. Distribution of profit under different scenarios for the different supply chain networks that are identified and presented in Figure 3.

algorithm. In reality, the planners may have a tentative plan that may be used for initializing the hybrid algorithm for further improvement and/or finding alternative plans. 4.2. Example 2 - Multimodal Optimization. A multiechelon SCN consisting of two material suppliers who supply five raw materials and one intermediate, four production facilities (three existing and one potential) producing one intermediate and five products, and five distribution centers (three existing and two potential) is considered. In addition, each production facility has two input and two output inventory facilities. The production facility that can potentially be disinvested and the one that can be newly invested in are in two different nations. Hence, the problem can be considered to be that of relocation from one nation to another. All features that are present in example 1 are present in example 2. However, the parameters are chosen such that multiple optimal supply chain networks are possible. The computational statistics are presented in Table 5. The MILP model implemented in GAMS had 2220 binary variables. In the hybrid algorithm, we have 1820 strategic and 400 tactical binary variables. The reformulation reduces the number of variables to only 108 variables. Furthermore, 2040 of the 2160 continuous (0-1) variables are also obtained from the hybrid algorithm and fixed in GAMS/CPlex. Thus, the hybrid algorithm requires only about 30-50 s of CPU time to reach a gap of 1% for the reduced MILP. In multimodal optimization, each potential optimal solution requires a subpopulation, from which it can evolve. Thus, the population for a multimodal optimization problem can be considered to have “n” subpopulations, where n is the number of optimal solutions. As we do not know the number of solutions a priori, we chose a larger population of 300 for multimodal optimization as compared to only 50 for the single-solution optimization. As mentioned earlier, a population size of only 50 (even for a large problem such as example 1) as opposed to the usual size of 5-10 times the string length is used for single-solution optimization. Thus, the performance of our algorithm does not appear to be severely affected by the population size. In view of this, we concluded that a population of 300 strings would suffice for multimodal optimization.

In addition to solving this example for multimodal optimization, we also solved it for single-network optimization. To this end, we used a population size of 50 in our hybrid algorithm, and removed the component of the selection operation (section 3.3.2) that helps us identify multiple solutions. This allows us to benchmark the hybrid algorithm with GAMS/CPlex, as both can solve the single-network optimization. The capacity profiles of the production facilities for the solution obtained from GAMS/CPlex and those obtained from the hybrid algorithm are given in Figures 2 and 3, respectively. In Table 6, we present the profit, standard deviation, and skewness of eight networks that the hybrid algorithm identified. We present only the capacity profiles in Figure 3. It is interesting to note that the single network obtained from GAMS/CPlex is better and different from the eight from the hybrid algorithm. The profits of the networks from the hybrid algorithm are only about 90-95% of that for the network from GAMS/CPlex. Figure 2 shows the single-solution network from the hybrid algorithm. Its profit is indeed higher than those of the eight from the multimodal optimization. It reaches >97% of the profit from GAMS/CPlex. To compare the eight networks from the standpoint of robustness, we now perform a scenario analysis with varying production costs, taxes, and demands. A total of 5000 scenarios are generated by assuming normally distributed parameters with means as the nominal values used in example 2. The scenarios are stored separately and are used to analyze the effect of uncertainty on each network. The distribution of profit across various scenarios is shown in Figure 4. Networks 1 and 4 have the highest mean profits compared to the others. One feature common to these two networks is that the capacity of facility 4 occurs only toward the end of the horizon and it is much lower than those in the other networks. This suggests that the investment in facility 4 should be delayed and its capacity should be kept to a minimum in order to maximize the profit. Next, we analyze the mean and skewness of each of the distributions obtained on performing the scenario analysis. One clear trend from Figure 3 and Table 6 is that a higher skewness is related to either higher investment in production capacity “4”

11106


(PF4) or existence of PF4 for long duration. That is, the investment in PF4 reduces the chances of lower profit (compared to the mean profit of the same network) due to uncertainty. However, considering the profit that a given network can make, network 1 where PF4 has only a small investment at the end of the planning horizon is a good option. Ideally, one would like to have one network that is superior in terms of all three criteria of profit, skewness, and standard deviation. However, no such network exits in this example and the decision maker must make a choice based on his preference of higher profit (skewness) versus lower risk (standard deviation). 5. Conclusions Redesign of large supply chains involves many decisions and a monolithic approach based on MILP fails. We have developed a novel hybrid algorithm that (1) combines the strengths of MILP and evolutionary algorithms by judiciously and intelligently reformulating constraints involving the binary and 0-1 continuous variables and (2) gives multiple supply-chain networks with different configurations. In a large example, our strategy reduced the number of binary variables from 15888 to 8800 binary variables in the reduced MILP and 417 decisions in the evolutionary algorithm. Our hybrid evolutionary algorithm outperforms the monolithic MILP approach for large-scale supply chain redesign. In another smaller example, our algorithm successfully identified eight different optimal networks with profits exceeding 90% of the best obtained from the MILP approach. For smaller supply chains and when only one network is desired, the hybrid algorithm attains >97% of the best profit from the MILP approach. When used in conjunction with a scenario analysis using uncertain data, our algorithm enables the selection of a robust network from several nearly optimal supply-chain networks. This is a significant contribution, as most business data, parameters, and projections are highly uncertain for most strategic planning scenarios in real life.

Binary Variables EXet ) 1, if a facility is sold and the capacity is no longer available at t EYet ) 1, if an added capacity is available for production at t EZet ) 1, if a facility is in operation at t TUet ) 1, if a facility undergoes a technology upgrade at t Continuous 0-1 Variables EEYet ) 1, if an facility exists at t ESYet ) 1, if an facility that existed before t is sold off at t or before UFet ) 1, if a facility has undergone a technology upgrade PositiVe Variables dij ) the total number of variables that are dissimilar between secondary strings i and j. EEYse ) time when a facility begins to exists (as a fraction of planning horizon H) EEYde ) the duration a facility exists (as a fraction of the period which it can exist after beginning to exist) rEYie ) time (as a fraction of the period where expansions can occur) at which the ith expansion occurs nci ) niche count, showing how different a given string is with others Shij ) sharing function, defining the similarity of two strings i and j. TUse ) time (as a fraction of the period where technology upgrade can occur) at which a technology upgrade occurs. σshare ) the fraction of the length of the secondary string used to measure similarity/ dissimilarity threshold Strings r1i ) a random string of zeros and ones r2i ) a random string of numbers between zeros and ones ci, c1i, c2i ) a child chromosome, where i is the ith element pi, p1i, p2i ) a parent chromosome, where i is the ith element Literature Cited

Acknowledgment We would like to acknowledge the financial support for this work under the project ICES/06-432001. We would like to thank Prof. Carlos Artemio Ceollo Coello, Instituto Polite´cnico Nacional, Mexico for his valuable suggestions. Nomenclature Note: Superscript L (U) denotes the lower (upper) limit or the lowest (highest) value in a parameter vector. Indices e ) an entity (raw material, supplier, inventory or production facility) t ) period Sets UP ) set of production facilities that may undergo a technology upgrade Parameters CLTe ) construction lead time for a facility NEe ) maximum number of times that a facility may expand PTpe ) interval in which routine maintenance must be done at least once for a facility STle ) maximum number of periods for which a facility may operate during PTpe (STle < PTpe ) TUPe ) initial interval during which no technology upgrade is allowed for a facility

(1) Naraharisetti, P. K.; Karimi, I. A.; Srinivasan, R. Supply Chain Redesign through Optimal Asset Management and Capital Budgeting. Comput. Chem. Eng. 2008, 32, 3153–3169. (2) Karimi, I. A.; Srinivasan, R.; Han, P. L. Unlocking Supply Chain Improvements through Effective Logistics. Chem. Eng. Progress 2002, 32– 38. (3) Bansal, M.; Karimi, I. A.; Srinivasan, R. Selection of Third-Party Service Contracts for Chemical Logistics. Ind. Eng. Chem. Res. 2008, 47, 8301–8316. (4) Oh, H. C.; Karimi, I. A. Global Multiproduct Production-Distribution Planning with Duty Drawbacks. AIChE J. 2006, 52, 595–610, 2. (5) Oh, H.-C.; Karimi, I. A. Regulatory Factors and Capacity-Expansion Planning in Global Chemical Supply Chains. Ind. Eng. Chem. Res. 2004, 43, 3364–3380. (6) Adhitya, A.; Srinivasan, R.; Karimi, I. A. Supply Chain Risk Management through Hazop and Dynamic Simulation. Presented at the 18th European Symposium on Computer Aided Process EngineeringsESCAPE 18, Lyon, France, June 1-4, 2008. (7) McDonald, C. M.; Reklaitis, G. V. Design and Operation of Supply Chains in Support of Business and Financial Strategies. Presented at the Foundations of Computer-Aided Process Design, FOCAPD, Princeton, NJ, July 11-16, 2004. (8) Bjorkqvist, J.; Roslof, J. Strategic Investment Planning in the Pulp and Paper Industry Using Mixed Integer Linear Programming. Presented at the AIChE Annual Meeting, Cincinnati, OH, Oct 30-Nov 04, 2005. (9) Bok, J.-K.; Lee, H.; Park, S. Robust Investment Model for LongRange Capacity Expansion of Chemical Processing Networks under Uncertain Demand Forecast Scenarios. Comput. Chem. Eng. 1998, 22, 1037– 1049. (10) Gatica, G.; Papageorgiou, L. G.; Shah, N. Capacity Planning under Uncertainty for the Pharmaceutical Industry. Trans. Inst. Chem. Eng. 2003, 81, 665–678.

Ind. Eng. Chem. Res., Vol. 48, No. 24, 2009 (11) Levis, A. A.; Papageorgiou, L. G. A Hierarchical Solution Approach for Multisite Capacity Planning under Uncertainty in the Pharmaceutical Industry. Comput. Chem. Eng. 2004, 28, 707–725. (12) Papageorgiou, L. G.; Rotstein, G. E.; Shah, N. Strategic Supply Chain Optimization for the Pharmaceutical Industries. Ind. Eng. Chem. Res. 2001, 40, 275–286. (13) Guillen, G.; Mele, F. D.; Bagajewicz, M. J.; Espuna, A.; Puigjaner, L. Multiobjective Supply Chain Design under Uncertainty. Chem. Eng. Sci. 2005, 60, 1535–1553. (14) Guillen, G.; Mele, F. D.; Espuna, A.; Puigjaner, L. Addressing the Design of Chemical Supply Chains under Demand Uncertainty. Ind. Eng. Chem. Res. 2006, 45, 7566–7581. (15) Arntzen, B. C.; Brown, G. G.; Harrison, T. P.; Trafton, L. L. Global Supply Chain Management at Digital Equipment Corporation. Interfaces 1995, 25, 69–93, 1. (16) McDonald, C. M.; Karimi, I. A. Planning and Scheduling of Parallel Semicontinuous Processes. 1. Production Planning. Ind. Eng. Chem. Res. 1997, 36, 2691–2700. (17) Wilkinson, S. J.; Cortier, A.; Shah, N.; Pantelides, C. C. Integrated Production and Distribution Scheduling on a Europe-wide Basis. Comput. Chem. Eng. 1996, 20, S1275–S1280. (18) Tsiakis, P.; Shah, N.; Pantelides, C. C. Design of Mulitechelon Supply Chain Networks under Demand Uncertainty. Ind. Eng. Chem. Res. 2001, 40, 3585–3604. (19) Gjerdrum, J.; Shah, N.; Papageorgiou, L. G. Transfer Prices for Multienterprise Supply Chain Optimization. Ind. Eng. Chem. Res. 2001, 40, 1650–1660. (20) Goetschalckx, M.; Vidal, C. J.; Dogan, K. Modeling and Design of Global Logistics Systems: A Review of Integrated Strategic and Tactical Models and Design Algorithms. Eur. J. Oper. Res. 2002, 143, 1–18. (21) Perea-Lopez, E.; Ydstie, B. E.; Grossmann, I. E. A Model Predictive Control Strategy for Supply Chain Optimization. Comput. Chem. Eng. 2003, 27, 1201–1218. (22) Ryu, J.-H.; Pistikopoulos, E. N. Design and Operation of an Enterprise-wide Process Network Using Operation Policies. 1. Design. Ind. Eng. Chem. Res. 2005, 44, 2174–2182.

11107

(23) Holland, J. H. Adaptation in Natural and Artificial Systems; University of Michigan Press: Ann Arbor, MI, 1975 (re-issued by MIT Press 1992). (24) Storn, R.; Price, K. Differential Evolution - A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. J. Global Optimiz. 1997, 11, 341–359. (25) Glover, F. Tabu SearchsPart I. ORSA J. Comput. 1989, 1, 190– 206. (26) Glover, F. Tabu SearchsPart II. ORSA J. Comput. 1990, 2, 4–32. (27) Glover, F. A Template for Scatter Search and Path Relinking. Artificial EVolution, Lecture Notes in Computer Science;, Hao, InL J.-K., Lutton, E., Ronald, E., Schoenauer, M., Snyers, D., Eds.; Springer: New York, 1998; Vol. 1363, pp 3-51. (28) Dorigo, M.; Gambardella, L. M. Ant Colony System: A Cooperative Learning Approach to the Traveling Salesman Problem. IEEE Trans. EVol. Comput. 1997, 1 (1), 53–66. (29) Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by Simulated Annealing. Science 1983, 220, 671–680. (30) Cerny, V. A Thermodynamical Approach to the Travelling Salesman Problem: An Efficient Simulation Algorithm. J. Optimiz. Theory Appl. 1985, 45, 41–51. (31) Ferris, M. C. Interfacing Optimization and Visualization Software; Computer Sciences Department, University of Wisconsin: Madison, WI, 1998. (32) Michalewicz, Z.; Janikow, C. Z. Handling Constraints in Genetic Algorithms. Proceedings of the Fourth International Conference on Genetic Algorithms, San Diego, CA, July, 1991. (33) Angira, R.; Babu, B. V. Optimization of Process Synthesis and Design Problems: A Modified Differential Evolution Approach. Chem. Eng. Sci. 2006, 61, 4707–4721. (34) Goldberg, D. E.; Richardson, J. Genetic Algorithms with Sharing for Multimodal Function Optimization. Proc. First Int. Conf. Genet. Algorithms Their Appl. 1987, 41–49.

ReceiVed for reView February 17, 2009 ReVised manuscript receiVed July 17, 2009 Accepted September 15, 2009 IE9002574

Supply Chain RedesignâMultimodal Optimization Using a Hybrid

Recommend Documents

Supply Chain RedesignâMultimodal Optimization Using a Hybrid