Optimization of a Dynamic Supply Portfolio Considering Risks and Discount ’ s Constraints

Purpose: Nowadays finding reliable suppliers in the global supply chains has become so important for success, because reliable suppliers would lead to a reliable supply and besides that orders of customer are met effectively. Yet, there is little empirical evidence to support this view, hence the purpose of this paper is to fill this need by considering risk in order to find the optimum supply portfolio. Design/methodology/approach: This paper proposes a multi objective model for the supplier selection portfolio problem that uses conditional value at risk (CVaR) criteria to control the risks of delayed, disrupted and defected supplies via scenario analysis. Also we consider discount’s constraints which are common assumptions in supplier selection problems. The proposed approach is capable of determining the optimal supply portfolio by calculating valueat-risk and minimizing conditional value-at-risk. In this study the Reservation Level driven Tchebycheff Procedure (RLTP) which is one of the reference point methods, is used to solve small size of our model through coding in GAMS. As our model is NP-hard; a meta-heuristic approach, Non-dominated Sorting Genetic Algorithm (NSGA) which is one of the most efficient methods for optimizing multi objective models, is applied to solve large scales of our model.


Introduction and motivation
Supply management deals with optimization of supply portfolio, which consists of supplier selection and order allocation.Optimizing the allocation of orders on the part of the buyer in a multi-supplier environment has become a major concern in supply chains (Moliné & Coves, 2013).To survive in the increasingly fierce global competition, firms require having cooperative and reliable relationship with their suppliers, so supplier selection is an important issue in supply management (Lee, 2009a).Supplier selection consists of recognition and assessment of suitable suppliers (Ravindran, Bilsel, Wadhwa & Yang, 2010).In the selection process it is not reasonable to consider the cost lonely, because cost-based approach cannot guarantee quality and on time delivery.In contemporary supply chain management, the performance of potential suppliers is evaluated against multiple criteria rather than considering a single factor-cost.For a suitable choice, many factors are involved such as price, quality, reliability, flexibility, and so on.These factors may contradict each other and may not be satisfied simultaneously (Sawik, 2011a;Ho, Xu & Dey, 2010).
In a risky make-to-order environment, customer-oriented manufacturers should be prepared to produce varieties of products to meet different customer needs.Each product is typically composed of many common and non-common (custom) parts that can be sourced from different approved suppliers with different supply capacities.An important issue is how to allocate the best order of parts among various suppliers so as to fulfill all customer orders of products and to achieve a high customer service level at a low cost and, in addition, to alleviate the impact of supply chain risks.The supply chain risk management has been extensively studied over the past decade (Sawik, 2011b).According to the literature, there are two risk levels: operational risks and disruption risks (Tang, 2006).Operational risks related to the uncertainty such as uncertain customer demand.Disruption risks could be originated from natural disasters (such as earthquakes, etc.) or manmade disasters (such as labor strikes, terrorist attacks, etc.) and caused interruption or failure of the activities.Disruption risks may be seriously detrimental to economy rather than operational risks (Sawik, 2011b).Knemeyer, Zinn and Eroglu (2009) considered a proactive planning for catastrophic events in supply chains, based on an innovative methodology which is used by the insurance industry to quantify the risk of multiple types of catastrophic events on the key supply chain locations.
Despite the stochastic nature of supplier selection problem, a few researches consider uncertainty and risk.Selecting the supply portfolio under variety of risks is a hard discrete stochastic optimization problem.Risk measurement is a key concept in risk management.Two popular methods of risk measurement, which have been used extensively in financial engineering, are value-at-risk (VaR) and conditional value-at-risk (CVaR).VaR estimates maximum loss when cost distribution is symmetric but CVaR determines the maximum loss for undesirable conditions when cost distribution is asymmetric (Sawik, 2011a).Wu and Olson (2010) combined DEA with VaR to improve a performance measurement system and presented a new approach for selection of vendors in enterprise risk management (ERM).
They used Monte Carlo Simulation for benchmarking their proposed method.Ravindran et al. (2010) developed multi-criteria supplier selection models with considering two types of risk measures, value-at-risk (VaR) and miss-the-target (MtT).Their model had two phases, in phase one; suppliers were ranked by AHP and key suppliers were selected.In phase two, the orders of quantities are allocated among different suppliers with using the presented model.Lee (2009a) proposed a fuzzy analytic hierarchy process model (FAHP) for evaluating buyersupplier relation with considering benefits, opportunities, costs and risks concepts.Due to the uncertainty in the real-life, fuzzy set theory was used in this model.Also Lee (2009b) suggested FAHP model for supplier selection and applied it for a TFT-LCD manufacturer as a case study.Wu, Zhang, Wu and Olson (2010) presented a fuzzy multi-objective programming for supplier selection with considering quantitative (costs, quality acceptance levels, and ontime delivery distributions) and qualitative (economic, environment and vendor rating) supplier selection risk factors.Xiao, Chen and Li (2012) studied a supplier selection problem under operational risk, for this purpose they extended fuzzy concepts.Sawik (2010) introduced a mixed integer model for supplier selection for make-to-order manufacturing.For controlling risks, the author set an upper limit for average defect rate and late delivery rate.Discount policies are usually considered between buyers and suppliers.Discount offered by the supplier is considered in the model provided by Sawik (2010).Also Sawik (2011c) controlled operational risk by Value-at-Risk (VaR) and Conditional Value-at-Risk (CVaR) for problem of supplier selection in make-to-order environment.Then in (Sawik, 2011b), the portfolio approach was enhanced to consider a single-period supplier selection and supply chain disruption risks.The author controlled this risk by VaR and CVaR concepts.Sawik (2012) considered disruption risks in supplier selection portfolio.He presented models that deal with the optimal selection and protection of part suppliers and order quantity allocation in a supply chain with disruption risks.The objective of model was to achieve a minimum cost of suppliers' protection, emergency inventory pre-positioning, parts ordering, purchasing, transportation and shortage so as to mitigate the impact of disruption risks by minimizing the potential worstcase cost.
The models were developed for supplier selection and order allocation could be either singleperiod models that do not consider inventory or multi-period models which consider the inventory; however it is obvious that in make-to-order and multi-period environment there is no need to consider inventory.Ustun and Demirtas (2008) studied a multi-objective and multiperiod supplier selection problem.For selecting best supplier and optimal allocation, they evaluated all suppliers with fourteen criteria via ANP technique and the multi objective model was solved by multi-objective optimization technique.Sawik (2011a) studied a multi-period supplier selection in the presence of disruption and delay risks.Moreover, the author used VaR and CVaR via scenario analysis.
Due to a major gap between the theoretical solutions of multi-product multi-period supplier selection and real world problems, we intend to provide a multi-product, multi-period model with discount's constraints as a portfolio problem that is more realistic than the previous.
Hence, this paper presents a mixed integer, multi-period and multiple objectives programming model for supplier selection for a dynamic supply portfolio and incorporates risks via scenario analysis.The proposed model will make a simple tool for decision maker (DM) for evaluating suppliers to optimize the supply portfolio.Moreover, time value of money due to interest rate is considered for the expected cost objective.A Non-dominated Sorting Genetic Algorithm (NSGA), which is one of the most efficient methods for optimizing multi objective models, is presented to solve large scales of our model.Also Reservation Level driven Tchebycheff Procedure (RLTP), which is one of the reference point methods, is used to solve small size of our model.To demonstrate the capability of our model numerical examples are generated.This paper is organized as follows.Related works are presented in the next section.Problem definition and the proposed model are discussed in sections 3 and 4 respectively.Solution methods are suggested in section 6. Numerical examples and some computational results are mentioned in section 7.And finally conclusions and future research directions are made in the last section.

