Addressing the Design of Chemical Supply Chains under Demand

7566

Ind. Eng. Chem. Res. 2006, 45, 7566-7581

Addressing the Design of Chemical Supply Chains under Demand Uncertainty Gonzalo Guille´ n, Fernando D. Mele, Antonio Espun˜ a, and Luis Puigjaner* Department of Chemical Engineering, UniVersitat Polite` cnica, de Catalunya, ETSEIB, AVda. Diagonal, 647 E-08028, Barcelona, Spain

This article addresses the design of chemical supply chains (SCs) under demand uncertainty. Given a design horizon consisting of several time periods in which expected but uncertain product demands materialize, the objective is to select the design that maximizes the expected profit. To tackle this problem, a multistage stochastic formulation is derived wherein the design decisions are made irrespective of the specific realization of the uncertain parameters. The resulting model is, in general, computationally intensive or even intractable, as the number of variables and equations is proportional to the number of possible nodes on the scenario tree, which explodes exponentially with the number of stages (periods). In view of this, a decomposition technique is introduced aiming at the objective of overcoming the numerical difficulties associated with the underlying large-scale stochastic mixed-integer linear problem (MILP). The resulting strategy combines genetic algorithms (GAs) and mathematical programming tools. Computational results indicate that the proposed approximation method provides solutions which are within a few percent of the optimal ones and also reduces the computation time required by the rigorous mathematical model. 1. Introduction The concept of the supply chain (SC), which first appeared in the early 1990s, has recently been the focus of much interest, as the possibility of providing an integrated management of a SC can reduce the propagation of unexpected/undesirable events throughout the network and can markedly improve the profitability of all the parties involved. Supply chain management (SCM) aims at the objective of integrating plants with their suppliers and customers so that they can be managed as a single entity and coordinating all input/output flows so that products are produced and distributed in the right quantities, to the right locations, and at the right time.1 The SCM problem may be considered at different levels depending on the strategic, tactical, and operational variables involved in decision making.2 Therefore, a large spectrum of a firm’s strategic, tactical, and operational activities are encompassed by SCM: (i) The strategic level concerns those decisions that will have a long-lasting effect on the firm. It is focused on SC design and entails determining the optimal configuration for an entire SC network, including the design of the embedded plants. (ii) The tactical level encompasses long/medium term management decisions, which are typically updated at a rate ranging between once every quarter and once every year. These include overall purchasing and production decisions, inventory policies, and transport strategies. (iii) The operational level refers to day-to-day decisions such as scheduling, lead-time quotations, routing, and lorry loading. The complexity of an SC is increased by the high degree of uncertainty brought about by external factors, such as continuously changing market conditions, customer expectations, and internal parameters, such as product yields, qualities, and processing times. The design and operation of a SC can both be posed as largescale dynamic decision problems.3 These problems are receiving increasing attention in both the operations research (OR) and the process systems engineering (PSE) communities, and as a * Corresponding author. E-mail: [email protected]. Tel.: +34 934 016 678. Fax: +34 934 010 979.

result of this effort, significant progress has been made in problem formulations and decomposition strategies. This work focuses on the strategic level of the SCM problem and addresses the SC design problem under demand uncertainty. In simple terms, this problem involves the identification of the combination of suppliers, producers, and distributors able to provide the right mix and quantity of products and services to customers in an efficient way.4 To tackle the aforementioned problem, a rigorous multistage stochastic mathematical formulation is derived and a decomposition strategy is also proposed, aiming at overcoming the numerical difficulties associated with the resolution of the underlying large-scale mixed-integer linear problem (MILP). The paper is organized as follows. Section 2 presents a critical review of previous work on SC modeling and design. In Section 3, a formal definition of the problem of interest is given, while in Section 4, the corresponding multistage stochastic mathematical formulation is derived. The decomposition strategy aiming at the reduction of the computational expenses associated with the rigorous model is described in Section 5, and its performance is illustrated through a number of examples in Section 6. Finally, some conclusions about this work are drawn in Section 7. 2. Literature Review A large body of literature exists on SC analysis and optimization. Many of the available works are deterministic approaches5-8 that assume nominal parameter values and do not consider uncertainty. This assumption is not realistic because many SC parameters, such as cost coefficients, production rates, demand, etc., cannot be perfectly known and should be considered as sources of uncertainty. Several approaches take into account these sources of technical and commercial uncertainty arising in the area of SCM. These methodologies can be roughly classified into simulation-based and optimization-based methods according to the features of the modeling approach applied in each case. 2.1. Simulation-Based Approaches. Simulation approaches in the area of SCM under uncertainty usually work at an operational/tactical level. Part of the effort has been oriented through control theory, in which the uncertainty is modeled as disturbances arriving to a dynamic model of the system.9-11

10.1021/ie051424a CCC: $33.50 © 2006 American Chemical Society Published on Web 09/28/2006

Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7567

Multiagent-based approaches have been also applied for simulating distributed systems with somewhat decentralized decision making (especially for short-term decisions).12-14 In all the cases, the different players in the SC are represented by agents who are able to make autonomous decisions based on the information that they have available and on the messages that they receive. The agents include warehouses, customers, plants, and logistics functions. 2.2. Optimization-Based Approaches. Concerning the solution approaches available in the literature to deal with optimization under uncertainty problems, in which a sequence of decisions must be computed under an uncertain environment, one can find two major equivalent methodologies from different research streams that formulate and solve them by using a probabilistic representation of the uncertainty: stochastic optimal control and multistage stochastic programming. 2.2.1. Multistage Stochastic Programming with Recourse. Multistage stochastic programming deals with problems that involve a sequence of decisions reacting to outcomes that evolve over time. At each stage, one makes decisions based on currently available information (i.e., past observations and decisions) prior to the realizations of future events. The multistage stochastic programming has been extensively applied to solve tactical/ strategic SCM problems.15-17 Contributions belonging to this group differ primarily in the selection of the decision variables and the way in which the expected value term, which involves a multidimensional integral accounting for the probability distribution of the uncertain parameters, is computed. In view of this, two distinct methodologies for representing uncertainty can be identified within probabilistic methods. These are the scenario-based and distribution-based approaches. In the former approach,18-21 the uncertainty is described by a set of discrete scenarios capturing how the uncertainty might play out in the future. Each scenario is associated with a probability level representing the decision maker’s expectation of the occurrence of a particular scenario. The scenario-based approach avoids the problem of multivariate integration when the random variables follow multidimensional continuous distributions. This is achieved by generating a finite set of scenarios, from sampling or a discrete approximation of the given distributions, to represent the probability space. Generating scenarios to approximate the underlying problems is a research subject in itself.22,23 With the scenarios or scenario tree specified, the stochastic program becomes a deterministic equivalent program. With integer requirements, the powerful features of convexity and duality that the original problem may have are lost. In general, those decomposition schemes that are efficient for linear models can no longer be formally justified. In those cases where a natural set of discrete scenarios cannot be identified and only a continuous range of potential futures can be predicted, a distribution-based approach must be used. By assigning a probability distribution to the continuous range of potential outcomes, the need to forecast exact scenarios is obviated. The distribution-based approach is adopted in several works,24,25 in which demand is modeled as normally distributed with a specified mean and standard deviation. However, the resulting formulation is solvable in a limited number of cases. 2.2.2. Stochastic Optimal Control. Stochastic optimal control, or Markov decision process, characterizes a sequential decision problem in which the decision makers choose an action in the state occupied at any decision epoch according to a decision rule or policy. Dynamic programming provides a framework for studying such problems, as well as for devising

algorithms to compute an optimal control policy. In the area of SCM, the stochastic optimal control has been mainly applied to solve operational/tactical problems.26,27 Finally, the most preferred nonprobabilistic method to cope with the uncertainty in SCM has been fuzzy programming.28 2.3. SC Design under Uncertainty. The multistage stochastic programming approach has been the prevalent method for addressing the SC design problem under uncertainty.21,25,29,30 In many cases, the stochastic models devised to take strategic decisions in SCM apply fairly simple assumptions concerning the SC structure to be optimized. Specifically, they usually determine the number of echelons of a production/distribution network and the connectivity between adjacent echelons, but they do not consider the detail configuration of the facilities that belong to the network. In other words, they obviate the scheduling problem inherent to the batch plants embedded in the SC. Moreover, they usually apply the concept of fixed “echelons”, i.e., they assume a given fundamental structure for the network in terms of the echelons involved (e.g., suppliers, manufacturing plants, warehouses, distribution centers, and customers). Thus, a rather rigid structure is imposed on the SC.31 These assumptions allow one to formulate the SC design problem as a multistage stochastic programming problem without scenario-dependent binary variables. Thus, the deterministic equivalent of the stochastic model gives rise to a largescale MILP problem that has only nonscenario-dependent binary variables that represent the existence and type of the different nodes of the network.16,17 The absence of scenario-dependent binary variables makes it easy to exploit the decomposable structure of the problem, thus reducing the computational effort required to solve it.29,32-34 Compared with the SCM approaches presented to date, the methodology proposed in this paper focuses on developing a model and a solution strategy for a broader problem that deals with the design of chemical SCs taking into account the detailed configuration of the nodes embedded in it. Moreover, this work does not assume a fixed SC structure in terms of its number and type of echelons and also considers the demand uncertainty under which the network operates. To tackle this problem, a scenario-based multistage stochastic programming model is derived. The resulting formulation exhibits scenario-dependent binary variables for representing the scheduling decisions to be taken as the information arrives and the uncertainty unveils over time. The incorporation of scheduling constraints within the model avoids the use of fairly simple representations of capacity and provides a good estimate of the production limits of the batch plants of the network. To overcome the numerical difficulties associated with this formulation, we also investigate a decomposition strategy that is based on the combined use of genetic algorithms (GAs) and mathematical programming tools. 3. Problem Statement Production-distribution networks that can be described through the state-task-network (STN) representation35 are considered in this work. Materials are manufactured in multipurpose batch chemical plants and are stored in warehouses prior to being sent to final markets where they become available to customers. The following data is assumed to be known in advance: (i) The set of raw materials and intermediate and final products to be manufactured, stored, and transported through the nodes of the network. (ii) The set of production recipes and prices of final products. (iii) The superstructure of the network (i.e., number of available nodes and their discrete capacities). Regarding the

7568

Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006

Figure 2. Scenario tree. Figure 1. Problem statement.

batch plants, it is necessary to provide the number of available equipment units, their discrete capacities, and their suitabilities for the labor tasks. (iv) The cost functions associated with the equipment units, raw materials, utilities consumption, holding inventory over the horizon, and lost demand due to inadequate production. The SC operation under demand uncertainty for the whole time horizon is as follows (see Figure 1). Sales of products are executed at the beginning of each of the periods of time in which the time horizon is divided. The demand associated with each of these periods cannot be perfectly forecasted, and its uncertainty is represented by a set of scenarios with a given probability of occurrence. Decisions regarding scheduling tasks are made in stages as information arrives and uncertainty unveils over time. Therefore, the scheduling decisions regarding the whole SC (number of tasks to be performed, batch sizes, and assignment and sequencing decisions) must be taken at the beginning of each period, that is, prior to the realization of the uncertainty. On the other hand, sales are computed once the random events take place at the end of each period (i.e., beginning of the next one). Thus, the problem contemplated involves a sequence of decisions, which in our case are the schedules, purchases of raw materials, and sales of final products, that react to outcomes (i.e., demand realization) that evolve over time. Let us mention that, unlike other approaches in the literature, this work assumes that some of the demand can actually be left unsatisfied because of limited production capacity considerations. 3.1. Modeling Assumptions. The time horizon is assumed to be divided in |K| periods, wherein the random variable (demand) takes a finite number of possible realizations. Each set of possible outcomes in each of the given periods is called a scenario. The description of scenarios when talking about multistage stochastic programming is often made on a tree such as that in Figure 2. Here, there are |MK| scenarios that are considered to be evident in the last stage K with a probability equal to PmK. In previous stages, there are a more limited number of possible realizations, which we call the stage k scenarios and represent by mk. Mk denotes the set of scenarios till the end of period (k - 1). For example, M1 consists only of the root node (1), M2 consists of the nodes at the end of period 1 (i.e., 2 and 3), and so on. Each of these period k scenarios is said to have a single ancestor scenario in stage k - 1 and one or several descendant scenarios in stage k + 1. ADk is defined as the set of binary combinations representing all possible ancestor-descendant

combinations at stage k of the scenario tree. For the example presented in Figure 2, AD3 ) {(2, 4),(2, 5),(3, 6),(3, 7)}, as is depicted in the same figure. DSmk represents the set of descendant scenarios of scenario mk. For example, for m ) 2 and k ) 3, DSmk ) {6, 7}, as is shown in the figure. Finally, ASmk represents the ancestor scenario of scenario mk. For instance, for m ) 3 and k ) 6, ASmk ) {3}, as is also depicted in the figure. Different scenarios at stage k may correspond to the same mk realizations and are only distinguished by differences in their ancestors.36 Moreover, the root node of the scenario tree represents the initial state of the system, while the leaves of the tree are the end of the time horizon. Furthermore, every unique path from the root node to a leaf node represents one unique timeline (a unique future) along which the scheduling decisions executed could live through. Therefore, TLmK represents the set of stage k scenarios belonging to the timeline that ends in scenario mK at final stage K. For the example presented before, we have TL11 ) {1, 2, 5, 11}. The problem then is to obtain a SC design that maximizes the expected profit, which can be achieved by matching demand as closely as possible, thus reducing inventories and minimizing the labor costs, the costs for underproduction, and the investment costs. 4. Multistage Stochastic Programming Approach In a multistage stochastic optimization approach, the uncertain model parameters are considered as random variables with an associated probability distribution and the decision variables are classified into several stages. The k-stage variables correspond to those decisions made after the uncertainty associated with the stages up to k - 1 is unveiled. After the decisions up to k - 1 are taken and the random events up to k - 1 are realized, the k-stage decisions are made subject to the restrictions imposed by the k-stage problem. Birge and Louveaux36 provide an overview of this kind of stochastic technique. Having in mind the concepts previously presented, it is possible to formulate the SC design problem under uncertainty as a (|K| + 1)-stage stochastic MILP whose solution involves the computation of the optimal SC design and the scheduling decisions to be implemented on each of the nodes of the scenario tree. Three types of constraints should be considered within this formulation: the assignment, the mass balance, and the capacity constraints. The mathematical formulation next derived is based on the state-task-network (STN) representation developed by Kondili et al.35 and implies the discretization of each time period k into |T| scheduling intervals of lower length. The mass balance


equations have been modified to properly model the transport of materials between the nodes of the SC. To decouple units from tasks, we use the following rule: if a task i can be performed in two units j and j′, then two tasks i (performed in unit j) and i′ (performed in unit j′) are defined.37,38 The involved sets of equations are next described in detail. 4.1. Assignment Constraints. The basic assignment constraint 1 has been taken from the work of Shah et al.,39 in which the authors reformulated the original assignment constraint used by Kondili et al.35 t

∑

∑

i∈Ij∩NTR t′)t-pti+1

Wit′kmk e 1 ∀j ∈ J, ∀t ∈ T, ∀k < |K|, ∀mk ∈ Mk (1)

This constraint implies that unit j cannot be assigned to tasks other than i during the interval [t - pi + 1, t] in every period k and scenario mk (i.e., ensures that at most one of these suitable tasks i is assigned to unit j during the defined interval in every period and scenario). It is only applied for production equipment (i ∈ Ij ∩ NTR) and not for transport ones (i ∈ Ij ∩ TR) (i.e., trucks). For transport equipment, it is necessary to take into account both the time required to transport materials from one node to another and the time spent in coming back to the origin once the transport task has been accomplished, as is stated in eq 2. t

interval of each of the periods of time in which the time horizon is divided. For final products (s ∈ FP), Purchskmk should be removed if the possibility of outsourcing is not contemplated. Moreover, for raw materials (s ∈ RM), Salesskmk should be removed, and for intermediate products (s ∈ IP), both should be removed. Equation 4 represents the mass balance for scheduling intervals beyond the first one. In this equation, we only consider the amount of materials in state s produced in scheduling interval t, thus omitting the purchases and sales terms. Constraint 5 is suitable for the first time interval of the last period of time, for which only sales decisions are computed. Finally, eqs 6 and 7 are applied to connect each scenario in stage k + 1 with its ascendant scenario in stage k through the initial amounts of materials in state s.

Equation 8 states that the sales can be lower than or equal to the demand (i.e., part of the demand can be left unsatisfied because of limited production capacity). I ) Bitkmk‚(-FIis) B istkm k

∀i ∈ NTR, ∀s ∈ SI′i, ∀t ∈ T, ∀k < |K|, ∀mk ∈ Mk (9)

O Bistkm ) Bitkmk‚FOis k

∀i ∈ NTR, ∀s ∈ SO′i, ∀t ∈ T, ∀k < |K|, ∀mk ∈ Mk (10)

∑ ∑ ) Wit′km e 1 i∈I ∩TR t′)t-2pt +1 k

j

i

∀j ∈ J, ∀t ∈ T, ∀k < |K|, ∀mk ∈ Mk (2)

When transport tasks are carried out by external suppliers, which turns out to be a common situation, it seems reasonable to model them as an external utility consumed by the network. In this case, the assignment constraints corresponding to the transport tasks can be omitted. 4.2. Mass Balance and Utility Consumption Constraints. The model considers that scheduling decisions are taken within each period of time and are linked to adjacent periods of the scenario tree through the initial amounts of materials. Hence, the mass balances are enforced via the following constraints:

Sst0kmk + Purchskmk ) Sstkmk +

∑

i∈SIs

I Bistkm k

+ Salesskmk

∀s ∈ S, ∀k ∈ |K|, t ) 1, ∀mk ∈ Mk (3) Sst-1kmk +

∑

i∈SOs

O Bist-pt ikmk

) Sstkmk +

∑

i∈SIs

I B istkm k

∀s ∈ S, ∀k ∈ |K|, ∀t ) 1, ∀mk ∈ Mk (4)

Sst0kmk + Sstkmk + Salesskmk

∀s ∈ S, k ∈ |K|, t ) 1, ∀mk ∈ Mk (5)

S0s ) Sst0kmk

∀s ∈ S, ∀k ∈ K, ∀mk ∈ Mk (8)

Salesskmk e Demskmk

∀s ∈ S, k ) 1, ∀mk ∈ Mk

(6)

Sstk-1mk-1 ) Sst0kmk′

∀s ∈ S, ∀k > 1, t ) |T|, ∀mk-1, m′k ∈ ADk (7)

That is, in eq 3, it is stated that the initial amount of state s plus the amount purchased must equal the holdup plus the amount consumed and the sales. This equation assumes that purchases and sales of products take place in the first scheduling

Equations 9 and 10 compute the amounts of materials produced O I (continuous variable B istkm ) and consumed (B istkm ) from the k k batch size of task i started in time interval t (Bistkmk) and scenario mk of period k. These equations are valid for nontransport tasks (i ∈ NTR) with constant values of mass fractions for consumption and production of states (FIis and FOis , respectively). I B istkm ) Bitkm ∑ s∈SI ′ k

i

k

∀i ∈ TR, ∀t ∈ T, ∀k < |K|, ∀mk ∈ Mk (11)

I ) B istkm k

∑

s′∈SOi′∩TSOs

O Bis′tkm k

∀i ∈ TR, ∀s ∈ SI′i, ∀t ∈ T, ∀k < |K|, ∀mk ∈ Mk (12) On the other hand, constraints 11 and 12 are applied to model the transport tasks (i ∈ TR), which can operate over a range of mass fractions rather than for a given fixed value, as occurred with the production tasks. Therefore, eq 11 constrains the I summation over s of the continuous variable B istkm to be equal k to the batch size of the transport task i (Bitkmk) started in time interval t in scenario mk of period k. Equation 12 links the input and output states of a transport task i by enforcing the total amount of output states s′ coming from state s (s′ ∈ TSOs) to equal the input flow of s. pti-1

Uutkmk )

∑ ∑ Rui‚Wit-t′km + βui‚Bit-t′km i∈I t′)1 k

k

j

∀u ∈ U, ∀j ∈ J, ∀t < T, ∀k < |K|, ∀mk ∈ Mk (13) 0 e Uutkmk e U max uk

∀u ∈ U, ∀t ∈ T, ∀k < |K|, ∀mk ∈ Mk (14)

7570


Finally, eq 13 is applied to compute the total consumption of utility u in time interval t within period k in scenario mk. In this constraint, fixed and variable coefficients for consumption of utility u by task i are introduced (Rui and βui). The amount of utility consumed in a given period is also constrained to be lower than the maximum available one Umax uk , as is stated in eq 14. 4.3. Design Constraints. Three types of equipment are considered in the formulation, namely, production (j), single storage tanks (l), and warehouses (n). Three binary variables (Xjdj, Yldl, and Zndn, respectively) are, therefore, defined in order to represent the existence of an equipment of discrete size d. These variables take the value of 1 if the equipment being represented is of discrete size d and the value of 0 otherwise.

Witkmk e

∑d Xjd

4.4. Objective Function. The presented model accounts for the maximization of the total expected profit (eq 22). In each period, revenues are obtained through sales of final products, while costs are due to underproduction (i.e., leaving part of the demand unsatisfied), purchases of raw materials, consumption of utilities, and holding of inventories, as is stated in eq 23. The investment cost is given by eq 24

E[Profit] ) ProfitmK )

∑ ∑ ∑ s∈S k∈K m′∈TL

∑ ∑ ∑ s∈S k∈K m′∈TL k

j

k

(22)

(Demskmk′ - Salesskmk′)‚UDCosts mK

k

Purchskmk′‚RMCosts -

mK

∑∑∑ ∑

Uutkmk′‚UCostu -

u∈U t∈T k∈K m′k∈TLmK

∑ ∑∑ ∑ s∈S t∈T k∈K m′∈TL

j

∀i ∈ Ij, ∀t ∈ T, ∀j ∈ J, ∀k < |K|, ∀mk ∈ Mk (16) Equation 15 forces all the binary variables representing the starting of the tasks performed in equipment j to be equal to 0 in the case that equipment j is not selected. Moreover, eq 16 limits the batch sizes according to the capacities of the equipment. In this equation, µi is the size factor of task i and ESizejdj represents the discrete size d of equipment j.

Yld ‚STSizeld ∑ d l

- IC

Salesskmk′‚Prices -

j

Sstkmk e

K

mK

∑ ∑ ∑ s∈S k∈K m′∈TL

j

Xjd ‚ESizejd ∑ d

K

K

∀i ∈ Ij, ∀t ∈ T, ∀j ∈ J, ∀k < |K|, ∀mk ∈ Mk (15)

Bitkmk·µi e

Pm ‚Profitm ∑ m

k

IC )

Sstkmk′‚ICosts (23)

mK

Xjd ‚ECostjd + ∑∑Yld ·STCostld + ∑j ∑ d l d j

j

j

l

l

l

Znd ‚WCostnd ∑n ∑ d n

n

(24)

n

Finally, the overall problem (model SCHEDMS) can be expressed as follows:

l

l

∀s ∈ SSl, ∀t ∈ T, ∀k < |K|, ∀mk ∈ Mk, ∀l ∈ L (17)

max E[Profit] subject to constraints 1-24

Equation 17 expresses the limits for the storage tanks. The amount of state s stored in storage tank l (Sstkmk) in any time interval and scenario is bounded by the capacity of the associated storage tank. Here, STSizeldl represents the discrete size d of storage tank l. In this work, the warehouses are modeled as common storage tanks that can store simultaneously more than one product. This assumption implies that the final products are packed in the plants before being sent to the warehouses. Thus, different products can share the same warehouse without being mixed.

The multistage stochastic formulation described above leads to a large-scale MILP in which the number of continuous and binary variables as well as equations is proportional to the number of nodes of the scenario tree, which explodes exponentially with the number of periods and scenarios. In view of this, the overall problem is, in general, computationally intensive or even intractable because of the large size of the scenario tree. A special numerical technique based on decomposition is, therefore, required to solve large-scale instances of the problem.

Znd ‚WSizend ∑ Sstkm e ∑ s∈SS d

5. Decomposition Strategy

k

n

n

n

n

∀t ∈ T, ∀k < |K|, ∀mk ∈ Mk, ∀l ∈ L (18)

Constraint 18 expresses the limits for the warehouses. The total amount of states s stored in warehouse n in any time interval and scenario is bounded by its capacity. Here, WSizendk represents the discrete size d of warehouse n. Finally, the selection of a single discrete size for each equipment item can be enforced by adding the following constraints (19-21) to the model:

Xjd e 1 ∑ d

∀j ∈ J

(19)

Yld e 1 ∑ d

∀l ∈ L

(20)

Znd ∑ d

∀n ∈ N

(21)

j

j

l

l

n

n

e1

As was mentioned before, the solution of the multistage stochastic program previously described could be addressed, in principle, by means of standard mathematical programming tools. However, this would lead to high computational requirements owing to its scale and complexity, which would make its implementation in a real scenario at present rather difficult. The approximation strategy presented in this article aims at overcoming the numerical difficulties associated with the resolution of this model. The proposed algorithm is based on the combined use of GAs and mathematical programming tools. 5.1. Genetic Algorithms (GAs) in SCM. The decomposition strategy presented in this work applies GAs, which constitute an optimization method of searching based on evolutionary processes. As opposed to many other optimization methods, GA works with a population of solutions instead of just a single solution. GA assigns a value to each individual in the population according to the specific objective function of the problem to be solved. A survival-of-the-fittest step selects individuals from


Figure 3. Proposed approach.

the old population, and a reproduction step applies operators such as crossover or mutation to these individuals to produce a new population that is expected to be fitter than the previous one. The main advantage of GAs over other metaheuristics approaches (e.g., simulated annealing) is that they are capable of exploring a larger area of the solution space with a smaller number of objective function evaluations. Since, in the context of our work, evaluating the objective function entails running a large number of two-stage stochastic models, being able to find high-quality solutions early in the search space is of critical importance. GAs have been used to optimize multimodal, discontinuous, and nondifferentiable functions. In the field of SCM, these approaches are supported by the research in GA combined with mathematical programming.40-42 Research in GA for SCM focuses on solving, for instance, combinatorial operations research problems such as the multistage facility location/ allocation problem in which the decisions to be determined are the facilities to be opened or the distribution network to be used. GAs have also been implemented for the optimal design of multiproduct batch chemical process.43,44 In these approaches, the GA is usually used as the search engine that seeks optimal values for the discrete variables of the problem and uses other mathematical models as a subroutine for evaluating them. This idea is applied in this work for decomposing the multistage stochastic programming model previously presented. 5.2. Proposed Approach. The overall SC design problem is then hierarchically decomposed into two levels (see Figure 3), a higher strategic level and a lower tactical/operational one. This leads to tractable problems and avoids inflexible monolithic formulations that additionally become impossible to solve in the case of large-scale SC problems. At the strategic level, the capacities of the warehouses and plants of the network (i.e.,

Figure 4. Inner loop: Two-stage stochastic programming models computed within a rolling horizon strategy.

number and size of equipment and storage tanks) are considered. At the lower level, the tactical and operational decisions, namely, the quantities manufactured at the factories, the quantities transported between nodes, and the sales, are computed by applying a two-stage stochastic SC scheduling model based on the structure of the network being assessed. Demand uncertainty is accommodated within this framework by running this twostage stochastic scheduling model in the so-called rollinghorizon mode.45 That is, the operating horizon is divided into a certain number of periods, and the model with suitable demand forecasts is solved to yield scheduling decisions for each period; only those belonging to the first period are implemented. At the end of the first period, the state of the system, including inventory levels, is updated and the cycle is repeated with the horizon advanced by one period (see Figure 4).

7572


The overall algorithm operates as follows. At the upper level (master problem), a GA, which provides an oriented search mechanism that decreases the computational effort required by rigorous optimization techniques, is applied, thus generating potential design alternatives. At the inner level (slave problem), the operational variables associated with each design candidate are optimized under the uncertain environment by running the corresponding two-stage stochastic SC scheduling model in a rolling-horizon mode. This procedure provides the scheduling decisions to be implemented in each of the nodes of the scenario tree. The two-stage stochastic scheduling formulation (SCHEDTS) can be easily derived from the multistage stochastic model (SCHEDMS) previously presented. In this case, only one period of time should be considered. Moreover, the binary variables representing the SC design must also be fixed according to the structure of the network being assessed. In fact, the values of the binary variables are provided in each iteration by the GA. With regard to the inner level of our algorithm, let us mention that the general scheme presented in this article could be easily adapted to work with other decomposition techniques that address scheduling under uncertainty. In this work, two-stage stochastic models computed within a rolling-horizon mode have been applied to optimize the scheduling variables under demand uncertainty. However, shrinking-horizon approaches46 and deterministic rolling-horizon approaches, among others, could also be used. The fitness functions of the individuals belonging to each population generated by the GA (i.e., of the SC design alternatives computed in each iteration of the GA) are determined as follows. First, a set of two-stage stochastic SC scheduling models based on the SC configuration being analyzed are run within a rolling-horizon mode for all the nodes of the scenario tree. For each period of time and scenario, a two-stage stochastic model is thus computed to yield scheduling decisions for the period of time being considered, as is depicted in Figure 4. This formulation does not allow modifications to the original schedule once this is fixed at the beginning of the period of time for which it has been defined. Moreover, the two-stage stochastic formulation considers only one period of time and takes into account all the descendant scenarios of the node (i.e., of the scenario) for which it is computed. Thus, in Figure 4, “TS” means two-stage stochastic programming and is placed above the arrow that covers the schedule that has been computed by applying a two-stage stochastic formulation. “Fixed” means that the schedule under the arrow has been already computed in previous periods and is “frozen”. For example, in the first node of the third time period (node 4), the schedules of the first two periods are fixed while the schedule of the third period is computed by applying a two-stage stochastic formulation. This two-stage stochastic model considers one period of time (period 3) and takes into account all the descendant scenarios of node 4 (i.e., scenarios 8 and 9). Second, the expected profit associated with each configuration is computed by calculating the average value of the profits achieved over the entire range of scenarios. This expected profit is finally used to assess the SC design alternatives. It should be mentioned that the proposed decomposition strategy also suffers from the “curse of dimensionality”, as the computation effort explodes with the number of scenarios. Nevertheless, as will be discussed in the case study, it is capable of providing solutions within a few percent of the optimal ones, incurring less CPU time than the rigorous multistage stochastic formulation. Moreover, for large-size problems, our algorithm

could be applied in conjunction with specific techniques22,23 that aim to reduce the number of scenarios required to model the uncertainty. 6. Computational Results Several examples are next presented to illustrate the performance of our algorithm. All the models were solved on an AMD Athlon 3000 computer using the MIP solver of CPLEX 7.0. 6.1. Implementation. Several aspects concerning the implementation of the proposed strategy which are considered to be of interest are next described in detail. 6.1.1. Software Implementation. The outer loop of the proposed strategy (GA) has been implemented in MATLAB,47 taking advantage of the GA toolbox available in version 7.0 of the software. The inner loop, which involves the computation of the two-stage stochastic scheduling SC formulation within a rolling-horizon mode, has been implemented in GAMS48 and has been solved with CPLEX 7.0. The interfacing between both programs is carried out by means of the software library developed by Ferris.49 6.1.2. Encoding of SC Configurations. In the GA applied in this work, a chromosome represents a specific SC configuration. Each discrete size of an equipment unit is represented through an integer number that is expressed in the chromosome in a binary basis. 6.1.3. Initial Population Creation. Although several methods exist for creating the initial population, the chosen technique consists of a random string initialization. This strategy guarantees a population various enough to explore several zones of the search space. 6.1.4. Fitness Evaluation. The optimization criteria considered for fitness evaluation involves the computation of the twostage stochastic scheduling formulation within a rolling-horizon mode. 6.1.5. Selection. The selection procedure specifies how the GA chooses parents for the next generation. The function implemented in this work lays out a line in which each parent corresponds to a section of the line, the length of which is proportional to its scaled value. The algorithm moves along the line in steps of equal size. At each step, the algorithm allocates a parent from the section it lands on. The first step is a uniform random number which must be lower than the step size. 6.1.6. Crossover. Crossover specifies how the GA combines two individuals, or parents, to form a crossover child for the next gene. In this work, the crossover creates a random binary vector and selects the genes where the vector is a 1 from the first parent and the genes where the vector is a 0 from the second parent. Then, it combines both sets of genes to form the child. 6.1.7. Mutation. Mutation is the process through which the GA makes small random changes in the individuals in the population to create mutation children. Mutation provides genetic diversity and enables the GA to search a broader space. The mutation function in this work is a Gaussian function which adds a random number taken from a Gaussian distribution to each entry of the parent vector. 6.1.8. Reproduction Options. Reproduction options specify how the GA creates children for the next generation. There is a parameter which specifies the number of individuals that are guaranteed to survive to the next generation which is called the elite count. The elite count is a positive integer lower than or equal to the population size. The crossoVer fraction is the other parameter to be fixed and specifies the fraction of the next generation, excluding the elite children, that are produced by crossover. The crossover fraction must be a fraction between

Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7573 Table 1. Processing Times of Tasks i for Case Study 1 (pti (t.u.)) task

pti

task

pti

R1-RI R2-RI R3-RIIa R3-RIIb R4-RIIa

2 2 6 6 6

R4-RIIb P1-IW1 IW1-Ret.1 IW1-Ret.2

6 0 0 4

Table 2. Mass Fractions for Consumption and Production of States by Task i for Case Study 1 (FIis and FOis (adim.)) state

Figure 5. STN representation of case study 1.

task

RM1

R1-RI1 R1-RI2 R2-RI1 R2-RI2 R3-RII R4-RII

-1 -1

RM2

-1 -1

I1

I2

I3

A

B

1 1 0.1 0.1

0.9 0.9 -0.9

-1 -0.1

1 1

Table 3. Price of State s for Case Study 1 (Prices (m.u./kg)) state

Prices

state

Prices

state

Prices

RM1 RM2

1500 2500

A(Ret.1) A(Ret.2)

3500 3500

B(Ret.1) B(Ret.2)

4500 4500

Table 4. Inventory Cost of State s for Case Study 1 (ICosts (m.u./kg‚t.u.)) state

ICosts

state

ICosts

state

ICosts

RM1 RM2 I1 I2 I3

7.5 12.5 12.5 12.5 7.5

A(PW1) B(PW1) A(IW1) B(IW1) A(Ret.1)

17.5 22.5 17.5 22.5 17.5

B(Ret.1) A(Ret.2) B(Ret.2)

22.5 17.5 22.5

Table 5. Penalization for Demand of State s Unsatisfied for Case Study 1 (UDCosts (m.u./kg)) Figure 6. SC superstructure of case study 1.

0 and 1. For example, if the population size is 20, the elite count is 2, and the crossover fraction is 0.8, the numbers of each type of child in the next generation is as follows: there are 2 elite children and 18 individuals other than elite children, so the algorithm rounds 0.8 × 18 ) 14.4 to 14 to get the number of crossover children. The remaining four individuals, other than elite children, are mutation children. 6.1.9. Stop Criterion. The stop criterion considered in this study is the maximum number of generations (maxgen) to reach. 6.2. Example 1a. In this example, we investigate the retrofit of an existing SC comprising one plant, one warehouse, and two retailers. The plant embedded in the SC (P1) manufactures two different products (A and B), which are stored in the final product warehouse of the plant (PW1) before being sent to an intermediate warehouse (IW1), from where they are transported to the retailers (Retailer1 and Retailer2) in which they become available to customers. The STN representation of this example is given in Figure 5, while the associated superstructure is given in Figure 6. The structure of the plant has been adapted from the case study introduced by Maravelias and Grossmann.50 There are two types of reactors available for the process (types I and II) with different numbers of corresponding units available: one reactor (RI) of type I and two reactors (RIIa and RIIb) of type II. Reactions R1 and R2 require a type I reactor, whereas reactions R3 and R4 require a type II reactor. The plant operates under a FIS policy (finite intermediate storage) for all the states. Four possible discrete sizes are available for all the reactors and storage tanks and also for the warehouses. Discrete size 1 represents the actual state of the network (i.e., of the equipment units, storage tanks, and warehouses), while discrete sizes 2, 3,

state

UDCosts

state

UDCosts

A(Ret.1) A(Ret.2)

700 700

B(Ret.1) B(Ret.2)

900 900

and 4 represent the option of increasing the actual size up to a certain level. The discrete capacities of the equipment as well as the costs associated with each of them are also given. The size factor is equal to 1 L/kg for all the tasks. The rest of the data of the case study can be found in Tables 1-7. The initial inventories equal zero for all the states. Two different types of utilities are considered: the first one, which is associated with the labor tasks, has a unitary cost of 0.25 m.u., while the second one, which corresponds to the transport services, costs 0.1 m.u./ unit. The fixed and variable coefficients for consumption of utilities (Rui and βui) are equal to 2000 units of utility/time unit and 20 units of utility/kg‚time unit, respectively, for all the labor tasks, and 50 units of utility/time unit and 5000 units of utility/ kg‚time unit, respectively, for the transport ones. The availability of both utilities is assumed to be unlimited. This example considers a time horizon of 40 time units, divided into two equally long time periods of 20 time units. Each time period contains 10 scheduling intervals of 2 time units each. Two possible events in each period of time (i.e., two possible values for the demand) are considered. This leads to four scenarios () 22) for the whole time horizon. The multistage stochastic formulation is able to find the optimal solution of the problem (i.e., a solution with a 0% optimality gap) in 514 CPU s. This solution has an expected profit of 3 662 794 m.u. The GA is able to compute the same solution in 665 CPU s. Let us note that the parameters of the GA, mainly the crossover fraction and the population size, have been tuned experimentally. In this case, we have fixed a

7574


Table 6. Discrete Sizes and Costs of Equipment Units and Warehouses for Case Study 1 ESizejdj (1)

ECostjdj (m.u.)

equipment

1

2

3

4

1

2

3

4

RI RIIa RIIb

500 200 200

1000 400 400

1500 600 600

2000 800 800

0 0 0

1 577 393 910 282 910 282

2 011 846 1 160 996 1 160 996

2 390 881 1 379 730 1 379 730

STSizeldl (1) storage tanks RM1 RM2 I1 I2 I3

STCostldl (m.u.)

1

2

3

4

1

2

3

4

500 500 200 100 100

1000 1000 400 200 200

1500 1500 600 300 300

2000 2000 800 400 400

0 0 0 0 0

946 436 946 436 546 169 360 337 360 337

1 207 108 1 207 108 696 598 459 583 459 583

1 434 529 1 434 529 827 838 546 169 546 169

WSizendn (1)

WCostndn (m.u.)

warehouses

1

2

3

4

1

2

3

4

plant (PW1) intermediate (IW1) Retailer1 Retailer2

1000 1000 1000 1000

2000 2000 2000 2000

3000 3000 3000 3000

4000 4000 4000 4000

0 0 0 0

239 088 239 088 239 088 239 088

304 939 304 939 304 939 304 939

362 390 362 390 362 390 362 390

Table 7. Demand of State s for Case Study 1A (Demandskmk (kg)) state

Demandskmk

state

k)1 A(Ret.1) A(Ret.2) B(Ret.1) B(Ret.2)

Demandskmk k)3

186 182 320 372

A(Ret.1) A(Ret.2) B(Ret.1) B(Ret.2)

(193, 280) (286, 289) (362, 489) (257, 332)

k)2 A(Ret.1) A(Ret.2) B(Ret.1) B(Ret.2)

(151, 360) (139, 392) (202, 478) (282, 341)

population size ) 8 and a crossover fraction ) 0.4. Figure 7 shows the encoding procedure applied in this example. Each size of an equipment unit, storage tank, and warehouse is encoded in two gens. For example, in the SC configuration given in the figure, the size of the tank that stores RM1 takes the value of the discrete size 1.

Figure 7. Encoding for case study 1.

In Figure 8, the expected profit value of the best solutions computed by the GA and the rigorous model have been plotted against the computation time. The figure shows how the points corresponding to the GA lie entirely bellow those associated with the rigorous model. In other words, for the same CPU times, the solutions computed by the multistage stochastic programming model are better than those determined by the GA (i.e., for the same CPU time, the multistage stochastic model provides a solution with an optimality gap lower than that given by the GA). The optimal decision of this example consists of keeping the original structure of the SC (i.e., not undertaking any retrofit activity in the network). 6.3. Example 1b. In this case study, the number of time periods and scenarios is increased. The aim is to demonstrate that the multistage stochastic formulation is likely to fail as the complexity of the model, which is mainly given by its number of constraints and binary and continuous variables, increases. Thus, this new case study assumes a time horizon of 90 time


Figure 8. Comparison between the multistage stochastic programming model and the proposed approach for case study 1a.

Figure 9. Comparison between the multistage stochastic programming model and the proposed approach for case study 1b.

Table 8. Demand of State s for Case Study 1b (Demandskmk (kg))

Table 9. Comparison for Case Study 1b

state

Demandskmk

state

Demandskmk


k)1 571 645 752 784


k)3 (536, 567, 668, 676) (554, 557, 588, 627) (570, 653, 796, 876) (633, 675, 791, 823)


k)2 (561, 568, 714, 735) (557, 661, 702, 767) (770, 812, 877, 928) (707, 791, 822, 957)


k)4 (576, 659, 713, 717) (626, 704, 726, 751) (691, 846, 883, 1014) (696, 713, 857, 921)

units, divided into three equally long time periods of 30 time units. Each time period contains 15 scheduling intervals of 2 time units each. Four possible events in each period of time (i.e., four possible values for the demand) are considered. This leads to 64 scenarios ()43) for the whole time horizon. The new demand is given in Table 8. As was expected, the computation time associated with both approaches dramatically increases in comparison with the previous example. In this case, the multistage stochastic formulation is not able to find the optimal solution of the problem after 50 000 s of computation time. Specifically, the best solution obtained by the rigorous approach is equal to 4 805 500 m.u., still within 58.0% of the best relaxed solution computed by the branch and bound when the time limit of 50 000 seconds was exceeded (in this case, 11 436 000 m.u.). On the other hand, the GA computes a solution equal to 9 935 466 m.u. for the same CPU time. This solution is within 13.1% of the best relaxed solution. Here, we have fixed a population size ) 10 and a crossover fraction ) 0.4. Again, in Figure 9, the expected profit value of the best solutions computed by the GA and the rigorous model have been plotted against the computation time. Let us notice that the points corresponding to the GA lie entirely above those given by the rigorous model. Thus, for CPU times below 50 000 s, the GA performs better than the multistage stochastic programming model (i.e., it provides better solutions than the rigorous model for the same computation times). Each solution of the SC design problem includes both an SC design itself and a set of scheduling decisions associated with its operation under uncertainty. Let us note that, since nonoptimal solutions are being compared (i.e., solutions with optimality gaps > 0%), the optimality of either their scheduling

approach

E[Profit]a (m.u.)

gapb (%)

E[Profit]c (m.u.)

relative errord (%)

genetic algorithm multistage model

9 935 466 4 805 500

13.1 58.0

10 138 000 10 137 000

2.0 52.6

a Expected profit of the solution computed by each approach in 50 000 CPU s. b Optimality gap of the solution computed by each approach in 50 000 CPU s. c Expected profit of the SC design computed by each approach and evaluated through the multistage stochastic formulation. d Relative error between the expected profit of the solution with optimal scheduling decisions and the solution with nonoptimal scheduling decisions.

decisions or their SC designs cannot be guaranteed. Moreover, the scheduling decisions of a nonoptimal solution have not been demonstrated to be the best ones for the SC design of the same solution. Under these circumstances, it is necessary to compute the optimal scheduling decisions associated with each SC design to perform a fair comparison between them. This computation can provide an accurate performance measure of their operation under uncertainty, thus facilitating the direct comparison between them. The optimal scheduling variables associated with a given design can be computed by fixing in the multistage stochastic formulation the binary variables representing the SC structure and optimizing the remaining decision variables. Table 9 shows the computational results associated with the evaluation of the SC designs provided by both approaches, the structures of which are given in Table 10, through the abovementioned procedure. The multistage stochastic formulations have been solved in both cases with an optimality gap of 0.25%. The numerical results show that, although the expected profit value of the best solution computed by the GA in 50 000 seconds has an optimality gap that is 45% lower than the one associated with the multistage stochastic formulation, both designs yield similar expected profit values when they are evaluated under uncertainty through the multistage stochastic formulation (10 138 000 m.u. for the GA and 10 137 000 m.u. for the multistage stochastic model). Results also indicate that our proposed approach provides a “good” estimate of the SC performance under uncertainty. In fact, the two-stage stochastic programming models run within a rolling-horizon strategy provide an expected profit within a few percent of the optimal one. In this case, the SC design computed by the proposed approach yields an expected profit of 9 935 466 m.u. when it is evaluated through the two-stage stochastic models run within

7576


Table 10. SC Configurations for Case Study 1b equipment units

storage tanks

warehouses

approach

RI

RIIa

RIIb

RM1

RM2

I1

I2

I3

plant (PW1)

int. (IW1)

Ret. 1

Ret. 2


1 1

1 1

3 4

2 2

3 4

3 3

1 1

1 1

1 1

1 1

2 2

3 2

Table 11. Processing Times of Tasks i for Case Study 2 (pti (t.u.)) task

pti

task

pti

task

pti

task

pti

Heating-H(P1) Reaction1-RI(P1) Reaction1-RII(P1) Reaction2-RI(P1) Reaction2-RII(P1) Reaction3-RI(P1) Reaction3-RII(P1) Separation-S(P1) Heating-H(P2) Reaction1-RI(P2)

2 2 2 6 6 2 2 6 4 2

Reaction1-RII(P2) Reaction2-RI(P2) Reaction2-RII(P2) Reaction3-RI(P2) Reaction3-RII(P2) Separation-S(P2) PW1-IW1 PW1-IW2 PW1-IW3 PW2-IW1

2 2 2 2 2 6 0 4 2 4

PW2-IW2 PW2-IW3 IW1-Ret.1 IW1-Ret.2 IW1-Ret.3 IW1-Ret.4 IW1-Ret.5 IW1-Ret.6 IW2-Ret.1 IW2-Ret.2

0 2 0 4 2 4 4 6 4 0

IW2-Ret.3 IW2-Ret.4 IW2-Ret.5 IW2-Ret.6 IW3-Ret.1 IW3-Ret.2 IW3-Ret.3 IW3-Ret.4 IW3-Ret.5 IW3-Ret.6

2 4 2 4 2 2 0 2 2 2

a rolling-horizon strategy and an expected profit of 10 138 000 m.u. when it is assessed through the multistage stochastic formulation, which represents a relative error of 2.0%. On the

Figure 10. SC structure of case study 2.

Figure 11. STN representation of case study 2 (I).

other hand, although the branch and bound technique applied to the multistage stochastic model provides “good” SC designs, it also leads to poor scheduling decision variables and, therefore, to poor estimates of the SC performance under uncertainty. In this example, the relative error of the scheduling decisions associated with the best SC design computed by the multistage stochastic model is 52.6%. 6.4. Example 2. Here we investigate the retrofit under uncertainty of a SC comprising several plants, warehouses, and retailers distributed in different locations. The plants embedded in the SC (P1 and P2) manufacture two different products (F and G), which are stored in the final product warehouses of the plants (PW1 and PW2) before being sent to three warehouses (IW1 to IW3), from where they are transported to six retailers (Retailer1-Retailer6) in which they become available to customers. The structure of the SC under study is given in Figure 10, while its STN representation is depicted in Figures 11 and 12. The SC superstructure is given in Figure 13. The associated data are listed in Tables 11-17. The plants of the network are assumed to be multipurpose batch chemical plants with a structure adapted from the case study proposed by Kondili et al.35 The initial inventories equal zero for all the states. Four possible discrete sizes are available for all the equipment units


Figure 12. STN representation of case study 2 (II). Table 13. Price of State s (Prices for Case Study 2 (m.u./kg)) state

Prices

state

Prices

state

Prices

F(Ret.1) F(Ret.2) F(Ret.3) F(Ret.4) F(Ret.5) A(P2)

4400 4400 4000 3800 4600 200

F(Ret.6) G(Ret.1) G(Ret.2) G(Ret.3) G(Ret.4) B(P2)

5000 6000 6000 7000 7500 400

G(Ret.5) G(Ret.6) A(P1) B(P1) C(P1) C(P2)

6000 7000 240 440 210 200

Table 14. Inventory Cost of State s for Case Study 2 (ICosts (m.u./kg‚t.u.))

Figure 13. Superstructure of case study 2. Table 12. Mass Fractions for Consumption and Production of States by Task i for Case Study 2 (FIis and FOis (adim.))

state

ICst*

state

ICst*

state

ICst*

state

ICst*

A(P1) B(P1) C(P1) D(P1) E(P1) F(PW1) G(PW1) HA(P1) I(P1)

1.2 2.2 10.5 10.5 21.5 21.0 30.5 1.2 1.6

A(P2) B(P2) C(P2) D(P2) E(P2) F(PW2) G(PW2) HA(P2) I(P2)

1.0 2.0 10.0 10.0 21.0 20.0 30.0 1.0 1.5

F(IW1) F(IW2) F(IW3) G(IW1) G(IW2) G(IW3) F(Ret.1) F(Ret.2) F(Ret.3)

20.0 20.0 20.0 80.0 80.0 80.0 22.0 22.0 20.0

F(Ret.4) F(Ret.5) F(Ret.6) G(Ret.1) G(Ret.2) G(Ret.3) G(Ret.4) G(Ret.5) G(Ret.6)

19.0 23.0 25.0 30.0 30.0 35.0 37.5 30.0 35.0

state task

A

B

C

D

E

F

G

HA

I

Heating-H -1 1 Reaction1-RI -0.9 -0.1 1 Reaction1-RII -0.9 -0.1 1 Reaction2-RI -0.6 0.9 -0.4 0.1 Reaction2-RII -0.6 0.9 -0.4 0.1 Reaction3-RI -0.9 1 -0.1 Reaction3-RII -0.9 1 -0.1 Separation-S -1 0.9 0.1

and storage tanks. Discrete size 1 represents the actual state of the network (i.e., of the equipment units, storage tanks, and warehouses), while discrete sizes 2, 3, and 4 represent the option of increasing the actual size up to a certain level. The capacities of the warehouses (3000 L) are already fixed and cannot be modified. Two different types of utilities are considered. The first one, which is associated with the labor tasks, has a cost of

0.25 m.u./unit of utility consumed, while the second one, which corresponds to the transport services, costs 0.1 m.u./unit. The fixed and variable coefficients for consumption of utility by tasks (Rui and βui) are equal to 2000 units of utility/time unit and 20 units of utility/kg‚time unit, respectively, for all the labor tasks and 50 units of utility/time unit and 5000 units of utility/ kg‚time unit, respectively, for transport ones. The availability of both utilities as well as the capacities of the storage of raw materials, intermediate, and final products are assumed to be unlimited. This example assumes a time horizon of 72 time units, divided into three equally long time periods of 24 time units. Each time period contains 12 scheduling intervals of 2 time units each. Again, four possible events in each period of time (i.e., four possible values for the demand) are considered. This leads to 64 scenarios ()43) for the whole time horizon.

7578


Table 15. Penalization for Demand of State s Unsatisfied for Case Study 2 (UDCosts (m.u./kg)) state

UDCosts

state

UDCosts

state

UDCosts

state

UDCosts

F(Ret.1) F(Ret.2) F(Ret.3)

44 44 40

F(Ret.4) F(Ret.5) F(Ret.6)

38 56 50

G(Ret.1) G(Ret.2) G(Ret.3)

60 60 70

G(Ret.4) G(Ret.5) G(Ret.6)

75 60 70

Table 16. Demand of State s for Case Study 2 (Demandskmk (kg)) Demandskmk (kg) state

k)1

k)2

k)3

k)4

F(Ret.1) G(Ret.1) F(Ret.2) G(Ret.2) F(Ret.3) G(Ret.3) F(Ret.4) G(Ret.4) F(Ret.5) G(Ret.5) F(Ret.6) G(Ret.6)

158 74 81 60 80 65 100 157 267 231 163 259

(254, 294, 252, 287) (142, 154, 124, 146) (138, 155, 108, 169) (144, 132, 125, 143) (113, 143, 122, 138) (110, 141, 123, 137) (176, 220, 171, 230) (351, 479, 339, 486) (359, 428, 374, 416) (414, 451, 380, 472) (364, 423, 357, 481) (354, 461, 414, 429)

(112, 182, 137, 199) (75, 106, 75, 115) (78, 105, 73, 93) (92, 95, 77, 132) (90, 92, 69, 98) (84, 116, 70, 101) (128, 167, 127, 150) (236, 311, 251, 313) (257, 304, 239, 289) (257, 269, 238, 288) (264, 310, 239, 335) (240, 284, 215, 311)

(247, 281, 246, 290) (133, 151, 115, 140) (120, 141, 112, 170) (129, 141, 110, 159) (119, 157, 116, 161) (127, 166, 120, 143) (151, 259, 177, 237) (390, 432, 371, 487) (394, 444, 351, 476) (337, 422, 368, 408) (354, 447, 358, 513) (330, 443, 345, 471)

Table 17. Discrete Sizes and Costs of Equipment Units and Warehouses for Case Study 2 ESizejdj (1)

ECostjdj (m.u.)

equipment

1

2

3

4

1

2

3

4

H(P1) RI(P1) RII(P1) S(P1) H(P2) RI(P2) RII(P2) S(P2)

75 100 100 75 75 100 100 75

150 200 200 150 150 200 200 150

225 300 300 225 225 300 300 225

300 400 400 300 300 400 400 300

0 0 0 0 0 0 0 0

505 353 600 562 600 562 505 353 505 353 600 562 600 562 505 353

644 539 765 972 765 972 644 539 644 539 765 972 765 972 644 539

765 972 910 282 910 282 765 972 765 972 910 282 910 282 765 972

STSizeldl (1)

STCostldl (m.u.)

storage tanks

1

2

3

4

1

2

3

4

A(P1) B(P1) C(P1) HA(P1) I(P1) E(P1) D(P1) A(P2) B(P2) C(P2) HA(P2) I(P2) E(P2) D(P2)

100 100 100 25 25 25 25 100 100 100 25 25 25 25

200 200 200 50 50 50 50 200 200 200 50 50 50 50

300 300 300 75 75 75 75 300 300 300 75 75 75 75

400 400 400 100 100 100 100 400 400 400 100 100 100 100

0 0 0 0 0 0 0 0 0 0 0 0 0 0

480 450 480 450 480 450 209 128 209 128 209 128 209 128 480 450 480 450 480 450 209 128 209 128 209 128 209 128

612 777 612 777 612 777 266 727 266 727 266 727 266 727 612 777 612 777 612 777 266 727 266 727 266 727 266 727

728 226 728 226 728 226 316 979 316 979 316 979 316 979 728 226 728 226 728 226 316 979 316 979 316 979 316 979

In this example, the best solution computed by the multistage stochastic formulation in 50 000 s is equal to 12 478 622 m.u., still within 35.2% of the best relaxed solution (in our case, 19 269 000 m.u.) computed by the branch and bound when the time limit of 50 000 s was exceeded. The GA computes a solution of 17 060 664 m.u. for the same CPU time, which is within 11.5% of the best relaxed solution computed by the branch and bound. The population size ) 15, and the crossover fraction ) 0.5. Again, Figure 14 shows the best solutions computed by both approaches over time. As occurred in the previous case, the points corresponding to the GA lie entirely above those computed by the rigorous model. Table 18 shows the computational results associated with the evaluation of the SC designs computed by each approach, the structures of which are given in Table 19, through the multistage stochastic formulation. The

Figure 14. Comparison between the multistage stochastic programming model and the proposed approach for case study 2. Table 18. Comparison for Case Study 2 approach

E[Profit]a (m.u.)

gapb (%)

E[Profit]c (m.u.)

relative errord (%)


17 060 664 12 478 622

11.5 35.2

17 252 000 17 574 000

1.1 29.0

a Expected profit of the solution computed by each approach in 50 000 CPU s. b Optimality gap of the solution computed by each approach in 50 000 CPU s. c Expected profit of the SC design computed by each approach and evaluated through the multistage stochastic formulation d Relative error between the expected profit of the solution with optimal scheduling decisions and the solution with nonoptimal scheduling decisions.

obtained results are very similar to the ones of case study 1b, and similar conclusions can, thus, be derived. Again, when the SC design provided by the rigorous approach is evaluated through the multistage stochastic formulation, a large improvement in the expected profit value of the solution is achieved. In fact, the SC designs computed by both approaches yield similar expected profit values when they are assessed through the stochastic model (17 252 000 m.u. for the GA and 17 574 000 m.u. for the multistage stochastic model). The relative errors associated with the scheduling decisions computed by both approaches are 1.1 and 29.0%, respectively.

Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7579 Table 19. SC Configurations for Case Study 2 equipment units

storage tanks

approach

H(P1)

RI(P1)

RII(P1)

S(P1)

H(P2)

RI(P2)

RII(P2)

S(P2)

A(P1)

B(P1)

C(P1)


1 2

1 4

4 1

4 4

1 1

3 1

1 3

4 3

1 1

1 1

4 4

storage tanks approach

HA(P1)

I(P1)

E(P1)

D(P1)

A(P2)

B(P2)

C(P2)

HA(P2)

I(P2)

E(P2)

D(P2)


1 1

4 4

1 1

4 4

1 1

3 2

4 4

1 1

4 4

4 1

4 1

7. Conclusions

Sets

In this work, a multistage stochastic formulation has been derived to address the design of chemical SCs under demand uncertainty. To overcome the numerical difficulties associated with the resulting large-scale stochastic MILP, a decomposition technique has been presented. This strategy combines GAs and mathematical programming tools. At the upper level (master problem), a GA is implemented for managing the design problem and proposing potential SC structures. At the inner level (slave problem), a two-stage stochastic SC scheduling model is computed within a rolling-horizon strategy to yield scheduling decisions for all the nodes of the scenario tree and provide an estimate of the expected profit achieved under demand uncertainty. The main advantages of the proposed approach have been highlighted through several case studies where a comparison with the underlying multistage stochastic formulation has been carried out. Results indicate that the proposed strategy provides solutions within a few percent of the optimal ones and reduces the computation time required by the multistage stochastic formulation. Moreover, the multistage stochastic formulation has been shown from numerical examples to be able to compute “good values” for the design variables but “poor values” for the scheduling ones. The performance of the proposed approach and the multistage stochastic model may well differ depending upon the specifics of the case studies being analyzed. In any case, the proposed approach seems to offer better results as the complexity of the problem (i.e., complexity of the SC superstructure and number of time periods and scenarios of the model) grows and the computation effort required to solve the underlying large-scale MILP explodes.

ADk ) binary tuples representing all possible ancestordescendant combinations at stage k of the scenario tree ASmk ) set of ancestor scenarios of mk DSmk ) set of descendant scenarios of mk FP ) set of states corresponding to final products Ij ) set of tasks that can be performed in equipment j IP ) set of states corresponding to intermediate products NTR ) set of nontransport tasks (production tasks) RM ) set of states corresponding to raw materials SIs ) set of tasks receiving material from state s SI′i ) set of states consumed by task i SOs ) set of tasks producing material for state s SO′i ) set of states produced by task i SSn ) set of states that can be stored in warehouse n SSl ) set of states that can be stored in storage tank l TLmK ) set of stage k scenarios belonging to the time line which ends in scenario mK at final stage |K| TR ) set of transport tasks TSOs ) set of states s′ coming from state s

Acknowledgment Financial support received from the Spanish Ministerio de Educacioń, Cultura y Deporte (FPU programs), Generalitat de Catalunya (FI programs), and European Community (project PRISM, MRTN-CT-2004-512233) is fully appreciated. Besides, financial support from GICASA-D (I0353) and OCCASION (DPI2002-00856) projects is gratefully acknowledged. Nomenclature Indices d ) discrete sizes i ) tasks j ) equipment k ) periods of time l ) storage tanks mk ) scenarios in period k n ) warehouses s ) states t ) scheduling intervals u ) utilities

Parameters Demskmk ) demand of material in state s in period k in scenario mk ESizejdj ) discrete size d of equipment j ECostjdj ) cost of the discrete size d of equipment j ICosts ) inventory cost of s Prices ) price of material in state s pti ) processing time of task i RMCosts ) cost of raw material in state s S0s ) initial amount of material in state s STSizeldl ) discrete size d of storage tank l STCostldl ) cost of the discrete size d of storage tank l Umax uk ) maximum consumption of utility u in period k UCostu ) utility cost of u UDCosts ) cost due to unsatisfied demand of material in state s WSizendn ) discrete size d of warehouse n WCostndn ) cost of the discrete size d of warehouse n Variables Bistkmk ) batch size of task i started in time interval t at period k in scenario mk I Bistkm ) amount of material in state s consumed by task i k started in time interval t at period k in scenario mk O Bistkm ) amount of material in state s produced by task i k started in time interval t at period k in scenario mk E[Profit] ) expected profit IC ) total inventory cost ProfitmK ) profit of scenario mK

7580


Purchskmk ) amount of material in state s purchased at period k in scenario mk Sst0kmk ) initial amount of material in state s at the beginning of period k in scenario mk Sstkmk ) amount of material in state s at the end of time interval t at period k in scenario mk Salesskmk ) amount of material in state s sold at period k in scenario mk Uutkmk ) consumption of utility u in scheduling interval t at period k in scenario mk Wistkmk ) binary variable (1 if task i is started in interval t at period k in scenario mk, 0 otherwise) Xjdi ) binary variable (1 if equipment unit j is of discrete size d, 0 otherwise) Yldl ) binary variable (1 if storage tank l is of discrete size d, 0 otherwise) Zndn ) binary variable (1 if warehouse n is of discrete size d, 0 otherwise) Greek Symbols Rui ) fixed coefficient for consumption of utility u by task i βui ) variable coefficient for consumption of utility u by task i µi ) size factor of task i FIis ) proportion input to task i from state s FOis ) proportion output of task i for state s PmK ) probability of scenario mK Scheduling Models SCHEDMS ) multistage stochastic scheduling formulation SCHEDTS ) two-stage stochastic scheduling formulation Literature Cited (1) Simchi-Levi, D.; Kamisky, P.; Simchi-Levi, E. Designing and Managing the Supply Chain. Concepts, Strategies, and Case Studies; Irwin McGraw-Hill: New York, 2000. (2) Fox, M. S.; Barbuceanu, M.; Teigen, R. Agent-oriented supply chain management. Int. J. Flexible Manuf. Syst. 2000, 12, 165-188. (3) Applequist, G. E.; Pekny, J. F.; Reklaitis, G. V. Risk and uncertainty in managing manufacturing supply chains. Comput. Chem. Eng. 2000, 24 (9-10), 2211-2222. (4) Talluri, S.; Baker, R. C. A multiphase mathematical programming approach for effective supply chain design. Eur. J. Oper. Res. 2002, 141, 544-558. (5) Geoffrion, A. M.; Graves, G. Multicommodity distribution system design by Benders Decomposition. Manage. Sci. 1974, 20, 822-844. (6) Brown, G. G.; Graves, G. W.; Honczarenko, M. Design and operation of a multicommodity production/distribution system using primal goal decomposition. Manage. Sci. 1987, 33, 1469-1480. (7) Yan, H.; Yu, Z.; Cheng, T. E. A strategic model for supply chain design with logical constraints: formulation and solution. Comput. Oper. Res. 2003, 30, 2135-2155. (8) Guilleń, G.; Badell, B.; Espunã, A.; Puigjaner, L. Simultaneous optimization of process operations and financial decisions to enhance the integrated planning/scheduling of chemical supply chains. Comput. Chem. Eng. 2006, 30, 421-436. (9) Bose, S.; Pekny, J. F. A model predictive framework for planning and scheduling problems: A case study of consumer goods supply chain. Comput. Chem. Eng. 2000, 24, 329-335. (10) Perea-Lo´pez, E.; Grossmann, I.; Ydstie, E.; Tahmassebi, T. Dynamical modeling and classical control theory for supply chain management. Comput. Chem. Eng. 2000, 24, 1143-1149. (11) Seferlis, P.; Giannelos, N. A two-layered optimization-based control strategy for multi-echelon supply chain networks. Comput. Chem. Eng. 2004, 28, 799-809. (12) Julka, N.; Karimi, I.; Srinivasan, R. Agent-Based Supply Chain Management. 1: Framework. Comput. Chem. Eng. 2002, 26, 1755-1769. (13) Guilleń, G.; Mele, F. D.; Urbano, F.; Espunã, A.; Puigjaner, L. An agent-based approach for supply chain retrofitting under uncertainty. In European Symposium on Computer Aided Process Engineering-15 (20 B); Puigjaner, L., Espunã, A., Eds.; Elsevier: Amsterdam, 2005.

(14) Mele, F. D.; Guilleń, G.; Espunã, A.; Puigjaner, L. A simulationbased optimization framework for parameter optimization of supply chain networks. Ind. Eng. Chem. Res. 2006, 45, 3133-3148. (15) Subrahmanyam, S.; Pekny, J. F.; Reklaitis, G. V. Design of batch chemical plants under market uncertainty. Ind. Eng. Chem. Res. 1994, 33, 2688. (16) Tsiakis, P.; Shah, N.; Pantelides, C. C. Design of multi-echelon supply chain networks under demand uncertainty. Ind. Eng. Chem. Res. 2001, 40, 3585-3604. (17) Guilleń, G.; Mele, F.; Bagajewicz, M.; Espunã, A.; Puigjaner, L. Multiobjective supply chain design under uncertainty. Chem. Eng. Sci. 2005, 60, 1535-1553. (18) Bonfill, A.; Bagajewicz, M.; Espunã, A.; Puigjaner, L. Risk management in scheduling of batch plants under uncertain market demand. Ind. Eng. Chem. Res. 2004, 43, 741-750. (19) Bonfill, A.; Espunã, A.; Puigjaner, L. Addressing robustness in scheduling batch processes with uncertain operation times. Ind. Eng. Chem. Res. 2005, 44, 1524-1534. (20) Guilleń, G.; Pina, C.; Espunã, A.; Puigjaner, L. Optimal Offer Proposal Policy in an Integrated Supply Chain Management Environment. Ind. Eng. Chem. Res. 2005, 44, 7405-7419. (21) Guilleń, G.; Bagajewicz, M.; Sequeira, S. E.; Espunã, A.; Puigjaner, L. Management of pricing policies and financial risk as a key element for short-term scheduling optimization. Ind. Eng. Chem. Res. 2005, 44, 557575. (22) Dupacova, J.; Consigli, G.; Wallace, S. W. Generating scenarios for multi-stage stochastic programs. Ann. Oper. Res. 2000, 100, 25-53. (23) Hoyland, K.; Wallace, S. W. Generating scenario trees for multistage decision problems. Manage. Sci. 2001, 47 (2), 295-307. (24) Gupta, A.; Maranas, C. D. A two-stage modeling and solution framework for multisite midterm planning under demand uncertainty. Ind. Eng. Chem. Res. 2000, 39, 3799-3813. (25) Gupta, A.; Maranas, C. D. Managing demand uncertainty in supply chain planning. Comput. Chem. Eng. 2003, 27, 1219-1227. (26) Cheng, L.; Subrahmanian, E.; Westerberg, A. W. Multiobjective decisions on capacity planning and inventory control. Ind. Eng. Chem. Res. 2004, 43, 2192-2208. (27) Chenga, L.; Duran, M. A. Logistics for worldwide crude oil transportation using discrete event simulation and optimal control. Comput. Chem. Eng. 2004, 28, 897-911. (28) Sakawa, M.; Nishizaki, I.; Uemura, Y. Fuzzy programming and profit and cost allocation for a production and transportation problem. Eur. J. Oper. Res. 2001, 131, 1-15. (29) Iyer, R. R.; Grossmann, I. E. A bilevel decomposition algorithm for long-range planning of process networks. Ind. Eng. Chem. Res. 1998, 37, 474-481. (30) Levis, 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. (31) Shah, N. Process industry supply chains: Advances and challenges. Comput. Chem. Eng. 2005, 29, 1225-1235. (32) Liu, M. L.; Sahinidis, N. V. Optimization in process planning under uncertainty. Ind. Eng. Chem. Res. 1996, 35, 4154-4165. (33) Ahmed, S.; Sahinidis, N. V. Analytical investigations of the process planning problem. Comput. Chem. Eng. 2000, 23 (11-12), 1605-1621. (34) Verweij, B.; Ahmed, S.; Kleywegt, A. J.; Nemhauser, G.; Shapiro, A. The sample average approximation method applied to stochastic routing problems: A computational study. Comput. Appl. Optim. 2003, 24, 289333. (35) Kondili, E.; Pantelides, C. C.; Sargent, R. W. A general algorithm for short-term scheduling of batch operations. Comput. Chem. Eng. 1993, 17, 211-227. (36) Birge, Z.; Louveaux, S. Principles on Stochastic Programming; Springer-Verlag: New York, 1997. (37) Ierapetritou, M. G.; Floudas, C. A. Effective continuous-time formulation for short-term scheduling. Part 1. Multipurpose batch processes. Ind. Eng. Chem. Res. 1998, 37, 4341-4359. (38) Ierapetritou, M. G.; Floudas, C. A. Effective continuous-time formulation for short-term scheduling. Part 2. Continuous and semicontinuous processes. Ind. Eng. Chem. Res. 1998, 37, 4360-4374. (39) Shah, N.; Pantelides, C. C.; Sargent, R. W. H. A general algorithms for short-term scheduling of batch operations. 2. Computational issues. Comput. Chem. Eng. 1993, 17, 229-244. (40) Syarif, A.; Yun, Y.; Gen, M. Study on multi-stage logistic chain network: A spanning tree-based genetic algorithm approach. Comput. Ind. Eng. 2002, 43, 299-314. (41) Vergara, F. E.; Khouja, M.; Michalewicz, Z. An evolutionary algorithm for optimizing material flow in supply chains. Comput. Ind. Eng. 2002, 43, 407-421.

Ind. Eng. Chem. Res., Vol. 45, No. 22, 2006 7581 (42) Zhou, G.; Min, H.; Gen, M. The balanced allocation of customers to multiple distribution centers in the supply chain network: A genetic algorithm approach. Comput. Ind. Eng. 2002, 43, 251-261. (43) Wang, C.; Quang, H.; Xu, X. Optimal design of multiproduct batch chemical process using genetics algorithms. Ind. Eng. Chem. Res. 1996, 35, 3560-3566. (44) Tan, S.; Mah, R. S. H. Evolutionary design of noncontinuous plants. Comput. Ind. Eng. 1998, 22 (1-2), 69-85. (45) Reklaitis, G. V. Review of scheduling of process operations. AIChE Symp. Ser. 1982, 78 (214), 119-133. (46) Balasubramanian, J.; Grossmann, I. E. Approximation to multistage stochastic optimization in multiperiod batch plant scheduling under demand uncertainty. Ind. Eng. Chem. Res. 2004, 43, 3695-3713. (47) Works, T. M. MATLAB 7.0 (R14), Simulink 6.0, Stateflow 6.0, User’s Manual; The Math Works Inc.: Natick, MA, 2004.

(48) Brooke, A.; Kendrik, D.; Meeraus, A.; Raman, R.; Rosenthal, R. E. GAMSsA User’s Guide; GAMS Development Corporation: Washington, DC, 1998. (49) Ferris, M. C. MATLAB and GAMS: Interfacing optimization and Visualization software. Technical Report; Computer Sciences Department, University of Wisconsin: Madison, WI, 1998. (50) Maravelias, C. T.; Grossmann, I. E. A hybrid MILP/CP decomposition approach for the continuous time scheduling of multipurpose batch plants. Comput. Chem. Eng. 2004, 28, 1921-1949.

ReceiVed for reView December 21, 2005 ReVised manuscript receiVed June 30, 2006 Accepted August 24, 2006 IE051424A

Addressing the Design of Chemical Supply Chains under Demand

Recommend Documents