• Editorial Board +
• For Contributors +
• Journal Search +
Journal Search Engine
ISSN : 1598-7248 (Print)
ISSN : 2234-6473 (Online)
Industrial Engineering & Management Systems Vol.17 No.3 pp.417-433
DOI : https://doi.org/10.7232/iems.2018.17.3.417

# A Multi-Stage Stochastic Mixed-Integer Linear Programming to Design an Integrated Production-Distribution Network under Stochastic Demands

Mohammad Derakhshi, Seyed Taghi Akhavan Niaki*, Seyed Armin Akhavan Niaki
Department of Industrial Engineering, Sharif University of Technology, Tehran, Iran
Department of Statistics, Eberly College of Arts and Sciences, West Virginia University, Morgantown
Corresponding Author, E-mail: Niaki@Sharif.edu
September 13, 2017 May 21, 2018 June 7, 2018

## ABSTRACT

Supply chain management has gained much interest from researchers and practitioners in recent years. Proposing practical models that efficiently address different aspects of the supply chain is a difficult challenge. This research investigates an integrated production-distribution supply chain problem. The developed model incorporates parties with a specified number of processes to obtain raw materials from the suppliers in order to convert them to semi and final products. These products are then distributed through warehouses to end-distributors having uncertain demands. This uncertainty is captured as a dynamic stochastic data process during the planning horizon and is modeled into a multi-stage stochastic mixed integer linear program using a scenario tree approach. For large-size instances, a hybrid exact-approximate algorithm is proposed, where its effectiveness is assessed via several numerical cases. Furthermore, the model is generalized to its bi-objective version by considering the accessibility of the products based on the safety stock policy of the companies involved. In the end, an existing algorithm is combined with the ε-constraint method to obtain an approximate Pareto front.

## 1. INTRODUCTION

An integrated process involving suppliers, manufacturers, warehouses, distributors, and retailers that acquire raw materials, converts them into specified final products, and delivers them to the retailers is called a supply chain (SC). The flows in this chain are a forward flow of materials and the backward flow of information (Beamon, 1998). Altering aspects such as increasing the productivity, quality, and accessibility of products, innovations, etc. make it more difficult to determine the optimal behavior of the flow for such companies that fierce competition forces them to improve the efficiency of their manufacturing operations. In this research, mathematical programs are used as a reliable tool to obtain a solution. Note that providing effective models that capture different aspects of a given SC, has always been of interest to researchers.

Production and distribution planning are the two major optimization problems in any given SC. In production planning problems, decisions pertaining to hiring and firing labor, setting regular and overtime production, subcontracting, and machine capacity levels are made for a definite planning horizon. In distribution planning problems, the aim is to determine which facility has the ability to satisfy the demand of the intended markets (Fahimnia et al., 2013). Production-distribution (P-D) networks provide an effective tool to model manufacturing and logistics activities of a company (Dasci and Verter, 2001).

In this paper, a multi-stage stochastic mixed integer linear model is proposed for a complex P-D problem under uncertain demand, where the demand is presented as a scenario tree which is taken out of the scenario reduction method. Then, a hybrid exact-approximate algorithm is used to obtain the solution. Next, the single-objective optimization model is converted to a bi-objective version by considering the safety stock policy; a case that has not been addressed in the literature. Another solution algorithm namely, the epsilon-constraint method, is utilized to solve the proposed bi-objective optimization problem.

The rest of the paper is organized as follows: relevant literature is reviewed in Section 2. Section 3 provides a detailed explanation of the problem. The mathematical formulation is presented and analyzed in Section 4. Section 5 demonstrates the solution method. Section 6 presents numerical results and Section 7 encompasses the conclusion and discusses possible grounds for future research.

## 2. LITERATURE REVIEW

A significant number of research studies have addressed the P-D planning problem. Comprehensive literature reviews of previous studies on P-D planning problem are given in Beamon (1998), Fahimnia et al. (2013). In this section, a number of recent and relevant works along with their main characteristics are reviewed and presented. The main focus of the review is on stochastic and multiobjective models, while few studies relating to deterministic and robust models are reviewed.

The model suggested in Mohamed (1999) is an integrated production and distribution model for multinational companies operating under different inflation and exchange rates. Various government regulations and many international uncertainties complicate decision making for companies operating in these situations. A linear programming model was developed by Kanyalkar and Adil (2005) for an integrated multi-item production, parallel multi-plant production, and dynamic distribution. They also demonstrated the capability of the proposed model in coordinating the production and distribution planning of a company with multi-location parallel manufacturing facilities in a consumer-goods company. Boudia et al. (2007) presented a model to minimize the sum of setup, inventory and transportation network costs to determine production and routing decisions of the network. In their work, a meta-heuristics algorithm was employed as the solution method and was improved on through a reactive mechanism alongside a path-relinking process. A global and multi- national production-distribution supply chain planning problem, which is affected by government regulations was presented by Das and Sengupta (2009). Uncertainty in demand and transportation time along with multinational regulation parameters such as exchange rate and border crossing costs are considered in this study. In addition, a Tabu search algorithm was proposed by Bard and Nananukul (2009) for solving a simple production and distribution planning problem.

Nishi et al. (2007) investigated an integrated production scheduling and distribution planning optimization problem for an aluminum rolling processing line. The problem was formulated as a mixed integer linear program; an augmented Lagrangian approach was used to decompose the problem into two subproblems: scheduling and warehousing sub-problems. The scheduling subproblem was solved using a simulated annealing method and the warehousing sub-problem was solved using the commercial mixed integer linear programming (MILP) solver of the CPLEX software. Note that this study did not consider end-users and demand fluctuations. A just-intime distribution supply chain network including several products and a single plant was considered by Farahani and Elahipanah (2008). They proposed a bi-objective mixed-integer linear model to minimize both the network costs and the sum of back-orders and surpluses of products in all periods, simultaneously. A hybrid non-dominated sorting genetic algorithm was utilized in their work to solve real-size instances. Hamedi et al. (2009) introduced a six-level natural gas supply chain problem. A mixed integer nonlinear program based on the concepts involved in a gas company and the characteristics of the distribution network was developed for this research. A hierarchical algorithm based on the relations involved in a natural gas network was developed to attain a near optimal solution in a reasonable time. A multi-agent system of a P-D supply chain was studied by Kazemi et al. (2009), where three genetic algorithms with different crossover operators were developed to solve the problem. They compared the solutions obtained by their algorithms with those obtained by a Lagrangian algorithm, based on which the efficiency of their proposed method was proven. A GA-based solution procedure was developed by Yimer and Demirli (2010) to solve the scheduling problem of a built-to-order supply chain. Moreover, a two-echelon SC problem based on the integration of the aggregate production and distribution plan was studied by Fahimnia et al. (2012). They developed a mixed integer program to formulate the problem and utilized a GA as the solution method. Although their proposed algorithm provided good performance in solving medium-size instances, it failed in solving large-size problems.

In a study by Kanyalkar and Adil (2010), a robust optimization model for an integrated multi-site problem that incorporated both model and solution robustness was proposed. In their work, product demands were considered stochastic and were encompassed in diverse scenarios. The effectiveness of their proposed methodology was validated using several numerical examples. Safaei et al. (2010) used a hybrid mathematical and simulation technique to solve an integrated multi-product, multiperiod, multi-site production and distribution planning problem. However, their mixed integer linear program model ignored many characteristics of the real-world (e.g. backlogging costs, production costs, and detail of production). Nasiri et al. (2014) proposed a hierarchical model to design an integrated production-distribution network of a three-echelon supply chain with stochastic demands. They incorporated a capacitated facility location model in order to determine the locations of multi-level warehouses and employed a Lagrangian relaxation approach for solving the problem. Yu et al. (2015) developed a model to integrate the location, production, and distribution problem of a multi-product, multi-echelon, multi-period supply chain. The applicability of their model was confirmed using small and medium size instances solved by LINGO. They concluded that the customer demand parameter had a major impact on the optimal solution in their study.

Liu and Papageorgiou (2013) considered responsiveness and the customer service level of a global supply chain. They proposed a multi-objective mixed-integer linear program and used both the ε-constraint and the lexicographic minimax methods to solve the problem. Yoshida et al. (2014) examined the use of quantity discounts in a multiperiod production planning with a single supplier and single retailer under uncertain demands in a global supply chain environment. The problem was formulated as a stochastic multi-period bi-level program. A multi-objective bi-level P-D planning model that incorporated equilibrium between supply and demand was investigated by Jia et al. (2014). They employed a hybrid genetic algorithm to solve the problem. Moreover, Sel et al. (2015) developed a mixed integer linear program and a decomposition heuristic approach for a dairy industry supply chain. Their proposed model conducted both tactical and operational decisions properly.

A multiperiod, multi-product production planning problem under uncertainty in demand and quality of raw materials using a hybrid scenario tree approach was formulated in Kazemi Zanjani et al. (2010). They proposed a multistage stochastic model to address the problem and employed it in a real-world sawmill. In addition, a scenario- based model of an aggregate P-D planning problem was studied by Niknamfar et al. (2015). Their model incorporated multiple production facilities, several distribution centers, and multiple customer zones with uncertain parameters. Most of the cost parameters such as sale prices, fixed costs, production costs, storage costs, and transportation costs were assumed to be uncertain fuzzy parameters. Moreover, Sarrafha et al. (2015) developed a biobjective mixed-integer nonlinear programming (MINLP) model to design a multiperiodic structure of a supply chain network involving suppliers, factories, distribution centers (DCs), and retailers. Minimizing the total costs as well as minimizing the average tardiness of the product to DCs were the two objectives of their work. They utilized a multi-objective biogeography-based optimization (MOBBO) algorithm with tuned parameters in order to obtain a near-optimum solution. Finally, using mixed-integer nonlinear programming, Esmaeilikia et al. (2016) proposed a tactical SC planning model that considered a range of flexibilities in sourcing, production, and logistics, to cope with uncertainty.

This paper differs from the published papers in the literature in the following three main aspects

• 1) It incorporates the demand uncertainty into the model using the scenario tree approach, while almost all of the published papers addressed this issue through the robust optimization formulations that are naturally based on the worst-based scenario

• 2) It provides a stochastic mixed integer linear program to the problem, while most of the published papers proposed non-linear mixed integer programs for which a general effective algorithm was not proposed for a solution yet

• 3) It tackles the problem using a meta-heuristic algorithm and evaluates the quality of the solutions using a few performance measures that have not been applied previously.

## 3. PROBLEM DESCRIPTION

This research investigates a P-D network in which a set of manufacturers denoted by ${ 1 , … , ​ M − 1 , M }$ cooperates together to produce i types of products by assembling p unique components during t time periods. The set of manufacturers is divided into two subsets: the sub-manufacturers incorporating the set of ${ 1 , … , ​ M − 1 }$, and the main manufacturer denoted by {M}. There are s raw-material suppliers who provide the primary requirements of k types of raw-materials of both subs- and main manufacturers. Each of the components is produced solely in a unique sub-manufacturer. Furthermore, each of the sub-manufacturers is capable of producing at least one component with respect to their production technology and constraints. All the manufacturing operations are performed in regular as well as overtime periods. Moreover, the main manufacturer outsources a specific component, if none of the sub-manufacturers is able to produce it.

In order to meet the demand, the main manufacturer delivers finished products to e end-distributors through w warehouses. At the beginning of the planning horizon, in order to expedite the demand of the end distributors, a number of warehouses are selected. Each of the end distributors in each time period may only be served by exactly one warehouse. All the manufacturers have a minimum production level to meet in order to keep continuity in their production lines.

A schematic overview of the proposed P-D problem is illustrated in Figure 1. Binary decision variables at the beginning of the planning horizon determine the selected warehouses and their corresponding end-distributors. The rest of the decision variables are as follows:

• - The total units of raw materials to be transmitted from the supplier to each plant,

• - The total units of components to be produced at each sub-manufacturer in regular and overtime periods,

• - The total units to be outsourced by the main manufacturer,

• - The total units of components to be conveyed from sub-manufacturers to the main manufacturer,

• - The total number of items to be produced by the main manufacturer in regular time and overtime,

• - The total units of products to be conveyed from the plants to warehouses,

• - The total units of products to be transported from warehouses to the end distributors,

• - The inventory level of raw materials held at each plant,

• - The components inventory held at each plant, and

• - The inventory level of finished products at the warehouse and end-distributors.

As it can be seen in Figure 1, the location of the warehouses and their flow between end-distributors are characterized by dashed lines to show their possible selections. In other words, after determining the warehouses to be constructed, each of the end-distributors is supported by one of them.