Related works
Meena, Sarmah and Sarkar (2011) studied a supplier selection problem under risks of supplier failure due to the catastrophic events disruption.An analytical model is developed to determine the optimal number of suppliers considering different failure probability, capacity, and compensation.Yu, Zeng and Zhao (2009) focused on evaluating the impacts of supply disruption risks on the choice between the famous single and dual sourcing methods in a twostage supply chain with a non-stationary and price-sensitive demand.Xu and Nozick (2009) formulate a two-stage stochastic program and a solution procedure to optimize supplier selection to hedge against disruptions.Their model allows for the effective quantitative exploration of the trade-off between cost and risks to support improved decision-making in global supply chain design.Tsai and Wang (2010) applied a mixed integer programming approach to solve the sourcing and order allocation problem with multiple products and multiple suppliers in a supply chain; also two schemes of quantity discounts are used to compare the influence upon the buying decisions.Xia and Wu (2007) developed an integrated approach of analytical hierarchy process improved by rough sets theory and multi-objective mixed integer programming in order to simultaneously determine the number of suppliers to employ and the order quantity allocated to these suppliers in the case of multiple sourcing, multiple products, with multiple criteria and with supplier's capacity constraints.Authors considered that suppliers may offer price discounts on total business volume, not on the quantity or variety of products purchased from them.Stadtler (2007) presented a linear mixed integer programming (MIP) model, which not only represents the all-units discount but also the incremental discount case.Meena and Sarmah (2013) investigated an order allocation problem of a manufacturer/buyer among multiple suppliers under the risks of supply disruption.A mixed integer non-linear programming (MINLP) model is developed for order allocation considering different capacity, failure probability and quantity discounts for each supplier.Sawik (2011a) considered the problem of a multi period supplier selection and order allocation in make to order environment in the presence of supply disruption.In this regard the author presented a mixed integer programming model in order to find the optimal dynamic supply portfolio.Sawik (2011a) has suggested that the proposed model can be enhanced for a discount environment, where the suppliers offer discounts based on quantity or business volume of ordered parts.Sawik (2010) formulated similar problem of allocating orders for custom parts among suppliers in make to order manufacturing as a mixed integer program model which considers the quantity and business volume discounts offered by suppliers.The proposed model in (Sawik, 2010) considered a single period supplier selection and order allocation problem which provides a static supply portfolio.Sawik (2010) has mentioned to extend the model in a dynamic case where orders arrive irregularly over time.
The proposed model in this paper is a combination of two models which developed by (Sawik, 2010) and (Sawik, 2011a).Indeed, our model is a multi-period supplier selection and order allocation problem which finds the optimal dynamic supply portfolio, and also considers the quantity and business volume discounts offered by suppliers as well.Moreover, time value of money due to interest rate is considered.The developed model and assumptions is described in detail next.

Problem description
In this paper, we consider a supply chain in which various types of products are assembled by a single producer and also satisfy customer orders by providing different part types from multiple suppliers.In this supply chain, we will focus on the supplier selection portfolio which is a tactical decision and is made in each period.Due to the multi-period nature of the model, our portfolio will be dynamic.It is not bad to describe a little about the differences between a dynamic and a static supply portfolio.A dynamic supply portfolio is the allocation of orders among the suppliers combined with the allocation of orders among the planning periods, but a static supply portfolio is the allocation of orders for parts among the suppliers without allocation of them among the planning periods (Sawik, 2011a).
It is worthy to assume that the delivery date and the corresponding reliability of on-time delivery of each supplier may randomly vary in different periods.With respect to requested dates, different suppliers may deliver the ordered parts with different delays and a different operational risk can be associated with each supply portfolio.So the disruption risk is incorporated utilizing the concepts of VaR and CVaR.The portfolio is optimized by calculating VaR and minimizing CVaR simultaneously.For a finite number of scenarios, CVaR allows the evaluation of worst-case costs and shaping of the resulting cost distribution through optimal supplier selection and order allocation decisions, i.e. the selection of optimal supply portfolio (Sawik, 2011a).Likewise, the quality of parts delivered by each supplier may randomly vary and each supplier may have different capacity with different discount offers.When the suppliers are selected the risk of defective and unreliable (late) deliveries can be considered using past observations.The supply delays may result in the shortage of required parts and the corresponding delay penalty costs of delayed customer orders will be incorporated into the model.

The proposed model
Let I={1,…, M} be the set of M suppliers and J={1,…, N} the set of N part types which are required for the products, this is known ahead of time and not varying over the time horizon.
Denote by dj the demand for each part type j є J.The planning horizon consists of H planning periods (e.g.days or weeks) and let T={1,…, H} be the set of planning periods.Assume that supplies are subject to random disruptions.So supplies will be delayed or disrupted (never delivered) in some periods.
The procedure of Sawik (2011a) for calculating delays is as follows.On-time delivery scheduled for period t occurs in that period, whereas the late delivery may occur in one of the remaining periods t+1,…, H, i.e.To deliver the ordered parts by the end of the planning horizon the delay cannot be longer than H-t periods.By convention, deliveries delayed until a dummy period t=H+1 represents disruptions of supplies, i.e. no delivery of ordered parts.
Summarizing, it is assumed that the longest delay of each delivery scheduled for period t cannot be greater than H-t periods, whereas the infeasible delay of H-t+1 periods (i.e.delivery in a dummy period t=H+1), by convention represents a disruption, that is, no delivery.
Let us call delivery scenarios by l є L = {1,…, L} so Δ l ijt which is an integer vector, Δ l ijt є {0,...,H-t+1}, this represents that delivery of part type j from supplier i є I scheduled for period t є T is either on-time, if Δ l ijt = 0, is delayed by Δ l ijt periods, if 1 ≤ Δ l ijt ≤ H-t or is disrupted, if Δ l ijt = H-t+1.The probability that scenario l for part type j and supplier i occurs in period t is denoted by П l ijt.We assume that the last scenario in the set of scenarios represents the probability of occurrence the disruption (i.e.scenario L).It is obvious that: Notations are shown in Table 1.
Now we propose a mixed integer multi objective programming model for selection of a supply portfolio in the presence of supply delays and disruptions with discounts constraints.As mentioned before, the DM needs to select a dynamic supply portfolio, i.e. the allocation of orders for parts among the suppliers and among the planning periods.A dynamic supply portfolio in non-discount environment is defined as: and 0 ≤ Fit ≤ 1 is the fraction of the total demand for parts (all parts) ordered from supplier i in period t.Δ l ijt Difference between scheduled delivery time (period t) and realized delivery time of part type j ordered by supplier i П l ijt Disruption and delay probability for part type j scheduled to be supplied from supplier i in period t ξijtw Price discount rate (percentage of discount) associated with discount interval w of supplier i for part type j in period t bijtw Upper limit on discount interval w of supplier i for part j in period t β Annual interest rate

Decision Variables
VaR The targeted cost of portfolio based on the ∝-percentile of costs, i.e. in 100∝% of scenarios, the costs cannot exceed VaR (value-at-risk)

Fit
The fraction of total demand for parts ordered from supplier i to be delivered in period t (portfolio variable) xijt The fraction of total demand for part type j ordered from supplier i to be delivered in period t (order allocation variable) yijt 1, if an order for part j to be delivered in period t is placed on supplier i; otherwise yijt=0 (supplier selection variables) zijt 1, if part type j is ordered from supplier i in period t and total business volume (or total quantity) of that part type purchased from this supplier falls on the discount interval w; otherwise Zijtw=0 (customer order allocation variable) Uijtw 1, if total business volume (or total quantity) part type j from supplier i falls on the discount interval w in period t; otherwise = 0 T' Tail cost for delivery scenario l (auxiliary variable to make CVaR equation linear)

Qijtw
Auxiliary variable for making purchasing cost term in a linear expression  Sawik (2011a) suggested that the dynamic supply portfolio should be checked over the time horizon against the highest acceptable defect rate, delay rate and disruption rate of supplies.Now, denote by q, r and s the maximal acceptable average defect rate, delay rate and disruption rate, respectively.We will trace this suggestion in our model.
In make to order environment, in which custom parts are typically ordered in small lot sizes, supplier may sometimes offer discounts that depends on total value of sales volume (business volume) or on total quantity of ordered parts.In the context of business volume, discount the quantity or variety of purchased parts does not affect the offered price, while for the quantity discount the price does not depend on the total amount of sales volume (Sawik, 2010).
The procedure of Sawik (2010) for considering discount constraints is as follows.Assume that each supplier i offers cumulative (all-units) price breaks having bijtw discount intervals according to the total business volume (or the total quantity), (bijtw-1, bijtw], where bijtw is upper limit on the wth business volume (or quantity) discount interval for supplier i and part j in period t and bijt0=0 for all iєI, jєJ and tєT.Let 0<ξijtw<1 be the discount rate (percentage of discount) associated with interval w of supplier i for part type j in period t.If total business volume (total quantity) from supplier i falls on interval w in period t, then the price of each part type j is ((1-ξijtw).pijt).The following discount interval selection variable needs to be added to enhance the supplier selection and to order allocation problem for the discount environment (Sawik, 2010): Uijtw=1 if total business volume (or total quantity) for part type j from supplier i falls on the discount interval w in period t; otherwise Uijtw=0.
Furthermore, the order assignment variable Zijwt needs to be clarified as follows: Zijtw = 1, if parts type j are ordered from supplier i and total business volume (or total quantity) of parts type j purchased from this supplier falls on the discount interval w in period t; otherwise Zijtw = 0.
It is obvious that in make-to-order manufacturing, no inventory of custom parts can be kept on hand and the parts are requisitioned with each customer order.So inventory balance constraints do not need to be considered in the model.In addition, we can assume that parts would be delivered within a time window [ej,fj] and derived from the customer requested due date to reduce any further inventory holding cost.The delivery for part type j cannot be earlier than the earliest date ej.On the other hand, the required parts type j can be delivered late, beyond the latest date fj, or not delivered at all.Then, delay or disruption penalty costs will be charged.
Let aj be per unit and per period penalty cost of delayed part type j caused by the late delivery, where deliveries which are not later than fj, are not penalized.Similarly, let Bj be per unit and per period cost of unfulfilled order for part type j, caused by supply disruptions for part type j.
For delivery scenarios form (l=1) to (l=L-1) and each part type in each period deliveries which are later than fj, are penalized by aj.We can assume that fj is always smaller than H.The total expected delay penalty cost overall delivery scenarios, given an order allocation vector Z, is then given by: Note that in equation ( 3) summations have two following limitations: Summation on t is only for those parts which their fj are smaller than the sum of the time scheduled for deliveries of them and its corresponding delay of concerned scenario, and Summation on l is only for scenarios that their corresponding delays are between ( 1) to (H-t).
i.e. those scenarios that show delays not disruptions.
Also it is obvious that if fj=H, then no delivery within the planning horizon can be delayed and penalized.However, any delay beyond H, by convention is a disruption, and hence penalized with the cost of unfulfilled orders (Bj).The total expected penalty cost of unfulfilled orders for parts due to supply disruptions (i.e.deliveries delayed until period H+1) is: The total business volume or total quantity of parts purchased from supplier i is, It is clear that these two expressions are nonlinear and will be made linear by fifth part of constraints in the next section.Notice that the total business volume or the total quantity purchased from each supplier i falls on exactly one discount interval of this supplier, and their phrases are nonlinear that will be reformed to a linear phrase.Now, the average purchasing and ordering cost for all demand (D) can be expressed as follows: It is clear that y, Z and x are decision variables.The second part of equation ( 5) is a nonlinear term so it will be changed to linear term by some constraints and will be mentioned in the next part.Here the performance of the dynamic portfolio with discount in expected costs section can be measured by the sum F(y, Z, x) of average cost per part of ordering, purchasing, defects, delays and disruptions, respectively, in which the cost of a defective parts is assumed to be identical with its regular price: Now the procedure of evaluating the performance of the dynamic supply portfolio in risk section could be clarified.First it is necessary to discuss a little about main general concepts of VaR and CVaR, so for more details interested readers are referred to (Rockafellar & Uryasev, 2002;Uryasev, 2000).
Let αє(0,1) represents the confidence level for the cost distribution across all scenarios.The following two risk measures are commonly used in portfolio selection problems.
VaR at a 100 α% confidence level is the targeted cost of the portfolio such that for 100 α% of the scenarios, the costs will not exceed VaR.In other words, VaR is a parameter (fixed by the DM) or a decision variable based on the α-percentile of costs, i.e. in 100(1-α) % of the scenarios, the costs may exceed VaR.For instance, 90%-VaR is an upper estimation of losses which is exceeded with 10% probability (Rockafellar & Uryasev, 2002;Uryasev, 2000).
CVaR at a 100 α% confidence level is the expected cost of the portfolio in the worst 100(1-α)% of the cases.In other words, we allow 100(1-α) % of the outcomes to exceed VaR, and the mean value of these outcomes is represented by CVaR.So α-conditional Value-at-Risk (α-CVaR) is the minimizing of "the expected value of the costs in the (1-α)100% worst cases" (Rockafellar & Uryasev, 2002;Uryasev, 2000).
VaR represents the maximum cost associated with a specified confidence level of outcomes (i.e. the likelihood that a given portfolio's cost will not exceed the defined amount of VaR).
However, VaR does not explain the magnitude of the cost when the VaR limit is exceeded.In other words, VaR is the acceptable cost level above which we want to minimize the number of realizations and CVaR considers those portfolio outcomes, where costs exceed VaR.CVaR also can be considered as a weighted measure of VaR and the costs above VaR, which may be extremely high costs.When using CVaR to minimize worst-case costs, CVaR is always not less than VaR (Sawik, 2011a).
So we can conclude that the CVaR is a superior than the VaR.Because it is able to quantify losses more efficiently than the discrete distributions, and moreover it is coherent against of the VaR.In addition, optimization of the CVaR leads to optimize of the VaR and also linear programming approaches can be used for minimization of the CVaR.In the dynamic supply portfolio selection problem discussed in this section α is assumed to be fixed by DM to control the risk of losses due to supply delays and disruptions.The DM is willing to accept only portfolios which the total probability of scenarios with costs greater than VaR is not greater than 1-α.Furthermore, the DM wants to minimize the worst-case costs exceeding VaR.
Now we introduce the scenario-based calculation of the CVaR which is adapted to our work by following (Rockafellar & Uryasev, 2002;Uryasev, 2000).Assume that positive values of g(x,β) represent losses and β has a finite discrete distribution with L realizations and corresponding probabilities given as π l for βl, l=1,…, L (l is representative of an individual scenario) with π l ≥ 0 and l = 1 l L π = 1 å .For g(x,β), the α-CVaR can be stated by the following minimization formula: Where (g(x,β)-VaR) + = Max{g(x,β)-VaR,0}.
Let α-CVaR for loss random variable g(x,β) be hα(x).So, the α-CVaR equation can be restated as follows: The above model is non-linear programming problem.And so, by considering additional variables Tl for representing Max{g(x,β)-VaR,0} for all (l=1,…, L), this non-linear equation can be transformed into a linear programming model as follows (Uryasev, 2000): T l ≥ g(x,βl)-VaR; lєL T l ≥ 0; lєL From another point of view T l can be the amount by which costs (g(x,βl)) in scenario l exceed VaR.Thus by the above model the portfolio will be optimized by calculating VaR and minimizing CVaR simultaneously.In our model to select dynamic supply portfolio due to equation ( 9) for minimizing expected worst-case cost per part we have (corresponding constraints will be mentioned in the next part): Now we summarize equations and phrases in the former section in order to exhibit the complete form of our model.
Objective #1: Minimizing the net present value of expected cost per part (the first objective is like equation ( 6) with additional coefficients to consider time value of money for each phrase): Where: Objective #2: Selection of dynamic supply portfolio to minimize expected worst-case costs: F2 (Var, T) = Equation (10). Constraints: 1. Order assignment constraints: All the demands for each part type in each period are supplied by exactly one supplier in one discount interval.
For each supplier the total quantity of ordered parts cannot exceed its capacity.
d j x ijt ≤ c ijt y ijt ; i∈ I , j∈ J, t∈ T (13) For each part type the total dedications to suppliers must be equal to 1 and corresponding deliveries must not be delivered earlier than the earliest and not later than the latest delivery date.

Portfolio selection constraints:
The portfolio definition constraint is: The supplier i is selected for delivery of parts in period t if at least one order for part type j is assigned to supplier i in period t.
The average defect rate of the portfolio cannot be greater than q.
The average disruption rate of the portfolio cannot be greater than r.
The average delay rate of the portfolio cannot be greater than s.

Risk constraints:
The tail cost of scenario l for all parts which is supplied from all suppliers is defined by the nonnegative amount which cost per part in scenario l exceeds VaR, T l ≥ Right hand side of equation ( 11) -VaR; I є L (24) It is worthy to say that the variable VaR does not need to be constrained of being nonnegative.Since Tl is constrained of being non-negative, the model tries to decrease VaR and have a positive impact on the objective function.However, large reduction in VaR may result in more scenarios with costs more than VaR (Sawik, 2011a).
Discount constraints are followed by Sawik (Uryasev, 2000) with modifications: For a business volume discount environment, the total business volume purchased from each selected supplier can be exactly in one discount interval.
(where Qijtw = xijt zijtw) or For a quantity discount environment, the total quantity purchased from each selected supplier can be exactly in one discount interval.Constraints for making the second term of equations ( 11), ( 26) and ( 27) are linear.
We define Qijtw = Zijtw.xijt as an auxiliary variable.Then we need to add the following constraints to our model.

Solution methods
After developing the bi-objective mixed integer program for the dynamic supplier selection portfolio in a discount environment, now we have to propose an appropriate method for solving it.Reviewing literature for the procedure of solving same models shows that most of them solved them in small to medium scales test problems by AMPL programming language and the CPLEX solver (Sawik, 2010;Sawik, 2011c;Sawik, 2011a;Sawik, 2011b;Sawik, 2012).So as our model is multi-objective, we must search methods for solving multi objective models.
For solving small to medium size of our model, we suggest the Reservation Level driven Tchebycheff Procedure (RLTP) that could be coded in common software packages like Lingo or GAMS.This method is announced as a strong method for generating non-dominated solutions by literature for multi objective mathematical models.RLTP is suitable for the proposed model in this paper; because an interactive procedure is made by the DM for determining the tradeoff between worst case costs (CVaR) and expected cost, hence the DM can select the best solution based on his/her risk preference through this method.In fact this approach uses reservation levels based upon DM responses for reducing of the objective space in each of its iterations.
Some other advantages of RLTP are that it produces only Pareto-optimal solutions and can solve problems with non-convex nature; also this method is suitable for mixed integer multiobjective models.
To the best of our knowledge, heuristics and meta-heuristics methods have developed rarely for the proposed type of supplier selection portfolio.As Sawik (2011a) and Sawik (2012) suggested developing heuristics approaches for solving large size problems, because the size of our model increases rapidly by the number of suppliers (M), part types (N), planning horizon (H), discount intervals (m) and specially scenarios (L).The total number of potential supplies in all periods is (MNHm) and we can define each scenario as a supply realization; hence the total number of potential scenarios (L) as well as the number of variables and constraints in the mixed integer programs grow exponentially in terms of (MNHm).Unfortunately for a greater number of suppliers, part types, planning horizon, discount intervals and scenarios the model could not be solved by GAMS in a reasonable time.As a result, for large scale test problems we suggest a common and strong meta-heuristic method for solving multiple objective problems, NSGA-II.This method is so appropriate for solving multi objectives models and has been applied sufficiently in literature of multiple objective models.This method produces Pareto-optimal solutions and is appropriate for mixed integer multi-objective mathematical models.
It should be mentioned that the purpose of this paper is to develop a new model and suggest solution algorithms for the small, medium and large size problems.Undoubtedly other exact and meta-heuristic algorithms may have better results and this requires more discuss that is out of the scope of this paper and could be investigated in the future researches.In what follows, we describe a little about our selected methods i.e.RLTP and NSGA-II.Interested readers can refer to references that will be mentioned.

The exact solution procedure
The methods of multi-objective optimization have been successfully applied in the literature.
Tchebycheff metric based approaches have become popular for sampling the set of nondominated solutions in an interactive search for the most preferred solution in multiple objective decision making situations.These approaches systematically reduce the set of nondominated solutions which remain available for identification and selection from one iteration to the next (Reeves & MacLeod, 1999).Hereafter, we will prepare a little definition about the Reservation Levels Tchebycheff metric Procedure.
Suppose we have K objective functions, in RLTP method first we have to specify the number of solutions to be presented to the DM at each iteration, (P ≥ K), also we need to compute a reference vector, z**, where zi**= max {fi(x) / x є S}+εi and the εi is small positive scalar for each objective.Set Rli=-∞, (i=1,…, K) where RLi is the reservation level for the ith objectives.
Second it is necessary to generate a group of 2P dispersed weight vectors, A={λ є R k / λiє(0, 1), ∑λi=1.Third, the associated Tchebycheff program for each weight vector must be solved (the following model): Min {α-ρ∑zi} (36) Where ρ is a small positive scalar and the suggested range for this parameter is from 0.0001 to 0.01.After solving the model ( 36), the P most different of the resulting objective vector must be presented to the DM.If the DM wishes to continue to search for an improved solution, proceed to step 4. Otherwise, have the DM selects his/her current most preferred solution and stop.Forth, DM must partition the current solutions into more and less preferred solution, and adjust the RLs.Now we have to return to step 2 and continuing continue the algorithm.
Additional explanation about RLTP can be seen in Reeves & MacLeod (1999).
Using the RLTP algorithm and considering our model constraints, we arrive at a singleobjective, mixed integer programming model in each iteration; which can be solved efficiently using a linear programming solver.

The proposed NSGA-II method
Optimization problems are involved in maximizing or minimizing of the objective value function of the problem in order to satisfy all sets of constraints.If we have only one objective, the problem is considered as single-objective and with the use of the most meta-heuristics algorithms such as tabu search (TS), differential evolution algorithm (DEA), genetic algorithm (GA), simulated annealing (SA), particle swarm optimization (PSO), etc. a near global optimum solution or even global optimum solution can be obtained (Coello, Lamont & Veldhuizen, 2007).But if it is desired to consider more than one objective simultaneously, the problem is considered as multi-objective optimization problem and former algorithms are not proper anymore because they only find one solution in each running, therefore, new algorithms must be designed for these type of problems and the non-dominated sorting genetic algorithm (NSGA) was one of the first proposed algorithms by Srinivas and Deb (1994) for solving multiobjective optimization problems.Highlights of this algorithm are presented as follows: • Solutions are sorted based on how many solutions are better than the current solution.
• Fitness for solution is assigned based on their rank and non-dominating of other solutions.
• Using fitness sharing method for near solutions in order to solutions' variations are set properly and solutions are distributed uniformly in search set.
But, there are three main criticisms in the implementation of NSGA: first, high computational complexity of non-dominated sorting.Second, lack of elitism, and finally, need for specifying the sharing parameter (Deb, Pratap, Agarwal & Meyarivan, 2002).As for these shortcomings, second version of NSGA, non-dominated sorting genetic algorithm II (NSGA-II), was proposed by Deb et al. (2002), which is based on Pareto domination.
Highlights of this algorithm are presented as follows: • Using binary tournament selection operator.
• Defining crowding distance as alternative property for fitness sharing.
• Creating and saving non-dominated solutions which are obtained from pervious algorithm steps.
In NSGA-II the same as all evolutionary algorithms (EV), the first step is to create initial population.A good population size is related to the problem, at one side by increasing the population size, the chance of the algorithm to find a better solution will be improved, but at the other side the computational time of the algorithm for searching solutions increases exponentially, in this paper, we consider 100 as the population size.At the second step the fitness functions of all solutions are evaluated which at this paper these values are calculated by equations ( 10) and ( 11).If there is only one fitness function we can easily sort solutions but when there is more than one fitness function, the solutions should be ranked and sorted based on all fitness functions, therefore, solutions are divided into two main categories, dominated and non-dominated solutions.In a minimization problem, consider two solutions, m1 and m2, the solution m1 dominates the solution m2 if and only if the following conditions are satisfied.
 i є {1, 2,..., n}: fi(m1) ≤ fi(m2) (a)  j є {1, 2,..., n}: fj(m1) ≤ fj(m2) (b) If both of the above conditions are satisfied, then the solution m1 dominates the solution m2, and m1 is a non-domination solution, but if each one of the conditions (a) or (b) is violated, the solution m1 does not dominate m2, and based on these conditions, the NSGA-II ranks solutions.At the fourth step, some solutions of each generation are selected by the binary tournament selection method.Binary tournament selection method selects randomly two solutions through population, and then compares them and finally the better is selected.
Selection criteria in NSGA-II primarily are based on solution rank and after that are based on crowding distance.It is more favorable whatever solution rank is lower and crowding distance is higher.Now, with repeating binary selection operator on population of each generation, a set of population is selected for participating in crossover and mutation.In this paper, the crossover probability (cp) and the mutation probability (mp) consider as 0.8 and 0.3, respectively.
Accordingly, based on (cp) each pair of chromosomes are selected and with the implementation of the crossover operation new offspring are created, furthermore, based on (mp) the mutation operation is employed for each chromosome to explore more solution space.The created population merges with main population.The new population is sorted in ascending order by their ranks.Members of population with same rank are sorted in descending order by crowding distance.As the same size as main population, members are selected from new population.
The selected members create next generation population, and this cycle continues until termination conditions are fulfilled.Figure 1 and Figure 2, show schematic illustration of crowding distance criterion and pseudo-code of the NSGA-II, respectively (more information about these pictures are available in their corresponding references which are noted below them).

Numerical experiments
Now, numerical examples are presented to demonstrate possible applications of the proposed approach for supplier selection portfolio in a discount environment.As we have suggested two solution methods, in this part a test is regarded for each of them.Small size test problems are solved in a reasonable computation time by GAMS.But as the dimension of the test problems grow, computation time is amplified exponentially in GAMS.Hence large size test problems require the meta-heuristic method to be solved.

Defining parameters
Our small size example includes 3 suppliers; 3 part types; 3 periods; and 3 scenario deliveries.
Also we only consider quantity discount constraints in this example with 3 discount intervals for all suppliers and all part types in all periods.For delivery scenarios of supplier i, part type j and period t we assumed that the delivery scheduled for that period occurs either in period t or in period H or not all; i.e. ordered supplies may occur on-time (Δ l ijt= 0; for all i, j, t) or with the longest delay (Δ l ijt = H-t; for all i, j, t) or never (Δ l ijt = H-t+1; for all i, j, t).The probabilities of different scenarios (πijtl) are shown in Table 2.It should be noted that we assumed that the probability of occurrence a special scenario in a particular period is independent from earlier periods and their scenarios.Also the source of random generation for other parameters is shown in Table 3, Table 4, Table 5 and Table 6.

Findings through RLTP
The proposed model was solved through coding in GAMS 23.5 and RLTP was implemented for test problems.We developed a mixed integer model (MIP) so the CPLEX solver on a laptop with Intel Core i7 processor running at 2 GHz and with 4GB RAM (DDR3) was used.In what follows, first, we describe the results through the steps of RLTP and after that the final results will be illustrated.

Implementing RLTP for the proposed model
As declared in section 5, to implement RLTP for the mentioned example, following steps are done: Iteration 1: As we have 2 (K=2) objectives so we consider presenting 4 (n=4) solutions to the DM at each iteration.Also the best solution for each individual objective is obtained by solving two SOP (Single Objective Programming) model with corresponding constraints, the results are: Reservation levels (RLs) for the first iteration will be considered as follows: RL1 1 = 32 (First objective reservation level in iteration 1) RL2 Here we solve the associated AWTP ("Augmented weight Tchebycheff procedure", for more explanation see Reeves & MacLeod, 1999) model for each λ and present the n (=4) most different (maximally dispersed) of the resulting objective vector to the DM.The results are presented in Table 8.Now if the DM is satisfied with the best values of each objective, the RLTP will stop; otherwise it is necessary to go the next step that means the next iteration will start.In the proposed example we suppose the DM is not satisfied with the presented values so next iteration will be as follows.
Iteration 2: Based on the opinion of the DM the results are divided in two clusters, the most preferred solutions and the least preferred solutions.This task is shown in Table 9.  23.486, 23.475 23.525, 23.468 Table 9.The most and the least preferred of objectives The RLs have to be adjusted and the following procedure will be used: For each objective in each iteration, the difference between the worst value of the current most preferred solution for that objective and the worst value for that objective over all current solutions was determined.A percentage of that difference, 50% in this example, was then subtracted from the current most preferred solution value to arrive at the RL for that objective for the next iteration.To illustrate, the worst value of the first objective in the initial most preferred solution was 23.568.The worst initial value for that objective over the set of all current solution was 23.693 for a difference of 0.125.Subtracting 50% of that difference from the worst solution value over the set of all current solution yields a RL for the second iteration of 23.6305.This procedure for the second objective will produce 23.475.Using a higher or lower percentage would decrease or increase the rate of objective space reduction, respectively.The results of this part are: RL1 2 = 23.6305(First objective reservation level in iteration 2) RL2 2 = 23.475(Second objective reservation level in iteration 2) After this adjustment, we go back again to step 2 of RLTP.
The generated vectors for second iteration are shown in Table 10.so by dedicating larger weights to the first objective we consider the DM as a risk seeker one while by dedicating larger weights for the second objective means considering the DM as a risk avoiding one.Figure 1 shows the Pareto front which was made through two steps of RLTP.In other words the tradeoff between the expected cost and the expected worst-case cost is clearly shown in Figure 1, where the convex efficient frontier of Mean Cost-CVaR model for the confidence level α =0.75 is presented.As mentioned in (Reeves & MacLeod, 1999), RLTP only produces non-dominated solutions and this claim is illustrated in Figure 3.In this section the details of the selected solution through RLTP (i.e.23.568 for objective#1 and 23.461 for objective#2) is intended to be described.The model is run and the results are given as follows.
Since the main goal of the model is to determine the portfolio's variables, we go straight ahead to this part of solution first.The optimal supply portfolios (fraction of total demand dispersed between suppliers in the planning periods) are shown in Figure 6.In Figure 3   It is worthy to make a sensitive analysis on confidence level of the model.In the example which was solved before, α was set at 75%, which means that focus is on minimizing the highest 25% of all outcomes scenarios.Now we solve again the example by α=20 %.The resulted portfolio is shown in Figure 8.By comparing Figure 6 and Figure 8 we can see when α increases, a more risk-averse DM increases the number of selected supplies to mitigate the impact of risks by diversification of the supply portfolio (the number of selected supplies is 3 for α=25 % and 4 for α=75 %).
Figure 6.Supply portfolio for α=25% Figure 7. Supply portfolio for α=75% with an increase in delay and disruption penalizes Another important sensitive analysis can be done on penalty cost of delays and disruptions for the last solution produced by RLTP; we increased these penalties as Table 12 and results for portfolio variables are depicted in Figure 9.It can be concluded form new results that the more penalty values for disrupted or delayed supplies, the more balanced portfolio will be.This idea is concluded from the caparison between Figures 6 and 9 where it shows by more penalty values DM makes to dedicate more leveled orders from different suppliers to disperse more risks.In fact by more leveled portfolio, if an individual supply is disrupted or delayed, and if that supply has a high value, a large penalty will be charged.While with more leveled portfolio maximum risk will be minimized.As it is obvious in Figure 9, the smaller are differences between the highest and the lowest fraction of allotted demand.19 and Table 20; also Figure 9, depicts the results of portfolio variables which are the most important results of the model.The same sensitive analysis like section 6.1.4can be developed for large scales, but for avoidance of repetition they are not motioned again here.

Conclusion and future research
In this paper we presented a scenario-based multiple objectives model for dynamic supplier selection portfolio problem.We demonstrated that the proposed dynamic portfolio approach and the mixed integer programming model with conditional value-at-risk as a risk measure can be applied for the selection of suppliers and the allocation of orders quantity in supply chains with risks and discount constraints.As we considered dynamic supply portfolio, not only orders are allotted between suppliers but also are dispersed in planning horizon; hence the risks are mitigated and dispersed much more than the static supply portfolio.Our multi objective model minimizes the net present value of expected total cost over all scenarios; and at the same time the model aims to minimize the mean outcomes of scenarios which their costs exceeds VaR.
Two solution methods proposed in this paper for solving the model, the first one, which was coded in GAMS and solved through CPLEX solver, is RLTP.This method solved small to medium size test problems.The other method is NSGA-II and solved large size test problems.The nondominated supply portfolios, which were made, emphasized the effect of varying cost/risk preference of the DM.To demonstrate the effectiveness of the proposed model, various test problems were developed and results were reported.
Sensitive analysis was done and limited computational experiments indicated that the most important results from sensitive analysis was that the disruption and delay rates of suppliers are the most effective parameters in finding the optimal portfolio.The suppliers associated with the lowest disruption and delay rates are allotted the largest fraction of total demand.On the other hand, experiments demonstrated that the more risk averse the supply portfolio is, i.e. the higher confidence level (α), the more supplies are selected over planning horizon.
Furthermore, the number of selected suppliers increases with the confidence level, which demonstrates that the impacts of risks are mitigated by diversification of the supply portfolios.
At the end, we make the following relevant suggestions for future research: (1) determining the occurrence probability of scenarios in a more precise way because in the proposed model deliveries beyond the planning horizon are not allowed.It is possible to consider also the other scenarios that allow for delayed deliveries beyond the planning horizon.
(2) embedding scenarios to demand of part types (3) formulating the model as a multi stage stochastic one (4) determining the rates of disruption, delay and defect in a more realistic way (5) it would be interesting to develop a new approach based, for example, on fuzzy logic, to deal with uncertainty of disruption risks (6) our multi objective model in small sizes can be solved by other approaches, like multi-choice goal programming with utility method for generating efficient solutions, and the results can be compared to RLTP (7) also other meta-heuristic methods can be developed for large scale test problems.

Figure 3 .
Figure 3. Pareto front which was made through two iterations of RLTP average of the defect rate, delay rate disruptions rate for all part types for each period and each supplier are shown.Comparison of Figure 6 and Figure 7 indicate the lowers are defect, disruption and delays rates, the greater is the fraction of total demand.It is also observed that the supplies (i,t) associated with the lowest average disruption and delay rate (1,3) are allotted the largest fraction of total demand.The results indicate that the average disruption and defect rates are key determinants in the decision of dynamic allocation of demand among the suppliers.

Figure 4 .Figure 5 .
Figure 4. Supply portfolio for α=75% Now in Figure 8, we can comprise RLTP and NSGA-II methods.As shown in this figure NSGA-II has made more dispersed and a lot better solutions than NSGA-II.Also RLTP has a lot of difficulty in finding the best solution of each objective individually as well as in finding the solution of RLTP objective function.Details of one solution through RLTP for instance are shown in Table

Figure 8 .
Figure 8. Solutions which are made through NSGA II and RLTP for the large scale test problem Wijt is the set of discount intervals of supplier i for part type j in period tBjPer unit and per period cost of unfulfilled order for part type j, caused by supply disruptions for part type j Capacity of supplier i for part type j in period t dj Demand for part type j ej The earliest delivery date of part type j fj The latest delivery date of part type j oijt Cost of ordering part type j from supplier i in period t pijt Unit price of part type j purchased from supplier i in period t qijt Expected defect rate of supplier i for part type j in period t rijt Expected disruption rate of supplier i for part type j in period t sijt Expected delay rate of supplier i for part type j in period t q The largest acceptable average defect rate of supplies (defined by DM) r The largest acceptable average disruption rate of supplies (defined by DM) s The largest acceptable average delay rate of supplies (defined by DM) α Confidence level (defined by DM)

Table 1 .
Notations -225-Journal of Industrial Engineering and Management -http://dx.doi.org/10.3926/jiem.880 Parts are ordered from supplier i, if supplier i is selected and assigned at least one part order.

Table 3 .
Values of all parameters

Table 4 .
Values of aj, Bj and dj

Table 5 .
Values of unit price of part type j purchased from supplier i in period t

Table 6 .
Values of expected defect rate, disruption rate and delay rate of supplier i for part type j in period t

Table 7 .
1 = 32 (Second objective reservation level in iteration 1) Now we generate a group of 8 (=2×n) dispersed weight vectors randomly to use them in the next step.The generated vectors are shown in Table 7: Weight vectors for each objective

Table 8 .
Objective values for selected λ

Table 10 .
Weight vectors for each objectiveSolving the associated model for each λ and presenting the 4 most different (maximally dispersed) of the resulting objective vector to the DM, is done in Table11.

Table 11 .
Objective values for selected λFinally in this example we can suppose that the DM is satisfied by solution that is produced by λ6, i.e. 23.568 for objective#1 and 23.461 for objective#2 as the most preferred solution and RLTP is stopped.It is important to emphasis that the RLTP procedure is an interactive way for finding Pareto front.Also RLTP makes it possible to embed the risk preference of the DM by considering different weights for each objective i.e. as the first objective minimizes the expected net present value of costs and the second objective minimizes the worst case costs,

6.2. Large size test problem with NSGA-II In
this part we develop a test problem for large size problems that includes 5 suppliers; 10 part types; 12 periods; 3 delivery scenarios and 5 discount intervals.This size of model can be considered as a real size that can be applicable in real word.To demonstrate the capability of our meta-heuristic algorithm first we tried to make some non-dominated solutions by RLTP and then by our NSGA-II method for this test problem.Hence, we can compare solutions and see how good NSGA-II method is.It should be noted that a lot of test problems were made by the provided code in MATLAB (R 2010a version) (MATLAB, 2010), but for a lot of test problems RLTP method which is coded by GAMS could not make any solutions after a long time, so the comparison with our NSGA-II was not possible.Parameters of the proposed model are shown in Table15.As we have 2 (K=2) objectives so we will consider to present 8 (n=8) solutions to the DM at each iteration.So we randomly generate a group of 16 (=2×n) dispersed weight vectors to use them in the next step.The generated vectors are shown in Table16.Table16and RLTP, are stopped after one iteration.It should be noted that the solutions for each λ had considerable gap for RLTP objective function that those gaps are shown in the last column of Table17.

Table 17 .
Objective values for selected λ -247-Journal of Industrial Engineering and Management -http://dx.doi.org/10.3926/jiem.880NSGA-II method made the following results for each objective.Results of iterations 57, 100 and 350 are illustrated in Table 18.

Table 18 .
Objective values for selected λ

Table 19 .
Objective's values and confidence level of NSGA-II for the proposed test problem

Table 20 .
Results of portfolio variables (Fit) resulted from NSGA-II Figure 9. Values of supply portfolio variables for α=75% in a large size test problem