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). Productiondistribution (PD) networks provide an effective tool to model manufacturing and logistics activities of a company (Dasci and Verter, 2001).
In this paper, a multistage stochastic mixed integer linear model is proposed for a complex PD 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 exactapproximate algorithm is used to obtain the solution. Next, the singleobjective optimization model is converted to a biobjective version by considering the safety stock policy; a case that has not been addressed in the literature. Another solution algorithm namely, the epsilonconstraint method, is utilized to solve the proposed biobjective 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 PD planning problem. Comprehensive literature reviews of previous studies on PD 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 multiitem production, parallel multiplant production, and dynamic distribution. They also demonstrated the capability of the proposed model in coordinating the production and distribution planning of a company with multilocation parallel manufacturing facilities in a consumergoods 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 metaheuristics algorithm was employed as the solution method and was improved on through a reactive mechanism alongside a pathrelinking process. A global and multi national productiondistribution 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 subproblems. The scheduling subproblem was solved using a simulated annealing method and the warehousing subproblem was solved using the commercial mixed integer linear programming (MILP) solver of the CPLEX software. Note that this study did not consider endusers and demand fluctuations. A justintime distribution supply chain network including several products and a single plant was considered by Farahani and Elahipanah (2008). They proposed a biobjective mixedinteger linear model to minimize both the network costs and the sum of backorders and surpluses of products in all periods, simultaneously. A hybrid nondominated sorting genetic algorithm was utilized in their work to solve realsize instances. Hamedi et al. (2009) introduced a sixlevel 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 multiagent system of a PD 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 GAbased solution procedure was developed by Yimer and Demirli (2010) to solve the scheduling problem of a builttoorder supply chain. Moreover, a twoechelon 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 mediumsize instances, it failed in solving largesize problems.
In a study by Kanyalkar and Adil (2010), a robust optimization model for an integrated multisite 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 multiproduct, multiperiod, multisite production and distribution planning problem. However, their mixed integer linear program model ignored many characteristics of the realworld (e.g. backlogging costs, production costs, and detail of production). Nasiri et al. (2014) proposed a hierarchical model to design an integrated productiondistribution network of a threeechelon supply chain with stochastic demands. They incorporated a capacitated facility location model in order to determine the locations of multilevel 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 multiproduct, multiechelon, multiperiod 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 multiobjective mixedinteger 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 multiperiod bilevel program. A multiobjective bilevel PD 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, multiproduct 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 realworld sawmill. In addition, a scenario based model of an aggregate PD 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 mixedinteger 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 multiobjective biogeographybased optimization (MOBBO) algorithm with tuned parameters in order to obtain a nearoptimum solution. Finally, using mixedinteger 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 worstbased scenario

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