## 4. METHODOLOGY

In this section, the uncertainty in demand is first modeled using the scenario tree approach. Then, a multistage stochastic mixed integer linear program is developed based on this model. Finally, the safety stock policy is added up to the model in order to extend it to its biobjective version.

### 4.1. The Scenario Tree Approach

The quality/suitability evaluation of scenario-generation methods was discussed by Kaut and Wallace (2003), where they reviewed the most common scenario generation methods. A common method for scenario tree generation is based on scenario reduction proposed by Heitsch and Römisch (2003). In another study, the stability of the multistage stochastic program in the scenario reduction method was proved by Heitsch and Römisch (2010). This method is utilized in this paper to generate different demand scenarios.

The scenario tree is an effective way to describe the uncertainty involved in mathematical programs. With the tree representation, decision vectors that are non-anticipative, i.e. decisions taken at any stage that is dependent on a past realization of the random data process and cannot depend on the future realization, can be taken on by any node of the tree. The current state is invoked by the root node of the tree. Next stage scenarios are represented by branches of the tree where a probability is assigned to express the probability of the corresponding state. Each node has a single ancestor and possibly several descendants. The set of possible future sequences of data outcomes is called scenario (Kazemi Zanjani et al., 2010).

While it is assumed that products are shipped from warehouses to end-distributors prior to observing the actual demand, the end distributors’ demands for products are considered random variables that follow a lognormal distribution with mean μi and variance $σ ε 2$ , with the assumption that all demand realizations are within the interval $( μ i − 3 σ e , μ i + 3 σ e ) .$. This interval is then partitioned to sub-intervals $( μ i − 3 σ e , μ i − 2 σ e ) , ( μ i − 2 σ e , μ i − σ e ) ,$ $( μ i − σ e , μ i + σ e ) , ( μ i + σ e , μ i + 2 σ e ) , ( μ i + 2 σ e , μ i + 3 σ e )$ which correspond to very low, low, moderate, high, and very high states, respectively. In addition, it is assumed that at least one of these states (events) is realized at any stage. Each node corresponds to one of these states and contains i×e data generated from the underlying distribution function with its related interval. A set of scenarios generated by this method is called a scenario fan (Eichhorn and Romisch, 2006). The number of scenarios grows exponentially by increasing the number of planning periods in the scenario fan. Thus, the P-D supply chain problem at hand would be intractable. However, the scenario reduction method that will be further implemented in GAMS (the general algebraic modeling system) software using the Scenred module reduces this inherent intractability (Heitsch and Römisch, 2003).

### 4.2. The Stochastic Programming Model

The absence of accurate information in real-world problems results in demand uncertainty. This poses a significant challenge for designing an appropriate production and inventory control model. The following notations are used to develop a multistage stochastic program for the proposed P-D problem under investigation. Note that the demand for products in the end distributors is assumed to be an uncertain parameter.

#### 4.2.1. Indices and Sets

The indices and sets are defined as:

• n Index of scenario tree nodes

• M Map between nodes and time periods ($ℳ = { ( t , n ) | t = 1 , … , T }$)

• f(n)  Index for the ancestor of node n in scenario tree

• p Index of components

• i Index of products

• m Index of submanufacturers

• M Index of the main manufacturer

• w Index of warehouses

• c Index of end-distributors

• s Index of raw material suppliers

• k Index of raw materials

• t Index of time periods

• A Set of components which can be produced in sub-manufacturers’ plants

• B Set of components which outsourced by the main manufacturer

• Cm Set of raw materials which are used in manufacturer m plant

• Ds Set of raw materials provided by supplier s

• $ℰ$ Set of suppliers that can supply raw material k

• $ℱ i$ Set of raw materials to be used in the assembly of product i

• $𝒢 𝓟$ Set of raw materials to be used in the production of component p

• $ℋ m$ Set of components that can be produced by sub-manufacturer m

#### 4.2.2. Parameters

The parameters involved in the mathematical formulation of the problem are:

• Dietn Demand for product i by end-distributor e during time period t and node n

• prn Occurrence probability of node n

• Fw Fixed costs of opening and operating warehouse w

• $H p t w$ Unit holding cost for component p during period t in sub-manufacturers’ plants

• $H p t w M$ Unit holding cost of component p during period t in the main manufacturer’s plant

• Hiwt Unit holding cost for product i in warehouse w during period t

• $H i e t '$ Unit holding cost for product i in enddistributor e during period t

• $H k m t R$ Unit holding cost for raw material k in manufacturer m plant during period t

• TRksmt Unit transportation cost for raw material k from supplier s to manufacturer m during period t

• Tiwt Unit transportation cost of product i from the main manufacturer to warehouse w during period t

• $T i w e t '$ Unit transportation cost of product i from warehouse w to end-distributor e during period t

• TVpt Unit transportation cost for component p from its unique sub-manufacturer to the main manufacturer during period t

• api Quantity of component p required per unit of product i

• τkp Amount of raw material k used to produce a unit of component p

• $τ k t '$ Unit cost of raw material k provided by raw material suppliers during period t

• $τ k i M$ Amount of raw material k used to produce a unit of product i in the main manufacturer’s plant

• $υ p$ Processing time to produce a unit of component p

• $σ p t$ Labor/hour cost for regular-time production of component p during period t

• $σ p t '$ Labor/hour cost for overtime production of component p during period t

• $υ i M$ Processing time to produce a unit of product i

• $σ i t M$ Labor/hour cost for regular-time production of product i during period t

• $σ i t ' M$ Labor/hour cost for overtime production of product i during period t

• Gpt Unit outsourcing cost of component p during period t

• Lit Minimum production quantity of product i during period t

• Kpt Minimum production quantity of component p during period t

• Qwt Holding capacity of warehouse w during period t

• $Q e t '$ Holding capacity of end-distributor e during period t

• ηi Empty space assigned to a unit of product i

• $η p '$ Space occupied by a unit of component p

• ρk Space occupied by a unit of raw material k

• $W m t M a x$ Holding capacity for produced components in sub-manufacturer m in period t

• $W t M a x M$ Holding capacity of products’ components in the main manufacturer plant in period t

• $R m t M a x$ Holding capacity of raw materials at manufacturer m plant during period t

• $C P m t M a x$ Transportation capacity of components produced by sub-manufacturer m during period t

• $Q R s t M a x$ Capacity of raw materials supplied by supplier s during period t

• $λ m t$ Hours capacity for regular-time production at manufacturer m during period t

• $λ m t '$ Hours capacity for overtime production at manufacturer m during period t

• $λ t M$ Hours capacity for regular-time production at the main manufacturer during period t

• $λ t ' M$ Hours capacity for overtime production at the main manufacturer during period t

• $γ p 1$ Inventory level of component p on its own sub-manufacturer plant at the beginning of the planning horizon(t = 1)

• $γ p 1 M$ Inventory level of component p in the main manufacturer's plant at the beginning of the planning horizon(t = 1)

• $θ i w 1$ Inventory level of product i in warehouse w at the beginning of the planning horizon (t = 1)

• $ε i e 1$ Inventory level of product i in end-distributor e at the beginning of the planning horizon (t = 1)

• $κ k m 1$ Inventory level of raw material k in manufacturer m at the beginning of the planning horizon(t = 1)

#### 4.2.3. Decision Variables

The decision variables are:

• Iitn The quantity of product i produced in regular time at the main manufacturer’s plant during period t at node n

• $I i t n '$ The quantity of product i produced in overtime at the main manufacturer’s plant during period t at node n

• Uptn The quantity of component p produced in regular-time during period t at node n

• $U p t n '$ The quantity of component p produced in overtime during period t at node n

• $U p t n ' '$ The quantity of component p outsourced by the main manufacturer during period t at node n

• Wptn The Inventory level of component p at the beginning of period t at node n

• $W p t n M$ The Inventory level of component p at the main manufacturer at the beginning of period t at node n

• $R k m t n$ The Inventory level of raw material k at manufacturer m at the beginning of period t at node n

• $V i w t n$ The Inventory level of product i at warehouse w at the beginning of period t at node n

• $V i e t n '$ The Inventory level of product i at enddistributor e at the beginning of period t at node n

• Xksmtn The amount of raw material k transmitted from supplier s to manufacturer m during period t at node n

• Jiwtn The quantity of product i transmitted from the main manufacturer to warehouse w during period t at node n

• $J i w e t n '$ The amount of product i transmitted from warehouse w to end-distributor e during period t at node n

• Yptn The amount of component p transmitted from its unique manufacturer during period t at node n

#### 4.2.4. The Multi-Stage Stochastic Programming Model

Using the indices, parameters and decision variables defined in Section 4.2.3, the multi-stage stochastic programming model of the problem is derived as:(5)(6)(29)(30)(31)(33)(34)(35)(36)(37)(38)(39)(40)(41)(42)

$M i n Z = ∑ w F w O w + ∑ ( t , n ) ∈ ℳ p r n [ ∑ m ∑ s ∑ k ∈ m ∩ s T R k s m t X k s m t n + ∑ i ( ∑ k ∈ ℱ i τ k i M τ k t ' + ν i M σ i t M ) I i t n + ∑ i ( ∑ k ∈ ℱ i τ k i M τ k t ' + ν i M σ i t ' M ) I i t n ' + ∑ p ∈ ( ∑ k ∈ p τ k p τ k t ' + υ p σ p t ) U p t n + ∑ p ∈ ( ∑ k ∈ p τ k p τ k t ' + υ p σ p t n ' ) U p t n ' + ∑ p ∈ ℬ G p t U p t n ' ' + ∑ p ∈ H p t w W p t n + ∑ p H p t w M W p t n M + ∑ p ∈ T V p t Y p t n + ∑ i ∑ w H i w t V i w t n + ∑ i ∑ e H i e t ' V i e t n ' + ∑ m ∑ k ∈ m H k m t R R k m t n + ∑ i ∑ w T i w t J i w t n + ∑ i ∑ w ∑ e T i w e t ' J i w e t n ' ]$
(1)

s.t.

$R k m t n = R k m ( t − 1 ) f ( n ) + ∑ s ∈ ℰ k X k s m ( t − 1 ) f ( n ) − ∑ p ∈ ℋ m τ k p ( U p ( t − 1 ) f ( n ) + U p ( t − 1 ) f ( n ) ' ) ∀ m ≠ M , k ∈ m , ( t , n ) ∈ ℳ ' '$
(2)

$R k M t n = R k M ( t − 1 ) f ( n ) + ∑ s ∈ ℰ k X k s M ( t − 1 ) f ( n ) − ∑ i τ k i M ( I i ( t − 1 ) f ( n ) + I i ( t − 1 ) f ( n ) ' ) ∀ k ∈ M , ( t , n ) ∈ ℳ ' '$
(3)

$∑ p ∈ ℋ m υ p U p t n ≤ λ m t ∀ m ≠ M , ( t , n ) ∈ ℳ '$
(4)

$∑ p ∈ ℋ m υ p U p t n ' ≤ λ m t ' ∀ m ≠ M , ( t , n ) ∈ ℳ '$
(5)

$∑ i υ i M I i t n ≤ λ t M ∀ ( t , n ) ∈ ℳ '$
(6)

$∑ i υ i M I i t n ' ≤ λ t ' M ∀ ( t , n ) ∈ ℳ '$
(7)

$∑ i J i w t n ≤ Q w t O w ∀ w , ( t , n ) ∈ ℳ '$
(8)

$∑ i J i w e t n ' ≤ Q w t O w e ' ∀ w , e , ( t , n ) ∈ ℳ '$
(9)

$O w e ' ≤ O w ∀ w , e$
(10)

$∑ w O w e ' = 1 ∀ e$
(11)

$∑ e O w e ' ≥ O w ∀ w$
(12)

$I i t n + I i t n ' ≥ L i t ∀ i , ( t , n ) ∈ ℳ '$
(13)

$∑ w J i w t n = I i t n + I i t n ' ∀ i , ( t , n ) ∈ ℳ '$
(14)

$V i w t n = V i w ( t − 1 ) f ( n ) + J i w ( t − 1 ) f ( n ) − ∑ e J i w e ( t − 1 ) f ( n ) ' ∀ i , w , ( t , n ) ∈ ℳ ' '$
(15)

$V i e t n ' = V i e ( t − 1 ) f ( n ) ' + ∑ w J i w e ( t − 1 ) f ( n ) ' − D i e t n ∀ i , e , ( t , n ) ∈ ℳ ' '$
(16)

$W p t n = W p ( t − 1 ) f ( n ) + U p ( t − 1 ) f ( n ) + U p ( t − 1 ) f ( n ) ' − Y p ( t − 1 ) f ( n ) ∀ p ∈ , ( t , n ) ∈ ℳ ' '$
(17)

$W p t n M = W p ( t − 1 ) f ( n ) M + Y p ( t − 1 ) f ( n ) − ∑ i a p i ( I i ( t − 1 ) f ( n ) + I i ( t − 1 ) f ( n ) ' ) ∀ p ∈ , ( t , n ) ∈ ℳ ' '$
(18)

$W p t n M = W p ( t − 1 ) f ( n ) M + U p ( t − 1 ) f ( n ) ' ' − ∑ i a p i ( I i ( t − 1 ) f ( n ) + I i ( t − 1 ) f ( n ) ' ) ∀ p ∈ ℬ , ( t , n ) ∈ ℳ ' '$
(19)

$U p t n + U p t n ' ≥ K p t ∀ p ∈ , ( t , n ) ∈ ℳ '$
(20)

$∑ k ∈ m ρ k R k m t ≤ R m t M a x ∀ m , ( t , n ) ∈ ℳ ' '$
(21)

$∑ m ∑ k ∈ m ∩ s X k s m t n ≤ Q R s t M a x ∀ s , ( t , n ) ∈ ℳ '$
(22)

$∑ p ∈ ℋ m η p ' W p t ≤ W m t M a x ∀ m ≠ M , ( t , n ) ∈ ℳ ' '$
(23)

$∑ p η p ' W p t n M ≤ W t M a x M ∀ ( t , n ) ∈ ℳ ' '$
(24)

$∑ p ∈ ℋ m Y p t n ≤ C P m t M a x ∀ m , ( t , n ) ∈ ℳ '$
(25)

$∑ i η i V i w t n ≤ Q w t O w ∀ w , ( t , n ) ∈ ℳ ' '$
(26)

$∑ i η i V i e t ' ≤ Q e t ' ∀ e , ( t , n ) ∈ ℳ ' '$
(27)

$R k m t n = κ k m t ∀ m , k ∈ m & ( t , n ) = ( t 1 , n 1 )$
(28)

$W p t n = γ p t ∀ p ∈ & ( t , n ) = ( t 1 , n 1 )$
(29)

$W p t n M = γ p t M ∀ p & ( t , n ) = ( t 1 , n 1 )$
(30)

$V i w t n = θ i w t ∀ i , w & ( t , n ) = ( t 1 , n 1 )$
(31)

$V i e t n ' = ε i e t ∀ i , e & ( t , n ) = ( t 1 , n 1 )$
(32)

$I i t n , I i t n ' ≥ 0 ∀ i , ( t , n ) ∈ ℳ '$
(33)

$U p t n , U p t n ' , Y p t n ≥ 0 ∀ p ∈ , ( t , n ) ∈ ℳ '$
(34)

$W p t n ≥ 0 ∀ p ∈ , ( t , n ) ∈ ℳ ' '$
(35)

$U p t n ' ' ≥ 0 ∀ p ∈ ℬ , ( t , n ) ∈ ℳ '$
(36)

$W p t n M ≥ 0 ∀ p , ( t , n ) ∈ ℳ ' '$
(37)

$V i w t n , V i e t n ' ≥ 0 ∀ i , w , e , ( t , n ) ∈ ℳ ' '$
(38)

$J i w t n , J i w e t n ' ≥ 0 ∀ i , w , e , ( t , n ) ∈ ℳ '$
(39)

$X k s m t n ≥ 0 ∀ s , m , k ∈ m ∩ s , ( t , n ) ∈ ℳ '$
(40)

$R k m t n ≥ 0 ∀ m , k ∈ m , ( t , n ) ∈ ℳ ' '$
(41)

$O w , O w e ' ∈ { 0 , 1 } ∀ w , e$
(42)

such that

$ℳ ' ( ⊂ ℳ ) = { ( t , n ) | t = 1 , … , T − 1 } ℳ ' ' ( ⊂ ℳ ) = { ( t , n ) | t = 2 , … , T }$

The objective function in (1) aims to minimize all associated costs of the investigated P-D supply chain problem, i.e. total production, outsourcing, inventory holding, and transportation costs. Constraint sets (2) and (3) represent the inventory balance for raw materials in both sub and main manufacturers. Manufacturing capacity constraints for regular and overtime production in the sub and main manufacturers are captured through constraint sets (4)-(7), respectively. Constraint (8) enforces that products are conveyed to warehouses only when their production is complete. Constraint (9) ensures that products are transported from a specific warehouse to enddistributors within the coverage range. Constraint (10) restricts warehouses to serve end-distributors when they are constructed. Constraint (11) indicates that each one of the end distributors is served by exactly one of the candidate warehouses. Constraint (12) ensures that every selected warehouse covers at least one end-distributor. Constraint (13) imposes the minimum production level in the main manufacturer’s plant. Constraint set (14) ensures that all products are transported directly to the warehouse. Constraint sets (15) and (16) ensure the inventory balance of products in warehouses and end-distributors, respectively. The balance inventory of components in the sub and main manufacturers is guaranteed by constraint sets (17) and (18), respectively. Constraint set (19) represents the balanced inventory of outsourced components in the main manufacturer. Constraint (20) ensures a minimum level of production in the sub-manufacturers. Constraint set (21) guarantees the holding capacity of raw materials in the manufacturers’ plants. Constraint set (22) is for the supplier’s capacity. Constraint sets (23) and (24) are the components holding capacity at manufacturers. The transportation capacity for components from sub-manufacturer to the main manufacturer is controlled by constraint set (25). Constraint sets (26) and (27) are the products holding capacity in warehouses and end-distributors, respectively. Constraints (28)-(32) limit the inventory level of raw material, components, and finished products at sub and main-manufacturers, warehouses, and end-distributors, respectively. Finally, the remaining constraints define the types of the variables used.

### 4.3. The Bi-Objective Model

In real-world situations, accurate predictions of uncertain parameters are rarely available (Chopra and Meindl, 2007). Thus, the safety stock policy is used to mitigate the effect of a sudden increase in demand. This issue is especially important in industries where the risk of product obsolescence is high. To incorporate this concept in the model developed in Section 0, an attempt is made in this paper to impose a safety stock to the inventory of end-distributors as much as possible. This guarantees high accessibility for customers and subsequently a higher customer service level, but increases the supply chain costs. Hence, the following objective function can be used to address this issue:(43)

$M i n Z 2 = ∑ ( t , n ) ∈ ℳ p r n ∑ i ∑ e ( | V i e t n ' − S S i e | ) ,$
(43)

where $| . |$ is the absolute value function. As seen, the above function is nonlinear but can be linearized using the following:(44)

$V i e t n ' − S S i e : = E X i e t n − S H i e t n .$
(44)

Note that the parameter SSie is the safety stock of product i at end-distributor e which is mandated by the main decision maker to the model developer. Based on the above definition, we are minimizing:(45)(46)

$Z 2 = ∑ ( t , n ) ∈ ℳ p r n ∑ i ∑ e ( α 1 E X i e t n + α 2 S H i e t n )$
(45)

$s . t . : E X i e t n , S H i e t n ≥ 0 ∀ i , e , ( t , n ) ∈ ℳ ' '$
(46)

where α1, α2 are the weights of EXietn and SHietn, respectively.

## 5. SOLUTION METHOD

In this section, a hybrid exact-approximate algorithm is proposed to solve the multistage stochastic programming problem modeled in Section 4.2.4. In addition, an algorithm proposed by Ehrgott (2006) is combined with the ε-constraint method to obtain approximate Pareto front of the bi-objective optimization problem modeled in Section 4.3.

### 5.1. The Hybrid Exact-Approximate Algorithm

The “matout.gms” and the “gdxmrw” modules in GAMS were introduced by Ferris (1998) and Ferris et al. (2011), respectively in order to be linked with MATLAB in solving multistage stochastic programming optimization problems. The key idea behind the solution procedure proposed in this research is to satisfy a certain constraint while efficiently handling the remaining constraints. The proposed solution procedure imputes binary decision variables as a cover. Then, in order to obtain an efficient linear model which can be solved with GAMS, a cover is randomly generated in MATLAB such that it satisfies constraints 10-12. In order to avoid infeasible solutions, an infeasible cover associated with a big penalty that is added to the objective function is utilized. The two general steps involved in the proposed solution procedure explained above are illustrated in Algorithms 1) and 2, respectively. In the end, a GA-based algorithm (Algorithm 3) is developed to obtain an upper bound solution of the problem.

### 5.2. The Ehrgott (2006) Algorithm

The general form of a bi-objective optimization problem is shown in (47)

(47)

where $Z k , k = 1 , 2$ are the objective functions and X is the set of feasible solutions. As one is often faced with the objectives that are incompatible, there is no solution that optimizes both simultaneously. Therefore, the Pareto optimality concept replaces the former definition of the optimum solution (Ehrgott, 2006).

Definition 5.1: A feasible solution, $x ^ ∈ ,$, is called the Pareto optimal solution if there is no $x ∈$ such that $f k ( x ) ≤ f k ( x ^ )$ for k = 1, 2 and at least for one $j ∈ { 1 , 2 } w e h a v e f j ( x ) < f j ( x ^ )$. If $x ^$is the Pareto solution, then $f ( x ^ )$ is called non-dominated point. The set of all Pareto solutions is denoted by $E$ and is called efficient set. The set of all nondominated points is denoted by YN and called the Pareto front or the nondominated set.

Definition 5.2: The point $y I = ( y 1 I , y 2 I )$ given by(48)

(48)

is called the Ideal point of the bi-objective optimization problem. Moreover, the point $y N = ( y 1 N , y 2 N )$ given by(49)

(49)

is called the Nadir point of the bi-objective optimization problem. Ehrgott (2006) presented an algorithm to obtain ideal and nadir points of the above bi-objective optimization problems.

#### 5.2.1. The ε-Constraint Method

The ε-constraint method, firstly introduced by Vira and Haimes (1983), is probably the best-known technique for solving multi-objective optimization problems. In this method, one of the original objectives is minimized, while the others are transformed to constraints (Ehrgott, 2006). The bi-objective optimization problem is substituted by the following ε-constraint problem:(50)

(50)

where $i , j ∈ { 1 , 2 }$ and ij.

Algorithm 4 obtains an approximate Pareto front for the bi-objective optimization problem (47) by solving ( ) i j P ε problems iteratively. In this algorithm, $Δ = Z 2 ∗ N − Z 2 ∗ I n$, where n is an appropriate natural number provided by the main decision maker.

## 6. COMPUTATIONAL RESULTS

The effectiveness of the multistage stochastic program and the bi-objective formulation, as well as the applications of both the solution methods (Algorithm 3 and 4), is demonstrated in this section via several numerical examples.

### 6.1. The Efficiency of the Exact-Approximate Solution Method

In order to demonstrate the efficacy of the proposed exact-approximate solution method, various problem instances are randomly generated with different sizes. Table 1 shows the size of each test problem along with their parameters. In addition, Tables A1 and A2 in Appendix A show the randomly generated parameters of each instance. Moreover, the demand quantities in the reduced scenario tree of the first test problem are shown in Table B1 (Appendix B). As it can be seen, the main manufacturer produces two products by assembling three components. The first component is outsourced while the second and the third component are produced by the first and the second sub-manufacturer, respectively.

The performance of the proposed algorithm is compared with the mixed-integer solver of GAMS in terms of percentage deviation defined as (Farahani and Elahipanah, 2008):(51)

$C . R . = Z G A − Z G A M S Z G A M S ,$
(51)

where ZGA is the objective function value found by Algorithm 3 and ZGAMS is the objective function value obtained from the GAMS software. All the test problems are formulated and solved by GAMS 24.1.2 on a laptop with an Intel (R) Core (TM) i7-4500U CPU @ 1.80GHz processor, with 6GB RAM. It is worth noting that in cases for which the GAMS solver could not reach an optimal solution due to the lack of memory size, the best incumbent solution is considered. Moreover, the detailed solution of the first test problem is summarized in Table 2. As it can be seen in Table 3 in large instances for which no optimal solution can be found by the GAMS solver, Algorithm 3 finds a near-optimum solution in a reasonable computational time. Besides, while the C.R. values for the small and medium size instances is either zero or very close to zero, Algorithm 3 performs significantly better than the solver in GAMS in terms of the elapsed computational time. Moreover, Figure 2 shows the initial and the reduced demand scenario trees.

### 6.2. The Efficiency of the Hybrid Ehrgott-ε- Constraint Algorithm

Figure 3 shows the Pareto solutions of the first problem (as a representative of all test problems that can be solved in a reasonable computational time, i.e. small and medium size instances) found by both the hybrid Ehrgott-ε- constraint (Algorithm 4) and the exact method proposed by Ehrgott (2006). In addition, the approximate Pareto fronts of all test problems are shown in Figure 4. It is worth noting that while both the exact and the approximate solutions of the first two test problems, found by Algorithms 3 and 4, respectively, are exactly the same, they slightly differ for the third and the fourth test problems.

Wu and Azarm (2001) discussed the quality of an observed Pareto solution set based on the two most important metrics including Hyper-area Difference and Accuracy of an Observed Pareto Frontier. For an observed Pareto frontier P, let HD(P), and AC(P) be the Hyperarea Difference and the Accuracy of an observed Pareto frontier respectively defined in Equations (52) and (53) ( Wu and Azarm, 2001):

$H D ( P ) = 1 − { ∑ r = 1 n p { ( − 1 ) r + 1 × [ ∑ k 1 = 1 n p − r + 1 ⋯ ∑ k l = k l − 1 + 1 n p − ( r − l + 1 ) + 1 ⋯ ∑ k r = k r − 1 + 1 n p ∏ i = 1 2 ( 1 − M a x j = 1 r ( Z ¯ i ( x k j ) ) ) ] } }$
(52)

$A C ( P ) = 1 H D ( P ) − S d o ( P )$
(53)

with(54)

$S d o ( P ) = ∑ r = 1 n p { ( − 1 ) r + 1 × [ ∑ k 1 = 1 n p − r + 1 ⋯ ∑ k l = k l − 1 + 1 n p − ( r − l + 1 ) + 1 ⋯ ∑ k r = k r − 1 + 1 n p ∏ i = 1 2 ( 1 − M i n j = 1 r ( Z ¯ i ( x k j ) ) ) ] }$
(54)

where $Z ¯ i ( x k ) = Z i ( x k ) − Z i I Z i N − Z i I$ and np is the number of observed Pareto solutions.

In general, Wu and Azarm (2001) concluded that a Pareto frontier with lower HD(P) and higher AC(P) values is considered to be better than the one with higher HD(P) and lower AC(P) values. Table 4 shows these quality metrics for small (test problems 1 & 2) and medium size (test problems 3 & 4) instances. Note how the values of the metrics pertaining to the Pareto solutions (both the exact and the approximate solution procedures) are close to each other. This is a good indicator of how well the approximate approach is performing. This result suggests that Algorithm 4 is outperforming other algorithms. Moreover, Table C1 (Appendix C) depicts the Pareto sets of small and medium size instances investigated for the bi-objective problem. Due to the memory and time limitation, although it is impossible to assess the specified metrics in other remaining test problems, the quality of the observed Pareto solution sets, which are related to the values of the binary decision variables, can be estimated using Algorithm 3.

## 7. CONCLUSION

In this paper, a multi-echelon P-D network problem that determines production and outsourcing plans, inventory plans, transportation plans, and distribution routing policies under uncertainty of products’ demands was investigated. A multi-stage stochastic mixed-integer programming model that incorporates several realistic constraints was developed. Furthermore, a safety stock policy was incorporated in the model, which transformed the model into a bi-objective optimization problem. A hybrid exact-approximate approach was first devised to solve the P-D problem instances in various sizes in the single objective case. The results suggested that the proposed algorithm resulted in better solutions compared to the solutions obtained by the GAMS solver in a shorter time. Then, a hybrid Ehrgott- ε-constraint algorithm was proposed to obtain an approximate Pareto front of the biobjective problem. The Pareto solutions of the small and medium size instances indicated that the proposed approximate solution approach performed very well in terms of two multi-objective metrics.

While this paper provides a framework for P-D planning of companies that are obliged to consider different market scenarios, the following topics are recommended for future studies:

• Developing risk-aversion and chance-constraint integrated P-D planning models,

• Taking transportation vehicles and their capacities into account, and

• Considering perishability of the products.

## ACKNOWLEDGMENT

The authors are thankful for constructive comments of the respected anonymous reviewers. Taking care of the comments improved the presentation significantly.

## Figure

Scheme of the proposed P-D problem.

Initial and reduced scenario tree of the first test problem.

Comparing the Pareto solutions of the first test problem solved Algorithms 4 (approximate) and the Ehrgott’s (2006) exact method.

Approximate Pareto solutions of all test problems.

## Table

The random cover generation algorithm

The cover evaluating algorithm

The hybrid exact-approximate GA-based algorithm to solve the multistage stochastic programming problem

The hybrid Ehrgott-ε-constraint algorithm

Test problem sizes

The detailed optimal solution of the first test problem

Comparing Algorithm 3 with GAMS solver to solve the multi-stage stochastic programming problem

Multi-objective quality metrics for small and medium size instances

Technological and cost parameters in test problems

Capacity parameters in test problems

The peak counts of waterbirds and raptors recorded at Geum River

Pareto solution sets of the small and medium size instances

## REFERENCES

1. Bard, J. F. and Nananukul, N. (2009), The integrated production- inventory-distribution-routing problem , Journal of Scheduling, 12(3), 257-280.
2. Beamon, B. M. (1998), Supply chain design and analysis: Models and methods , International Journal of Production Economics, 55(3), 281-294.
3. Boudia, M. , Louly, M. A. O. , and Prins, C. (2007), A reactive GRASP and path relinking for a combined production-distribution problem , Computers & Operations Research, 34(11), 3402-3419.
4. Chopra, S. and Meindl, P. (2007), Supply chain management, Strategy, planning & operation. Das SummaSummarum des Management, Book Chapter, Springer, 265-275.
5. Das, K. and Sengupta, S. (2009), A hierarchical process industry production-distribution planning model , International Journal of Production Economics, 117(2), 402-419.
6. Dasci, A. and Verter, V. (2001), A continuous model for production-distribution system design , European Journal of Operational Research, 129(2), 287-298.
7. Ehrgott, M. (2006), Multicriteria optimization (2nd Ed.), Springer Science & Business Media, Berlin.
8. Eichhorn, A. and Romisch, W. (2006), Mean-risk optimization models for electricity portfolio management , Proceedings of the 2006 IEE International Conference on Probabilistic Methods Applied to Power Systems, Stockholm, Sweden.
9. Esmaeilikia, M. , Fahimnia, B. , Sarkis, S. , Govindan, K. , Kumar, A. , and Mo, J. (2016), A tactical supply chain planning model with multiple flexibility options: Anempirical evaluation , Annals of Operations Research, 244(2), 429-454.
10. Fahimnia, B. , Luong, L. , and Marian, R. (2012), Genetic algorithm optimisation of an integrated aggregate production-distribution plan in supply chains , International Journal of Production Research, 50(1), 81-96.
11. Fahimnia, B. , Farahani, R. Z. , Marian, R. , and Luong, L. (2013), A review and critique on integrated production- distribution planning models and techniques , Journal of Manufacturing Systems, 32(1), 1-19.
12. Farahani, R. Z. and Elahipanah, M. (2008), A genetic algorithm to optimize the total cost and service level for just-in-time distribution in a supply chain , International Journal of Production Economics, 111(2), 229-243.
13. Ferris, M. C. (1998), MATLAB and GAMS: Interfacing optimization and visualization software, Mathematical Programming Technical Report 9819. http://www.cs.wisc.edu/math-prog/matlab.html.
14. Ferris, M. C. , Jain, R. , and Dirkse, S. (2011), Gdxmrw: Interfacing gams and matlab, Online: http://www.gams.com/dd/docs/tools/gdxmrw.pdf.
15. Hamedi, M. , Farahani, R. Z. , Husseini, M. M. , and Esmaeilian, G. R. (2009), A distribution planning model for natural gas supply chain: A case study , Energy Policy, 37(3), 799-812.
16. Heitsch, H. and R misch, W. (2003), Scenario reduction algorithms in stochastic programming , Computational Optimization and Applications, 24(2-3), 187-206.
17. Heitsch, H. and R misch, W. (2010), Stability and scenario trees for multistage stochastic programs, Stochastic Programming, Springer, Berlin, 139-164.
18. Jia, L. , Wang, Y. , and Fan, L. (2014), Multiobjective bilevel optimization for production-distribution planning problems using hybrid genetic algorithm , Integrated Computer-Aided Engineering, 21(1), 77-90.
19. Kanyalkar, A. and Adil, G. (2005), An integrated aggregate and detailed planning in a multi-site production environment using linear programming , International Journal of Production Research, 43(20), 4431-4454.
20. Kanyalkar, A. P. and Adil, G. K. (2010), A robust optimisation model for aggregate and detailed planning of a multi-site procurement-production-distribution system , International Journal of Production Research, 48(3), 635-656.
21. Kaut, M. and Wallace, S. W. (2003), Evaluation of scenario- generation methods for stochastic programming, Stochastic Programming E-Print Series (SPEPS), Julie L. Higle; Werner R misch; Surrajeet Sen, Editors, online: http://edoc.hu-berlin.de/series/speps/2003-14/PDF/14.pdf.
22. Kazemi, A. , Zarandi, M. F. , and Husseini, S. M. (2009), A multi-agent system to solve the productiondistribution planning problem for a supply chain: A genetic algorithm approach , The International Journal of Advanced Manufacturing Technology, 44(1-2), 180-193.
23. Kazemi Zanjani, M. , Nourelfath, M. , Ait-Kadi, D. (2010), A multi-stage stochastic programming approach for production planning with uncertainty in the quality of raw materials and demand , International Journal of Production Research, 48(16), 4701-4723.
24. Liu, S. and Papageorgiou, L. G. (2013), Multiobjective optimisation of production, distribution and capacity planning of global supply chains in the process industry , Omega, 41(2), 369-382.
25. Mohamed, Z. M. (1999), An integrated production-distribution model for a multi-national company operating under varying exchange rates , International Journal of Production Economics, 58(1), 81-92.
26. Nasiri, G. R. , Zolfaghari, R. , and Davoudpour, H. (2014), An integrated supply chain production-distribution planning with stochastic demands , Computers and Industrial Engineering, 77, 35-45.
27. Niknamfar, A. H. , Niaki, S. T. A. , and Pasandideh, S. H. R. (2015), Robust optimization approach for an aggregate production-distribution planning in a threelevel supply chain , The International Journal of Advanced Manufacturing Technology, 76(1-4), 623-634.
28. Nishi, T. , Konishi, M. , and Ago, M. (2007), A distributed decision-making system for integrated optimization of production scheduling and distribution for aluminum production line , Computers & Chemical Engineering, 31(10), 1205-1221.
29. Safaei, A. , Moattar Husseini, S. , Farahani, R. Z. , Jolai, F. , and Ghodsypour, S. (2010), Integrated multi-site production-distribution planning in supply chain by hybrid modeling , International Journal of Production Research, 48(14), 4043-4069.
30. Sarrafha, K. , Rahmati, S. H. A. , Niaki, S. T. A. , and Zaretalab, A. (2015), A bi-objective integrated procurement, production, and distribution planning problem of a multi-echelon supply chain network design: A new tuned MOEA , Computers and Operations Research, 54, 35-51.
31. Sel, C. , Bilgen, B. , Bloemhof-Ruwaard, J. , and van der Vorst, J. (2015), Multi-bucket optimization for integrated planning and scheduling in the perishable dairy supply chain , Computers & Chemical Engineering, 77, 59-73.
32. Wu, J. and Azarm, S. (2001), Metrics for quality assessment of a multiobjective design optimization solution set , ASME Journal of Mechanical Design, 123(1), 18-25.
33. Vira, C. and Haimes, Y. Y. (1983), Multiobjective decision making: Theory and methodology, North-Holland/ Elsevier Science Publishing Company, Inc., New York.
34. Yimer, A. D. and Demirli, K. (2010), A genetic approach to two-phase optimization of dynamic supply chain scheduling, Computers & Industrial Engineering, 58(3), 411-422.
35. Yoshida, O. , Nishi, T., and Zhang, G. (2014), Analysis of quantity discounts for multi-period production planning for single supplier and retailer under uncertain demands, Proceedings of the 2014 IEEE International Conference on Industrial Engineering and Engineering Management, Bandar Sunway, Malaysia.
36. Yu, V. F. , Normasari, N. M. E., and Luong, H. T. (2015), Integrated location-production-distribution planning in a multiproducts supply chain network design model, Mathematical Problems in Engineering, Article ID 473172, 13 pages.