3) It tackles the problem using a metaheuristic 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 PD network in which a set of manufacturers denoted by $\{1,\hspace{0.17em}\dots ,\hspace{0.17em}\text{}M1,\hspace{0.17em}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 submanufacturers incorporating the set of $\{1,\hspace{0.17em}\dots ,\hspace{0.17em}\text{}M1\}$, and the main manufacturer denoted by {M}. There are s rawmaterial suppliers who provide the primary requirements of k types of rawmaterials of both subs and main manufacturers. Each of the components is produced solely in a unique submanufacturer. Furthermore, each of the submanufacturers 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 submanufacturers is able to produce it.
In order to meet the demand, the main manufacturer delivers finished products to e enddistributors 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 PD problem is illustrated in Figure 1. Binary decision variables at the beginning of the planning horizon determine the selected warehouses and their corresponding enddistributors. 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 submanufacturer in regular and overtime periods,

 The total units to be outsourced by the main manufacturer,

 The total units of components to be conveyed from submanufacturers 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 enddistributors.
As it can be seen in Figure 1, the location of the warehouses and their flow between enddistributors are characterized by dashed lines to show their possible selections. In other words, after determining the warehouses to be constructed, each of the enddistributors 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 scenariogeneration 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 nonanticipative, 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 enddistributors 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 ${\sigma}_{\epsilon}^{2}$ , with the assumption that all demand realizations are within the interval $({\mu}_{i}3{\sigma}_{e},\hspace{0.17em}{\mu}_{i}+3{\sigma}_{e}).$. This interval is then partitioned to subintervals $({\mu}_{i}3{\sigma}_{e},\hspace{0.17em}{\mu}_{i}2{\sigma}_{e}),\hspace{0.17em}({\mu}_{i}2{\sigma}_{e},\hspace{0.17em}{\mu}_{i}{\sigma}_{e}),$ $({\mu}_{i}{\sigma}_{e},\hspace{0.17em}{\mu}_{i}+{\sigma}_{e}),\hspace{0.17em}({\mu}_{i}+{\sigma}_{e},\hspace{0.17em}{\mu}_{i}+2{\sigma}_{e}),\hspace{0.17em}({\mu}_{i}+2{\sigma}_{e},\hspace{0.17em}{\mu}_{i}+3{\sigma}_{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 PD 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 realworld 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 PD 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 ($\text{\mathcal{M}}=\{(t,\hspace{0.17em}n)t=1,\hspace{0.17em}\dots ,\hspace{0.17em}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 enddistributors

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 submanufacturers’ plants

B Set of components which outsourced by the main manufacturer

C_{m} Set of raw materials which are used in manufacturer m plant

D_{s} Set of raw materials provided by supplier s

${\text{\mathcal{E}}}_{\text{}}$ Set of suppliers that can supply raw material k

${\text{\mathcal{F}}}_{\text{}}{}_{i}$ Set of raw materials to be used in the assembly of product i

${\text{\mathcal{G}}}_{\text{\U0001d4df}}$ Set of raw materials to be used in the production of component p

${\text{\mathcal{H}}}^{m}$ Set of components that can be produced by submanufacturer m
4.2.2. Parameters
The parameters involved in the mathematical formulation of the problem are:

D_{ietn} Demand for product i by enddistributor e during time period t and node n

pr_{n} Occurrence probability of node n

F_{w} Fixed costs of opening and operating warehouse w

${H}_{pt}^{w}$ Unit holding cost for component p during period t in submanufacturers’ plants

${H}_{pt}^{wM}$ Unit holding cost of component p during period t in the main manufacturer’s plant

H_{iwt} Unit holding cost for product i in warehouse w during period t

${H}_{iet}^{\text{'}}$ Unit holding cost for product i in enddistributor e during period t

${H}_{kmt}^{R}$ Unit holding cost for raw material k in manufacturer m plant during period t

TR_{ksmt} Unit transportation cost for raw material k from supplier s to manufacturer m during period t

T_{iwt} Unit transportation cost of product i from the main manufacturer to warehouse w during period t

${T}_{iwet}^{\text{'}}$ Unit transportation cost of product i from warehouse w to enddistributor e during period t

TV_{pt} Unit transportation cost for component p from its unique submanufacturer to the main manufacturer during period t

a_{pi} Quantity of component p required per unit of product i

τ_{kp} Amount of raw material k used to produce a unit of component p

${\tau}_{kt}^{\text{'}}$ Unit cost of raw material k provided by raw material suppliers during period t

${\tau}_{ki}^{M}$ Amount of raw material k used to produce a unit of product i in the main manufacturer’s plant

${\upsilon}_{p}$ Processing time to produce a unit of component p

${\sigma}_{pt}$ Labor/hour cost for regulartime production of component p during period t

${\sigma}_{pt}^{\text{'}}$ Labor/hour cost for overtime production of component p during period t

${\upsilon}_{i}^{M}$ Processing time to produce a unit of product i

${\sigma}_{it}^{M}$ Labor/hour cost for regulartime production of product i during period t

${\sigma}_{it}^{\text{'}M}$ Labor/hour cost for overtime production of product i during period t

G_{pt} Unit outsourcing cost of component p during period t

L_{it} Minimum production quantity of product i during period t

K_{pt} Minimum production quantity of component p during period t

Q_{wt} Holding capacity of warehouse w during period t

${Q}_{et}^{\text{'}}$ Holding capacity of enddistributor e during period t

η_{i} Empty space assigned to a unit of product i

${\eta}_{p}^{\text{'}}$ Space occupied by a unit of component p

ρ_{k} Space occupied by a unit of raw material k

${W}_{mt}^{Max}$ Holding capacity for produced components in submanufacturer m in period t

${W}_{t}^{MaxM}$ Holding capacity of products’ components in the main manufacturer plant in period t

${R}_{mt}^{Max}$ Holding capacity of raw materials at manufacturer m plant during period t

$C{P}_{mt}^{Max}$ Transportation capacity of components produced by submanufacturer m during period t

$Q{R}_{st}^{Max}$ Capacity of raw materials supplied by supplier s during period t

${\lambda}_{mt}$ Hours capacity for regulartime production at manufacturer m during period t

${\lambda}_{mt}^{\text{'}}$ Hours capacity for overtime production at manufacturer m during period t

${\lambda}_{t}^{M}$ Hours capacity for regulartime production at the main manufacturer during period t

${\lambda}_{t}^{\text{'}M}$ Hours capacity for overtime production at the main manufacturer during period t

${\gamma}_{p1}$ Inventory level of component p on its own submanufacturer plant at the beginning of the planning horizon(t = 1)

${\gamma}_{p1}^{M}$ Inventory level of component p in the main manufacturer's plant at the beginning of the planning horizon(t = 1)

${\theta}_{iw1}$ Inventory level of product i in warehouse w at the beginning of the planning horizon (t = 1)

${\epsilon}_{ie1}$ Inventory level of product i in enddistributor e at the beginning of the planning horizon (t = 1)

${\kappa}_{km1}$ 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:

I_{itn} The quantity of product i produced in regular time at the main manufacturer’s plant during period t at node n

${I}_{itn}^{\text{'}}$ The quantity of product i produced in overtime at the main manufacturer’s plant during period t at node n

U_{ptn} The quantity of component p produced in regulartime during period t at node n

${U}_{ptn}^{\text{'}}$ The quantity of component p produced in overtime during period t at node n

${U}_{ptn}^{\text{'}\text{'}}$ The quantity of component p outsourced by the main manufacturer during period t at node n

W_{ptn} The Inventory level of component p at the beginning of period t at node n

${W}_{ptn}^{M}$ The Inventory level of component p at the main manufacturer at the beginning of period t at node n

${R}_{kmtn}$ The Inventory level of raw material k at manufacturer m at the beginning of period t at node n

${V}_{iwtn}$ The Inventory level of product i at warehouse w at the beginning of period t at node n

${V}_{ietn}^{\text{'}}$ The Inventory level of product i at enddistributor e at the beginning of period t at node n

X_{ksmtn} The amount of raw material k transmitted from supplier s to manufacturer m during period t at node n

J_{iwtn} The quantity of product i transmitted from the main manufacturer to warehouse w during period t at node n

${J}_{iwetn}^{\text{'}}$ The amount of product i transmitted from warehouse w to enddistributor e during period t at node n

Y_{ptn} The amount of component p transmitted from its unique manufacturer during period t at node n

${O}_{w}=\{\begin{array}{l}1\text{ifwarehouse}w\text{isopen}\\ \text{0}\text{otherwise}\end{array}$
${O}_{we}^{\text{'}}=\{\begin{array}{l}1\text{ifenddistribuotr}e\text{issupportedbywarehouse}w\\ 0\text{otherwise}\end{array}$
4.2.4. The MultiStage Stochastic Programming Model
Using the indices, parameters and decision variables defined in Section 4.2.3, the multistage stochastic programming model of the problem is derived as:(5)(6)(29)(30)(31)(33)(34)(35)(36)(37)(38)(39)(40)(41)(42)
s.t.
such that
$\begin{array}{l}{\text{\mathcal{M}}}^{\text{'}}(\subset \text{\mathcal{M}})=\{(t,\hspace{0.17em}n)t=1,\hspace{0.17em}\dots ,\hspace{0.17em}T1\}\\ {\text{\mathcal{M}}}^{\text{'}\text{'}}(\subset \text{\mathcal{M}})=\{(t,\hspace{0.17em}n)t=2,\hspace{0.17em}\dots ,\hspace{0.17em}T\}\end{array}$
The objective function in (1) aims to minimize all associated costs of the investigated PD 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 enddistributors 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 enddistributor. 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 enddistributors, 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 submanufacturers. 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 submanufacturer to the main manufacturer is controlled by constraint set (25). Constraint sets (26) and (27) are the products holding capacity in warehouses and enddistributors, respectively. Constraints (28)(32) limit the inventory level of raw material, components, and finished products at sub and mainmanufacturers, warehouses, and enddistributors, respectively. Finally, the remaining constraints define the types of the variables used.
4.3. The BiObjective Model
In realworld 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 enddistributors 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)
where $\left.\right$ is the absolute value function. As seen, the above function is nonlinear but can be linearized using the following:(44)
Note that the parameter SS_{ie} is the safety stock of product i at enddistributor e which is mandated by the main decision maker to the model developer. Based on the above definition, we are minimizing:(45)(46)
where α_{1}, α_{2} are the weights of EX_{ietn} and SH_{ietn}, respectively.
5. SOLUTION METHOD
In this section, a hybrid exactapproximate 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 biobjective optimization problem modeled in Section 4.3.
5.1. The Hybrid ExactApproximate 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 1012. 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 GAbased 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 biobjective optimization problem is shown in (47)
where ${Z}_{k},\hspace{0.17em}k=1,\hspace{0.17em}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, $\widehat{x}\in \text{},$, is called the Pareto optimal solution if there is no $x\in \text{}$ such that ${f}_{k}(x)\le {f}_{k}(\widehat{x})$ for k = 1, 2 and at least for one $j\in \{1,\hspace{0.17em}2\}\hspace{0.17em}we\hspace{0.17em}have\text{}\hspace{0.17em}{f}_{j}(x)<{f}_{j}(\widehat{x})$. If $\widehat{x}$is the Pareto solution, then $f(\widehat{x})$ is called nondominated point. The set of all Pareto solutions is denoted by ${\text{}}_{E}$ and is called efficient set. The set of all nondominated points is denoted by Y_{N} and called the Pareto front or the nondominated set.
Definition 5.2: The point ${y}^{I}=({y}_{1}^{I},\hspace{0.17em}{y}_{2}^{I})$ given by(48)
is called the Ideal point of the biobjective optimization problem. Moreover, the point ${y}^{N}=({y}_{1}^{N},\hspace{0.17em}{y}_{2}^{N})$ given by(49)
is called the Nadir point of the biobjective optimization problem. Ehrgott (2006) presented an algorithm to obtain ideal and nadir points of the above biobjective optimization problems.
5.2.1. The εConstraint Method
The εconstraint method, firstly introduced by Vira and Haimes (1983), is probably the bestknown technique for solving multiobjective optimization problems. In this method, one of the original objectives is minimized, while the others are transformed to constraints (Ehrgott, 2006). The biobjective optimization problem is substituted by the following εconstraint problem:(50)
where $i,j\in \{1,2\}$ and i ≠ j.
Algorithm 4 obtains an approximate Pareto front for the biobjective optimization problem (47) by solving ( ) i j P ε problems iteratively. In this algorithm, $\Delta =\frac{{Z}_{2}^{\ast N}{Z}_{2}^{\ast 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 biobjective 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 ExactApproximate Solution Method
In order to demonstrate the efficacy of the proposed exactapproximate 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 submanufacturer, respectively.
The performance of the proposed algorithm is compared with the mixedinteger solver of GAMS in terms of percentage deviation defined as (Farahani and Elahipanah, 2008):(51)
where Z_{GA} is the objective function value found by Algorithm 3 and Z_{GAMS} 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) i74500U 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 nearoptimum 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 Hyperarea 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):
with(54)
where ${\overline{Z}}_{i}({x}_{k})=\frac{{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 biobjective 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 multiechelon PD network problem that determines production and outsourcing plans, inventory plans, transportation plans, and distribution routing policies under uncertainty of products’ demands was investigated. A multistage stochastic mixedinteger 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 biobjective optimization problem. A hybrid exactapproximate approach was first devised to solve the PD 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 multiobjective metrics.
While this paper provides a framework for PD planning of companies that are obliged to consider different market scenarios, the following topics are recommended for future studies: