A multi-period multi-commodity lot-sizing problem with supplier selection, storage selection and discounts for the process industry

European Journal of Operational Research 279 (2019) 393–406

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Production, Manufacturing, Transportation and Logistics

A multi-period multi-commodity lot-sizing problem with supplier selection, storage selection and discounts for the process industry Thomas Kirschstein a,∗, Frank Meisel b a b

School of Economics and Business, Martin Luther University Halle-Wittenberg, Gr. Steinstr. 73, Halle 06108, Germany Department of Business Administration, Christian-Albrechts-University zu Kiel, Wilhelm-Seelig-Platz 1, Kiel 24098, Germany

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 20 February 2018 Accepted 29 May 2019 Available online 13 June 2019

In this paper we consider the sourcing process of a production facility from the process industry. The production facility sources raw materials from a set of suppliers that offer different discount schemes. For storing raw materials, heterogeneous storage facilities can be used, which vary with respect to their associated costs and capacities. We formulate the corresponding planning problem of selecting suppliers and storage facilities as well as determining order quantities and transport flows under the discount schemes offered by the suppliers. For solving large problem instances, a heuristic based on kernel search is proposed and evaluated in a real world case study. The case study reveals that supplier selection and storage selection are highly interdependent decisions. By integrating both perspectives, significant savings can be generated. It turns out that single-sourcing is the dominant strategy for most raw materials. However, the selection of the optimal suppliers depends on both the prices offered by the suppliers as well as the associated logistics costs for transportation and stock-holding.

Keywords: Multi-commodity lot-sizing Discounts Supplier selection Storage selection Kernel search heuristic

© 2019 Elsevier B.V. All rights reserved.

1. Introduction In recent years, sourcing and production processes have become more and more complex and globalized in many industries. As a consequence, supply chains involve an increasing number of suppliers and logistics service providers, which are interdependent regarding material, information, and financial flows. Since procurement is a fundamental requirement for any production activity, the integration and coordination of suppliers and logistics service providers is a core issue for persistent economic success of producing firms. This incurs the need to integrate the capabilities and capacities of the logistics system and the suppliers in order to benefit from an overall improvement regarding procurement cost, logistics cost as well as supply reliability and flexibility. In this paper, we introduce a multi-period, multi-commodity lot-sizing problem with supplier selection, storage facility selection, and quantity discounts. Thereby, a single production facility producing multiple final products demands a set of raw materials, which are to be sourced from multiple suppliers. Such production systems most prominently occur in the process industry where a portfolio of final products is produced by multi-purpose machinery

∗

Corresponding author. E-mail addresses: [email protected] (T. Kirschstein).

[email protected],

https://doi.org/10.1016/j.ejor.2019.05.039 0377-2217/© 2019 Elsevier B.V. All rights reserved.

using numerous raw materials and various recipes. Typically, each supplier offers one or more raw materials whereby the charged unit prices follow discount schemes. We consider here three different discount schemes for the suppliers. Each supplier applies either a single quantity discount (SQD), a total quantity discount (TQD), or a hybrid quantity discount (HQD). Single quantity discounts are granted separately for each material depending on its order quantity. Total quantity discounts are calculated based on the total order quantity over all materials which are part of an order. An HQD combines single and total quantity discounts, i.e., the unit price for a certain raw material which is part of a multi-material order depends on the material’s order quantity as well as on the total order quantity of all ordered materials. For storing the raw materials, different storage facilities are available at the production site but also in the neighborhood of the production asset. Furthermore, materials can be transferred between storage facilities if required. Each storage facility consists of a limited number of (technically identical) storage items. Such storage items are, for example, silos or containers. The resulting lot-sizing problem is N P -complete. We model this problem as a mixed-integer linear program (MIP) and propose a tailor-made kernel search heuristic to solve large problem instances for a time horizon of up to six months. A computational study based on real-world instances shows the applicability and the performance of the proposed solution methods. The outline of this paper is as follows. In Section 2, we briefly review relevant literature on supplier selection with quantity

394

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

discounts. Afterwards, in Section 3, the planning problem is described in detail and a mathematical model is formulated. In Section 4 we describe the heuristic solution method. In Section 5, we present a real-world inspired case study from the chemical industry to highlight the applicability of the proposed model and the heuristic. The paper is concluded in Section 6. 2. Literature review Lot-sizing problems with supplier selection and discounts have found fertile ground in the scientific literature for a long time. Particularly, lot-sizing problems with quantity discounts have a long tradition, see e.g. Munson and Rosenblatt (1998) for an early overview about quantity discount policies. Basically, two types of discount policies are distinguished: all-unit discounts (i.e., discounts that apply for all units of an order) and incremental discounts (i.e., discounts that apply only for those units that exceed the threshold of the corresponding discount interval). All-unit discount models and incremental discount models have been investigated since the seminal works of Hadley and Whitin (1963) and Low and Waddington (1967). Broad literature reviews on discount problems can be found in Benton and Park (1996), Munson and Rosenblatt (1998), and Munson and Jackson (2015). For the case of one commodity with deterministic demand over a finite number of discrete periods and two discount intervals, Federgruen and Lee (1990) proposed optimal dynamic programming algorithms which were later improved by Xu and Lu (1998). Chaudhry, Forst, and Zydiak (1993) extended the work of Federgruen and Lee (1990) by considering multiple suppliers with either incremental or all-unit discount schemes in a single-period problem. A literature overview on joint lot-sizing and supplier selection problems can be found in Aissaoui, Haouari, and Hassini (2007). Tempelmeier (2002) designed a powerful heuristic for the case of multiple suppliers employing either all-unit or incremental discount regimes with multiple intervals. Additionally, their model includes time-varying purchasing prices as well as delivery periods and lead times. Similar heuristics are proposed by Parsa, Khiav, Mazdeh, and Mehrani (2013) and Mazdeh, Emadikhiav, and Parsa (2015). Stadtler (2007) proposed a generic MIP model for multiple suppliers and periods where the discount type is defined by the input data. The proposed MIP is strengthened by multiple valid inequalities. Lee, Kang, Lai, and Hong (2013) formulated a MIP for multi-period, multi-supplier settings where all-unit and incremental discounts are considered jointly. The problem is solved by a genetic algorithm. Recently, Ghaniabadi and Mazinani (2017) extended the model by allowing backlogs. For the case of incremental discounts, a dynamic programming approach is proposed to optimally solve the problem whereas for the case of all-unit discounts, commercial MIP solvers are used. In case of multiple commodities, two sub-schemes of all-unit discount models can be distinguished: Single-quantity discounts (SQD) and total quantity discounts (TQD). In the SQD-scheme, unit prices charged for a particular commodity depend on the order quantity of this very commodity only, i.e., discounts are independent for the ordered commodities. In contrast, in the TQD-scheme, discounts apply for all commodities involved in an order depending on the total quantity ordered. The single-period TQD problem has been shown to be N P -complete (Goossens, Maas, Spieksma, & Van de Klundert, 2007). Solution approaches range from tailormade branch-&-cut algorithms (Goossens et al., 2007; Manerba & Mansini, 2012) to MIP-based heuristics (Manerba & Mansini, 2014). Mansini, Savelsbergh, and Tocchella (2012) extended the problem by incorporating transportation cost and capacities of the suppliers. The model is solved by an iterative rounding heuristic. This TQD model is structurally very similar to the multi-period all-unit discount models studied in Lee et al. (2013) and the so-called

business-volume discounts studied in Balakrishnan and Natarajan (2014). In all model variants the fundamental similarity is a knapsack-like structure of the lot-sizing decisions, which creates cost benefits when a sufficient quantity is ordered from a supplier. This knapsack property is the root cause for the complexity of the problems and the reason why tailor-made branch-&-cut procedures work well for these problems. Crama and Torres (2004) presented a case study from the chemical industry by proposing a tactical multi-site production planning problem with total quantity discounts and recipe selection where neither transportation nor inventory holding is considered. Instances with two plants, eight suppliers, 25 ingredients, and 104 recipes are solved using a standard solver. Manerba, Mansini, and Perboli (2018) extended the capacitated single-period TQD problem by considering fixed cost for choosing a supplier (so-called activation cost) together with uncertainty of demands and unit prices. The proposed problems are solved by an exact solution approach that is inspired by Manerba and Mansini (2012). Manerba and Perboli (2019) extended this work by proposing a heuristic and a Benders approach for solving the capacitated single-period TQD problem with demand uncertainty. A similar study on incremental and total quantity discounts with stochastic supplier failures can be found in Bohner and Minner (2017). The effects of quantity discounts on limited and jointly used resources has been studied recently by Jackson and Munson (2016) and Shahsavar, Zoraghi, and Abbasi (2018). Jackson and Munson (2016) investigated a static multi-material lot-sizing problem with a TQD or an SQD scheme under a shared and limited resource that represents, for example, a storage space. Thereby, the resource capacity can be extended. The proposed model minimizes the sum of purchasing cost, holding cost, and capacity extension cost whereby the latter is modelled as a linearized capacity cost function to support maximum flexibility. The problems are solved by iteratively testing segments of the discount and capacity cost functions. Shahsavar et al. (2018) propose a resource-constrained project scheduling problem where each project is associated with material and resource demands. A TQD scheme is applied for purchasing the materials. Various procedures based on a genetic algorithm are proposed for solving the problem. This paper contributes to the literature by proposing two new problem aspects: First, different storage facilities with limited capacity have to be used for intermediate storage in a multi-period lot-sizing problem with discounts. The storage facilities differ regarding cost structures, capacities of the storage items, and transportation requirements. Second, we model a heterogeneous set of suppliers regarding the offered discount schemes. This means that each supplier offers either an SQD, or a TQD, or a hybrid quantity discount (HQD). Additionally, for solving large problem instances, we propose a MIP-based heuristic that uses a kernel search. In contrast to existing heuristics for lot-sizing problems with discounts, the proposed heuristic also decides on storage and transportation while still exploiting structural information about promising material-supplier combinations that is deduced from the linear programming relaxation. 3. Dynamic multi-commodity lot-sizing problem with discounts, supplier selection and storage selection 3.1. Problem description We consider a procurement system with multiple suppliers and storage facilities that supplies a multi-purpose production asset with the necessary raw materials. Such systems are typical for the chemical industry. A schematic illustration of the considered system is shown in Fig. 1. The notation used for formally describing the planning problem is summarized in Table 1. The overall task is to cost-efficiently supply the production asset with raw materials

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

Fig. 1. Structural overview of the procurement and storage system.

Table 1 Notation of sets and parameters. Sets: M Ms T S Sm S SQD/TQD/HQD SQD/TQD/HQD Sm Em I Im Bsm Ds Parameters: amt y¯ i l¯i lmi0 x¯i j oij qsmb q¯ smb qsd q¯ sd cisetup inv−var cm ciinv− f ix − f ix citrans j −var citrans j price csmb price csmd price csmbd

Set of all considered materials Set of materials supplied by supplier s Set of time periods Set of all suppliers Set of suppliers that can deliver material m Set of suppliers with SQD, TQD, and HQD schemes Set of suppliers with SQD, TQD, and HQD schemes that can deliver material m Set of transport relations Em ⊆ (S × I ) ∪ (I × I ) that are relevant for material m Set of all storage facilities Set of storage facilities usable for material m Set of SQD-discount intervals of supplier s ∈ S SQD ∪ S HQD for material m Set of TQD-discount intervals of supplier s ∈ S TQD ∪ S HQD Demand of material m in period t (measured in tons, like all further quantities) Number of available storage items of facility i Inventory capacity of a storage item of facility i Initial stocks of material m in storage facility i Transport capacity per transport vehicle on relation (i, j) Transportation lead time for transport relation (i, j) (measured in periods) Minimum order quantity of SQD-interval b ∈ Bsm of supplier s for material m Maximum order quantity of SQD-interval b of supplier s for material m Minimum order quantity of TQD-interval d ∈ Ds of supplier s Maximum order quantity of TQD-interval d of supplier s Setup cost per storage item of facility i Stock-holding cost rate per period and ton of material m Leasing cost per period and storage item of facility i Transport cost per transport vehicle on relation (i, j) Transport cost per ton shipped on relation (i, j) Purchasing price per ton in for material m Purchasing price per ton in for material m Purchasing price per ton in TQD-interval d of hybrid

SQD-interval b of supplier s TQD-interval d of supplier s SQD-interval b and supplier s for material m

395

such that a given production plan can be executed on time. The given production plan implies demands amt for raw materials m ∈ M and periods t ∈ T of a multi-period planning horizon. Like all further quantities, demands are expressed in tons in this paper. Each supplier s ∈ S can deliver a non-empty subset Ms ⊆ M of the raw materials. Vice versa, Sm ⊆ S denotes the non-empty set of suppliers for material m. Each supplier offers either an SQD, a TQD, or an HQD scheme. Correspondingly, the set of suppliers S is disjunctively partitioned into S SQD , S TQD , and S HQD . If a material is ordered from a supplier, it is shipped to a storage facility i ∈ I. For technical reasons, a material m can be stored in a non-empty subset Im ⊆ I of storage facilities only. Each storage facility i ∈ I consists of a limited number y¯ i of homogeneous storage items, like, for example, a number of identical silos or a number of container slots. Each storage item of facility i has a storage capacity of l¯i tons. A specific storage facility is the preproduction silo farm 0 ∈ I that is coupled directly to the production facility. It is important to know that all materials have to go through these pre-production silos in order to be fed into the production system. For the pre-production silos, we assume y¯ 0 < |M|, i.e., at least some of its silos have to be assigned to different materials during the time horizon. Holding stocks in different storage facilities incurs different costs. This is because some storage facilities are managed and owned by the company itself while others are offered by logistics service providers. In the first case, only setup and capital commitment costs have to be considered. In the second case, also leasing costs are charged for the time that the silos, containers, etc. are setup rented. More formally, ci denotes the setup costs for preparing

inv−var an item of storage facility i to store a particular material, cm expresses the stock-holding cost rate per ton and period for mateinv− f ix rial m, and ci expresses the leasing cost rate per period and item of storage facility i. For storage facilities owned by the proinv− f ix ducing firm, ci is set to zero. It is furthermore possible to relocate stored materials from one storage facility to another for reasons of capacity utilization or to shift inventory to freed storage facilities of lower cost. For transportation between suppliers and storage facilities as well as between storage facilities, different transport vehicles are used depending on the storage facilities involved (e.g. silo or container trucks). Furthermore, we denote by Em ⊆ (S × I ) ∪ (I × I ) the set of links (edges) that are available for transporting material m in the considered network. Each vehicle that uses a link (i, j ) ∈ Em has a transport capacity of x¯i j tons and requires oij time units for performing the transport on this link. Astrans− f ix

sociated transport costs ci j

are charged per chartered vehicle

−var citrans j

and costs are charged per transported ton. The purchase price charged by supplier s depends on the supplier’s discount-based pricing system. For suppliers providing SQDdiscounts, Bsm denotes the set of discount levels for supplier s and material m. If a certain discount interval b ∈ Bsm is chosen, the order quantity must fall into the interval [qsmb , q¯ smb ] and the correprice

sponding unit price csmb is charged. For suppliers providing TQDdiscounts, Ds denotes the set of discount levels of supplier s. For discount interval d ∈ Ds holds that the total quantity of all materials ordered from this supplier in a period must fall into the interprice val [qsd , q¯ sd ]. In this case, the unit price is csmd . In case of hybrid discounts, both restrictions have to hold simultaneously. Then, if an order of material m is placed in intervals b and d, the HQDprice supplier s charges a unit price csmbd . The task is then to determine order quantities and a storage plan (consisting of inventory levels and an assignment of materials to storage items) that minimize total cost over a multi-period time horizon. Thereby, the material quantities ordered from a supplier in a period can be split up and delivered to multiple storage facilities. The total cost consist of purchasing cost, transportation cost,

396

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406



Table 2 Notation of decision variables.

α smbt

Binary, 1 if material m is ordered from supplier s in SQD-interval b in period t Binary, 1 if an order is placed at supplier s in TQD-interval d in period t Continuous in [0, 1], share of amount q¯ sb of material m ordered in discount interval b from SQD-supplier s in period t Continuous in [0, 1], share of amount q¯ sd of material m ordered in discount interval d from TQD-supplier s in period t Continuous in [0, 1], share of amount q¯ sb of material m ordered in discount interval b from HQD-supplier s in period t Continuous, quantity of material m ∈ M shipped on link (i, j ) ∈ Em in period t ∈ T Integer, number of transport vehicles for material m on link (i, j) in period t Continuous, stock level of material m in storage items of facility i in period t Integer, number of storage items of facility i used for material m in period t Integer, number of storage items of facility i set up for material m in period t

β sdt zsmbt zsmdt zsmbdt xmijt nmijt lmit ymit smit

inventory holding cost, as well as of leasing and setup cost of the storage items.

xms jt =

(s, j )∈Em

∀t ∈ T , m ∈ M, s ∈ SmSQD , b ∈ Bsm

(5)

αsmbt ∈ {0, 1}

∀t ∈ T , m ∈ M, s ∈ SmSQD , b ∈ Bsm

(6)

If an order is placed in interval b ∈ Bsm , the corresponding binary variable α smbt takes value 1 and constraints (2) and (3) enforce the  order quantities zsmbt · q¯ smbt to be in the selected discount inb∈Bsm

terval [qsmb , q¯ smb ]. Otherwise, zsmbt is bound to be zero. Constraints (4) ensure that the order quantities are shipped to suitable storage facilities. Domain declarations are formulated in (5) and (6). For the suppliers offering total quantity discounts (S TQD ), purchasing cost function and ordering constraints are slightly different. Here, binary variable β sdt indicates whether an order in interval d ∈ Ds is placed in period t with supplier s or not. The purchasing cost function changes as one single discount interval is selected which applies to all ordered material units such that

3.2. Mathematical model

   

t∈T m∈M s∈S SQD b∈Bsm m

price csmb · zsmbt · q¯ smb

(1)

where the quantity of material m ordered from supplier s in pe riod t is given by b∈Bsm zsmbt · q¯ smb . The corresponding constraints SQD associated with ordering from suppliers Sm are as follows

zsmbt ≤ αsmbt

zsmbt ≥ αsmbt ·

qsmb q¯ smb

   

price csmd · zsmdt · q¯ sd .

(2)

∀t ∈ T , m ∈ M, s ∈ SmSQD , b ∈ Bsm

(3)

(7)

The ordering constraints for the TQD suppliers are as follows



zsmdt ≤ βsdt

∀t ∈ T , s ∈ S TQD , d ∈ Ds

(8)

m∈Ms



zsmdt ≥ βsdt ·

m∈Ms



xms jt =

(s, j )∈Em



qsd

∀t ∈ T , s ∈ S TQD , d ∈ Ds

q¯ sd

zsmdt · q¯ sd

∀t ∈ T , m ∈ M, s ∈ SmTQD

(9)

(10)

d∈Ds

zsmdt ∈ [0, 1]

∀t ∈ T , m ∈ M, s ∈ SmTQD , d ∈ Ds

(11)

βsdt ∈ {0, 1}

∀t ∈ T , s ∈ S TQD , d ∈ Ds

(12)

Once the buyer has chosen a discount level in the total discount model, the buyer is free to distribute the order quantities for all materials as long as the total order quantity fits into the chosen discount interval’s limits, which is ensured by (8) and (9). Then, (10) to (12) balance the transport flows and declare the domains of the variables. Modeling ordering decisions for suppliers offering a hybrid discount scheme (S HQD ) requires to apply both lsingle discount and total discount decisions. Like in the SQD case, Bsm and Ds denote the sets of lsingle and total discount levels, respectively. Likewise, α smbt and β sdt are binary variables indicating the lsingle discount level and the total discount level chosen in period t with supplier s. The corresponding purchasing cost function for suppliers with hybrid discounts is

PC HQD (z ) =

    

t∈T m∈M s∈S HQD b∈Bsm d∈Ds m

∀t ∈ T , m ∈ M, s ∈ SmSQD , b ∈ Bsm

(4)

zsmbt ∈ [0, 1]

t∈T m∈M s∈S TQD d∈Ds m

We formulate here the outlined problem as a MIP. The relevant decision variables are summarized in Table 2 and described as follows. Binary decision variables α smbt take value 1 if material m is ordered from supplier s in SQD-interval b in period t, 0 otherwise. Similarly, βsdt = 1 if the total material quantity ordered from supplier s in period t falls into TQD-interval d, 0 otherwise. For SQDorders, continuous variable zsmbt ∈ [0, 1] indicates the share of the discount interval limit q¯ sb of material m ordered from SQD-supplier s in period t. zsmdt fulfills the same purpose for TQD-intervals. For HQD-suppliers, zsmbdt indicates the share of q¯ sb of material m ordered in single discount interval b and total discount interval d in period t. Transport decisions are formalized by the continuous variable xmijt , which refers to the quantity of material m ∈ M that is shipped on link (i, j ) ∈ Em in period t ∈ T and by integer variable nmijt , which is the number of transport vehicles for material m on link (i, j) in period t. The continuous variable lmit denotes the stock level of material m in period t in storage items of facility i. Integer variable ymit denotes the number of storage items used for material m at facility i in period t and variable smit denotes the number of storage items that are newly set up for material m at storage facility i in period t. The objective is to minimize total cost over the time horizon. For the sake of readability, we introduce auxiliary cost functions for purchasing cost, transportation cost, and inventory-related costs. We begin with the purchasing cost function, where we formulate an individual function for each discount scheme. In the SQD case, the purchasing costs are calculated by

∀t ∈ T , m ∈ M, s ∈ SmSQD

zsmbt · q¯ smb

b∈Bsm

PC TQD (z ) =

PC SQD (z ) =



price csmbd · zsmbdt · q¯ smb .

(13)

Note that the quantity of an order of material m placed in pe  riod t with supplier s is expressed as zsmbdt · q¯ smb . That is, b∈Bsm d∈Ds

zsmbdt is positive if the total order quantity over all materials is in interval d and the order quantity of material m is in interval b. As lsingle and total order quantities are linked, ordering constraints

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

for hybrid suppliers are a bit more complex. They are formulated as follows



zsmbdt ≤ αsmbt

∀t ∈ T , m ∈ M, s ∈ SmHQD , b∈ Bsm

(14)

∀t ∈ T , m ∈ M, s ∈ SmHQD , b∈ Bsm

(15)

d∈Ds



zsmbdt ≥ αsmbt ·

d∈Ds

qsmb q¯ smb

  zsmbdt · q¯ smb ≤ βsdt · q¯ sd

∀t ∈ T , s ∈

HQD Sm ,d

∈ Ds

(16)

∀t ∈ T , s ∈ SmHQD , d ∈ Ds

(17)

∀t ∈ T , m ∈ M, s ∈ SmHQD

(18)

m∈Ms b∈Bsm

  zsmbdt · q¯ smb ≥ βsdt · qsd

m∈Ms b∈Bsm



  xms jt = zsmbdt · q¯ smb

(s, j )∈Em

b∈Bsm d∈Ds

zsmbdt ∈ [0, 1]

∀t ∈ T , m ∈ M, s ∈ SmHQD , d ∈ Ds , b ∈ Bsm

(19)

αsmbt ∈ {0, 1}

∀t ∈ T , m ∈ M, s ∈ SmHQD , b ∈ Bsm

(20)

βsdt ∈ {0, 1}

∀t ∈ T , s ∈ SmHQD , d ∈ Ds

(21)

Constraints (14) and (15) determine the single order interval selection whereas (16) and (17) assure that the proper total discount   interval is chosen. Note that zsmbdt · q¯ smb expresses the m∈Ms b∈Bsm

total order quantity over all materials. Constraints (18) to (20) balance transport flows and define domains of variables. Furthermore, for an HQD-supplier, a combination of an single discount interval b and a total discount interval d is infeasible if qsmb > q¯ sd . To forbid such infeasible combinations, the following constraint can be added



∀m ∈ M, s ∈ SmHQD , b ∈ Bsm , d ∈ Ds : qsmb > q¯ sd .

zsmbdt = 0

t∈T

(22) While (1)–(21) model the ordering processes of the suppliers, the subsequent constraints (23)–(32) model the transport processes and material holding at the storage facilities.

lmit = lmi(t−1) +



xm ji(t−o ji ) −

( j,i )∈Em :o ji
lm0t = lm0(t−1) +



xmikt

(i,k )∈Em

 xm ji(t−o ji ) − amt

∀t ∈ T , m ∈ M, i ∈ Im \0 ∀t ∈ T , m ∈ M

(23)

(24)

( j,i )∈Em :o ji
 xm ji(t−o ji )

l¯i · ymit ≥ lmi(t−1) +

∀t ∈ T , m ∈ M, i ∈ Im

(25)

( j,i )∈Em :o ji


ymit ≤ y¯ i

∀t ∈ T , i ∈ I

(26)

∀t ∈ T , i ∈ I, t > 1

(27)

m∈M

smit ≥ ymit − ymi(t−1) xmi jt ≤ nmi jt · x¯i j smit , ymit ∈ N

∀t ∈ T , m ∈ M, (i, j ) ∈ Em ∀t ∈ T , m ∈ M, i ∈ Im

(28)

(29)

397

nmi jt ∈ N

∀t ∈ T , m ∈ M, (i, j ) ∈ Em

(30)

lmit ∈ R+

∀t ∈ T , m ∈ M, i ∈ Im

(31)

xmi jt ∈ R+

∀t ∈ T , m ∈ M, (i, j ) ∈ Em

(32)

Constraints (23) and (24) constitute inventory balances. The latter refers to the stocks of the pre-production silos and guarantees a fulfillment of demands amt . Initial stocks lmi0 at the beginning of the planning horizon are given as parameters. Constraints (25) determine the required number of storage items ymit for material m at facility i in period t. Note that the maximum stock level at the beginning of period t is used here instead of the stock level at the end of the period as defined by (23) and (24). This is because the maximum stock level determines the number of storage items to be set up or leased. With (26) the storage item limit for each storage facility is obeyed while (27) calculates the number of storage items that are newly set up for material m at facility i in period t. Constraints (28) calculate the required number of transport vehicles to realize a transport of xmijt tons of material m on link (i, j) in period t. Finally, (29) and (32) constitute the remaining domain declarations. Furthermore, cost functions for transportation and inventory holding are introduced as follows:



T C (n, x ) =

 

m∈M (i, j )∈Em t∈T

IC (y, l, s ) =

 

− f ix −var citrans · nmi jt + citrans · xmi jt j j



(33)

inv−var ciinv− f ix · ymit + cm · lmit + cisetup · smit



m∈M i∈Im t∈T

(34) The transportation cost function (33) summarizes transport costs trans− f ix ci j of vehicles used on transport relations (i, j) and the han−var dling costs citrans of transport quantities xmijt over all periods. j inv− f ix

Inventory costs (34) comprise leasing costs ci

charged per

inv−var cm

storage item and period, stock-holding cost per ton of masetup terial m in stock, and setup cost ci for storage items. Finally, the MIP minimizes the total cost

min → PC SQD (z ) + PC TQD (z ) + PC HQD (z ) + T C (n, x ) + IC (y, l, s ) subject to constraint sets (2)–(6), (8)–(12), (14)–(20), and (23)–(32). 4. Heuristic solution approach Basically, the model at hand can be categorized as a multiperiod lot-sizing problem with all-unit discounts as well as supplier and storage selection. As a single-period all-unit discount model is known to be N P -hard (Goossens et al., 2007), the same holds for the problem studied here. Nonetheless, literature on single-period all-unit discount models shows that such problems are fairly easy to solve by commercial solvers (Manerba & Mansini, 2012). One reason for this is the good-natured structure of these problems, which allows to convert all-unit discount models into network flow problems where the suppliers with their associated discount intervals form the source nodes of the network flows. This effects that sophisticated cut procedures for network flow problems can be applied, which are nowadays part of most commercial solvers. Furthermore, a solution of the LP relaxation of the network flow problem can be transferred into a feasible solution quite easily by rounding [see exact methods proposed by](Goossens et al., 2007; Manerba & Mansini, 2012). Heuristics proposed for allunit discount models, therefore, primarily rely on MIP-based approaches that round LP-solutions and restrict the solution space in

398

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

an elegant way (e.g. by preselecting promising suppliers for each material). The problem at hand, however, is more complex as also integer transport flows and capacitated storage facilities are considered, which incur fixed cost. Thus, standard lbranch-&-cut procedures are less effective in finding optimal solutions. Moreover, since the storage capacities are limited (particularly for the pre-production storage), finding a general feasible solution is more difficult, too. Therefore, we propose a MIP-based heuristic for solving the problem under investigation. The heuristic is composed of five steps and follows roughly the ideas of a kernel search (Angelelli, Mansini, & Speranza, 2007; 2010; 2012; Guastaroba, Savelsbergh, & Speranza, 2017): 1. The LP relaxation is solved to identify promising materialsupplier assignments. More precisely, a subset of the binary variables α and β is identified and called kernel. Those binary variables that are not in the kernel are grouped according to their associated reduced costs in a number of sets calledbrk buckets. 2. A partially relaxed MIP problem is defined by applying an LPrelaxation to the integer variables and by restricting the binary variables in all buckets to 0. 3. The obtained (partially relaxed) solution is then used as the initial solution for the original MIP. The initial solution is transferred into a general feasible solution by using the repair routine of Gurobi and fixing all 1-valued binaries from the initial solution. This way the set of free binaries is further reduced. 4. Subsequently, the buckets identified in the first step are iteratively added to the kernel in order to find better solutions. The procedure is stopped if the (relative) improvement obtained in an iteration is below a threshold  and the accepted number of non-improving iterations (nacc ) is exceeded. 5. Finally, the best found solution is improved by solving the original problem with binary variables α and β fixed according to the best solution found so far. To describe the heuristic more precisely, let S denote the solution of a problem instance consisting of the values for continuous (X ), integer (Y), and binary (Z) variables. Note that all integer variables (i.e., n, y, and s) are associated with transportation and storage processes while binary variables (α and β ) address ordering decisions. The total cost of a solution is expressed by C(S). Let RC(SLP ) denote the reduced costs of the solution of the LP relaxation SLP . Furthermore, Z>0 (SLP ) denotes the subset of binary variables that have a value strictly greater than zero in SLP . We will distinguish two types of subsets of binary variables: The kernel K and the buckets Fi where the kernel subsumes free binary variables used in the optimization process while the buckets contain binaries which are initially fixed to 0. The index i denotes here the ith out of n buckets. The kernel consists of all binaries that are associated to non-zero order quantities in the LP solution SLP . I.e., for any selected combination of a supplier, a period, and a material, all associated binaries are added to the kernel, which encompasses all corresponding discount intervals. The binary variables are grouped into buckets according to their potential to be part of the optimal solution. As a proxy for this potential, we use the reduced costs of the associated z-variables. Finally, let MIP (X , Y, Z |Z

= 0 ) denote the MIP solved with variable sets X , Y, and Z while variables Z

are fixed to zero. If integer or binary variables are relaxed, this is indicated by Y˜ and Z˜ , respectively. If an initial solution is provided   to a MIP model, this is indicated by MIP X , Y, Z |Z

= 0, S . The whole heuristic procedure is shown in Algorithm 1. Building the buckets of binary variables relies on analyzing the reduced cost of the continuous ordering variables z. Note that each binary variable is associated to one or more z-variables as formalized in the ordering constraints (2)–(3), (8)–(9), and (14)–(17). As

Algorithm 1: Adapted kernel heuristic.

1 2

3

4 5

6 7 8

9 10

Input : Problem instance Output: Feasible solution S obtain LP solution SLP = MIP (X , Y˜ , Z˜ ); kernel K ⊃ Z>0 (SLP ) consists of all non-zero binaries Z>0 (SLP ) plus all binaries of associated discount intervals; buckets Fi are obtained by sorting and grouping zero-valued binaries Z=0 (SLP ) according to the reduced costs RC(SLP ) of the associated ordering variables z where n is the number of obtained buckets; obtain (partially relaxed) solution S = MIP (X , Y˜ , K ); obtain generally feasible solution Sbest = MIP (X , Y, K|Z=1 (S ) = 1, S ); i ← 1, lk ← 0,  ← ∞,  ← 0.005, nacc ← 0; while lk ≤ nacc ∨  ≥  do determine total cost of current best solution C old = C (Sbest ); update kernel K = K ∪ Fi by adding bucket i;  −1 update Sbest = MIP (X , Y, K|Z=0 (Sbest ) ∩ ij=1 F j = 0, Sbest ); (S )  = C −C ; C old lif  <  then k ← k + 1; old

11 12 13 14 15

16

best

i ← i + 1; end improve Sbest = MIP (X , Y, Z |Z=0 (Sbest ) = 0 ∧ Z=1 (Sbest ) = 1, Sbest ); return Sbest ;

z-variables appear in the purchasing cost functions, the associated reduced costs RC(SLP ) can be calculated. Our procedure for dynamically constructing the buckets works as follows. First, the vector of reduced costs of all zero-valued z-variables is obtained. As the reduced costs RC (SLP ) are often widely spread, we take their logarithm log(RC (SLP )) to avoid heavy tails in the distribution. To determine suitable buckets, we group variables with similar values of log(RC (SLP )). We first sort the variable vector according to descending values log(RC (SLP )). We then identify groups of variables that are characterized by large first-order differences  log(RC (SLP )) to their neighboring values. For this purpose, we determine a threshold for sufficiently high differences that lead to separate groups of variables. For determining thresholds, we assume that the vector of differences follows some distribution, whereby we are only interested in extreme observations. We choose the Gumbel distribution as the most-prominent extreme-value distribution, which is useful to model the distribution of the maximum of a vector of normally or exponentially distributed observations according to extreme value theory (see e.g. Coles, Bawa, Trenner, and Dorazio, 2001, for an introduction in extreme value theory). More precisely, we define the threshold as the 95%-quantile of Gumbel distribution whereby mean and variance are estimated based on  log(RC (SLP )). Observations exceeding the threshold indicate a new group, i.e., a new bucket. Fig. 2 illustrates the steps of this procedure for instance # 1 which is introduced in the next section. The problem contains 525 binaries variables in total. We obtain a total of 15 buckets (Fig. 2c) that contain between 1 and 196 binary variables. Eventually, we identify the binary variables associated to the z-variables grouped in each such bucket. To avoid overfragmentation, buckets of size less than 5% of the total number of binary variables are merged. Considering the example in Fig. 2c, after merging buckets with less than 0.05 · 525 ≈ 26 variables, we obtain a final set of seven buckets. Table 3 shows the size of the buckets before and after merging.

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

399

Fig. 2. Example for determining the buckets of test instance # 1. Table 3 Size of buckets before and after merging for instance # 1. before merging after merging

A 2

B 15

C 2 26 A

D 7

E 10

F 3

G 1 46 B

H 32

I 28

J 28

K 196

28 C

28 D

196 E

L 14

M 29 43 F

N 28

O 1 29 G

By subdividing the set of binary variables into kernel and buckets, it is to be expected to find a good initial solution quickly. This is based on the assumption that the supplier selection for each material can be determined by analyzing the LP relaxation of the model. Ordering the bucket variables according to their reduced costs is seen as a proxy for their potential to be part of the optimal solution. For the subsequent experimental study, we use the following setting for the parameters of the heuristic: The minimum improvement  is set to 0.5%, the accepted number of non-improving iterations nacc is set to 0, the minimum bucket size is set to 5% of the total number of binary variables. For the MIPs that are solved in steps 5 and 15 of Algorithm 1, different time limits are set and an optimality gap threshold of 1% is accepted. Time limits are set depending on the problem size. As a default, we set time limits of 60 seconds for small instances and 120 seconds for large instances. For the MIP in the kernel search of line 10 of Algorithm 1, time limits are set to 30 and 60 seconds depending on the size of an instance. We also conduct sensitivity analyses for variations of the heuristic’s parameters. Fig. 3. Structure of the case study network.

5. Real world case study 5.1. Case description In this section, we present a case study inspired by a real-world problem from the chemical industry. Fig. 3 shows the structure of the considered system. It consists of a multi-purpose chemical production facility that produces different types of plastic products based on Ethylene and various additives like, for example, plastic granulates. The recipes of the final products vary and determine the demands of Ethylene and the additives. The production plan of the facility is determined for a planning horizon of several weeks, where a period within this horizon corresponds to one day. Accordingly, the demand of each additive is known for each day of the planning horizon, too. While Ethylene is constantly supplied via a pipeline as it is locally produced from Naphtha, the additives are purchased from suppliers and transported to the production facility by truck. In the following, eight different problem instances are considered. The instances vary with respect to their time horizon |T | = {15, 31, 62, 182} and the number of additives considered |M | = { 6 , 9 } . The additives are supplied by three suppliers, each offering one of the discount schemes. Hence, suppliers are labeled SSQD , STQD ,

and SHQD in the following. Each supplier offers a set of additives but no supplier offers all additives. Table 4 shows in columns 2 to price 4 the ranges of purchasing prices csmb/d offered by the suppliers, or ‘—’ if the additive is not offered. The additives are named ‘A1’  to ‘A9’. The further columns show the total demand t∈T amt of each additive m in each instance. The second-last row in this table shows the average daily demand over all additives. The last row indicates the planning horizon length of each instance. Table 5 shows the upper bound of the order quantity for each discount interval offered by each supplier. The lower bound of an interval corresponds to 0 for the first interval and to the upper bound of the preceding interval otherwise. Note that the same interval bounds apply for all additives offered by a supplier. There are four storage facilities I0 to I3 for which Table 6 summarizes the related parameters. The additives used for production in the upcoming period are stored in the pre-production storage facility I0 that consists of four silos with a capacity of 15 tons each. Furthermore, the additives can be stored ahead of production in four nearby silos with a capacity of 100 tons each that are owned by the company (I3 ). The transportation of additives is done by a logistics service provider. This service provider also offers two stor-

400

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406 Table 4 Ranges of purchasing prices and total demand per additive and instance. additive

purchasing prices SSQD

instance #

STQD

SHQD

1

2

3

4

5

6

7

8

A1 60–75 — A2 40–50 39–52 A3 55–63 — A4 76–85 72–90 A5 — — A6 — 27–35 A7 — — A8 86–95 82–99 A9 42–53 — average overall demand per day: number of periods |T |:

62–80 — 62–80 — 38–45 34–40 38–45 — 49–57

78 51 100 100 44 46 — — — 28 15

156 133 186 202 151 102 — — — 30 31

148 107 134 152 111 563 — — — 20 62

195 336 175 384 440 1702 — — — 19 182

72 35 100 94 44 41 — 90 16 34 15

150 115 186 155 121 97 44 128 32 33 31

194 139 176 123 133 110 239 314 108 25 62

340 536 459 332 462 341 557 870 377 24 182

m

Table 5 Discount intervals for the three suppliers considered in the case study. discount

supplier SSQD

supplier STQD

supplier SHQD

interval b/d

q¯ sb

q¯ sd

q¯ sb

q¯ sd

15 40 200 —

60 120 10 0 0 —

1 2 3 4

25 50 75 10 0 0

50 100 150 10 0 0

Table 6 Parameters of all four storage facilities.

the default parameter setting described in Section 4 (indicated by “Heur. (def.)” in subsequent tables) as well as with an extended runtime limit that is twice the limit of the default setting (indicated by “Heur. (ext.)”). All computations are carried out on a single core of a 2.7 gigahertz power linux cluster with 2 gigabyte dedicated RAM. Table 7 shows for each obtained solution, the total costs as well as the split of costs for purchasing additives, transportation, leasing storage, inventory holding, and setting up storage items along with the runtimes and gaps. Gaps are calculated with respect to a reference value Cref as gap =

storage facility i

item type

y¯ i

l¯i

I0 I1 I2 I3

silo silo container silo

4 20 250 4

15 200 25 100

ciinv− f ix

cisetup

0 50 32 0

50 100 10 250

age facilities: A silo farm with 20 silos and 200 tons capacity each (I1 ) as well as a container yard with 250 slots and 25 tons capacity each (I2 ). A daily leasing cost rate of 50 and 32 € is charged per used storage item of storage facilities I1 and I2 . The stock-holding inv−var cost rate cm equals 0.1 € per ton and day for each additive. For storage items in facilities I0 to I3 , setup costs of 50, 100, 10, and 250 € occur when some additive is newly stored in one such item. The additives are shipped via silo truck or container truck with a payload capacity of 25 tons. The interior of the used containers is under pressure and isolated with a special film to avoid contamination of the additives. Additives delivered in containers can be transshipped into silos but it is not possible to refill the containers, which is one reason why transport relations are directed and not complete among all facilities, see Fig. 3. However, the service provider unloads containers only into its own silos (I1 ) or the pre-production silos (I0 ). Fig. 3 shows the resulting transport links where labels attached to arcs indicate the transportation lead times oi,j measured in days. For the sake of brevity further parameters such as cost rates, material demands, and initial stocks are not displayed here but provided in the complete instance files available at https://prodlog.wiwi.uni-halle.de/forschung/research_data/ instanzen_mp_mc_lsp_ss/. 5.2. Computational results 5.2.1. Comparison of solution approaches We first solve all instances by a standard MIP solver (Gurobi 7.5.1) and by the heuristic. The solver is given a runtime limit of 24 hours per instance. The heuristic is coded in R using the gurobi package for interacting with the solver. We test the heuristic with

C−Cre f Cre f

. We use as reference value,

the lower bound obtained from Gurobi when calculating MIP gaps (indicated by “gap (LB)”), the best solution’s objective value delivered by Gurobi (indicated by “gap (grb)”), as well as the heuristic’s final objective value (indicated by “gap (heur.)”). We observe that Gurobi can provide a proven optimal solution for instance # 1 only whereas it consumes the whole allowed runtime for all other instances. In contrast, the heuristic solves all instances within at most a few minutes, delivering very good solutions with a maximum gap of at most 2.1%. For instance # 4, the heuristic even found a better solution within less than seven minutes of runtime than Gurobi within 24 hours. The average gap of the heuristic solutions is merely 0.6%. Some further details of the solution process are presented in Table 8. It shows the gaps between the final solutions of the heuristic and its first feasible solution as well as the kernel search solution (see lines 5 and 10 in Algorithm 1). Additionally, the gap to the solution found by Gurobi at the time when the heuristic terminated is displayed. The first solutions found by the heuristic show costs that are up to 4.8% above the cost of the final solutions. The kernel search reduces this gap to at most 1.2%. Particularly for the larger instances, checking the buckets is a worthwhile feature to improve solution quality. The solutions delivered by Gurobi at the time when the heuristic terminated are similar for instances # 1 and # 5. For instances # 2 and # 6, the heuristic’s solutions are 4.2% and 1.7% better, respectively. For all remaining instances, Gurobi had not found a feasible solution at the time when the heuristic terminated. Hence, solving the monolithic model with Gurobi is reasonable for small instances only whereas the heuristic quickly delivers very good solutions in all cases. Looking at the cost split shown in Table 7, it appears that about 70% of the total costs are for purchasing the additives. The HQD-supplier shows the smallest contribution to total cost with 5% to 12%, whereas the TQD-supplier accumulates most of the purchasing cost in almost all instances with up to 48%. Transportation costs are another important cost category with shares ranging from 19% to 23%. Also inventory cost are substantial, whereby leasing costs (ICy ) clearly dominate stock-holding costs (ICl ). The split of purchasing cost, transportation cost, inventory cost and setup cost

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

401

Table 7 Performance measures of solutions from Gurobi and the heuristic for all case study instances. inst.

total cost

relative shares of total costs per cost category

C

PCSQD

PCTQD

PCHQD

runtime

gap(LB),

1st solution Gurobi

TC

ICy

ICl

ICs

[second]

gap(grb)

26,990 27,004 27,004

31% 31% 31%

30% 30% 30%

time [second]

gap(LB)

7% 7% 7%

19% 19% 19%

9% 10% 10%

1% 1% 1%

4% 4% 4%

7,250 150 300

0.0% 0.1% 0.1%

<1

25.3%

Gurobi Heur. (def.) Heur. (ext.)

65,524 66,143 65,638

27% 27% 27%

32% 32% 32%

10% 9% 10%

21% 21% 21%

7% 8% 7%

1% 1% 1%

2% 2% 2%

86,400 150 300

5.0% 0.9% 0.1%

<1

40.2%

62

Gurobi Heur. (def.) Heur. (ext.)

71,436 71,738 71,589

20% 20% 20%

38% 37% 38%

6% 6% 6%

23% 23% 23%

8% 8% 8%

1% 1% 1%

4% 4% 4%

86,400 361 721

9.0% 0.4% 0.2%

503

39.0%

6

182

Gurobi Heur. (def.) Heur. (ext.)

193,448 192,434 190,695

11% 12% 13%

48% 43% 42%

9% 9% 9%

22% 23% 23%

6% 9% 9%

1% 1% 1%

3% 3% 3%

86,400 366 605

14.0% −0.5% −1.4%

5,279

43.4%

5

9

15

Gurobi Heur. (def.) Heur. (ext.)

29,034 29,153 29,153

22% 22% 22%

32% 32% 32%

5% 5% 5%

22% 22% 22%

11% 12% 12%

1% 1% 1%

7% 7% 7%

86,400 180 360

1.4% 0.4% 0.4%

<1

36.8%

6

9

31

Gurobi Heur. (def.) Heur. (ext.)

66,967 67,328 67,328

23% 23% 23%

31% 31% 31%

8% 8% 8%

21% 21% 21%

11% 12% 12%

1% 1% 1%

4% 4% 4%

86,400 181 300

4.9% 0.5% 0.5%

9

29.4%

7

9

62

Gurobi Heur. (def.) Heur. (ext.)

118,136 120,646 119,919

20% 19% 20%

32% 31% 31%

12% 11% 11%

21% 20% 20%

10% 12% 11%

1% 1% 1%

5% 5% 5%

86,400 364 724

7.9% 2.1% 1.5%

1,836

32.3%

8

9

182

Gurobi Heur. (def.) Heur. (ext.)

365,896 368,285 366,573

19% 18% 19%

36% 34% 35%

11% 10% 11%

21% 19% 20%

8% 13% 10%

1% 1% 1%

4% 5% 4%

86,400 393 743

11.7% 0.7% 0.2%

37,051

32.8%

#

|M|

|T |

1

6

15

Gurobi Heur. (def.) Heur. (ext.)

2

6

31

3

6

4

method

Table 8 Relative gaps of intermediate heuristic solutions and Gurobi’s solution (at the time when the heuristic terminated) to the final heuristic solutions (base setting). inst.

Relative gaps to

#

|M|

|T |

1st sol.

after kernel search

Gurobi (at heur. time)

1 2 3 4 5 6 7 8

6

15 31 62 182 15 31 62 182

1.3% 1.0% 4.8% 2.3% 0.4% 0.9% 1.9% 2.7%

0.0% 0.3% 1.2% 0.9% 0.0% 0.3% 0.6% 1.0%

0.1% 4.2% — — −0.2% 1.7% — —

9

is similar for all solutions. For instances # 5-8, the shares of leasing cost (ICy ) and setup cost (ICs ) increase compared to instances # 1-4 as more raw materials have to be stored. The shares of purchasing cost for the three suppliers vary mostly according to the demand of the supplied raw materials (compare Fig. 6 and Table 4). To study the effect of the bucket handling, we vary the minimum bucket size within [5%, 10%, 20%, 30%] of the total number of binary variables. The accepted number of non-improving iterations nacc is set to [0%, 25%, 50%, 75%, 100%] of the maximum number of iterations. Fig. 4a shows the median and spanwidth of the gaps to the Gurobi solutions over all instances for all tested configurations. We observe here that increasing the number of accepted non-improving iterations tends to lower gaps (better solutions). In contrast, increasing the minimum bucket size increases gaps in tendency. However, the spread over the considered instances is too large to derive reliable statements for both parameters’ effects. Fig. 4b shows the corresponding average runtimes. Here, the effects are as expected: increasing the minimum bucket size decreases average runtimes as less iterations need to be conducted whereas increasing the accepted number of non-improving iterations leads to more iterations, and, thus, larger runtimes.

5.2.2. Structural analysis of solutions This section analyzes the structural properties of the solutions obtained by the heuristic for the instances of the case study. Fig. 5 displays the transport flows for instances # 1 and # 8. We observe that primarily the storage facilities I1 and I2 of the logistics service provider are used in the solutions. Although using the company’s own silos I3 incurs no renting cost, the transport cost rate between I3 and I0 is exceedingly higher than from I1 or I2 to I0 . Also relocating from one storage facility to another is seldom used. Table 9 sheds more light on the storage decisions. It shows the average stock level per stored additive (∅l), the average number of storage items occupied per stored additive (∅y), and the total number of setups over the whole planning horizon ( s). The ∅y columns indicate that if an additive is stored in a storage facility, most often just one storage item is used. Highest average stock levels are observed in I1 with up to 63.7 tons. This effects a fairly small utilization rate given a capacity of 200 tons per silo. On the other side, average stock levels in I3 are even smaller although the capacity of 100 tons per item is also quite large for this facility. In the container facility I2 , capacity utilization is much higher as in the pre-production silo. As the pre-production silos are scarce, more setups take place compared to the other storage facilities. Likewise, the number of setups increases with the number of additives considered (compare  s for instances # 1–4 and instances # 5–8). Furthermore, the supplier selection depends on the demand pattern and discount schemes offered by the suppliers. Fig. 6 shows the shares of total order quantities categorized by instance, supplier, and additive. It appears that for each additive an optimal supplier exists that does not change across the instances. However, the assignments of additives A2 and A8 to STQD indicate that strategic ordering takes place. Although supplier SSQD offers lower prices for small order quantities of additives A2 and A8, both additives are supplied from STQD to achieve higher total quantity discounts. Fig. 7 shows the shares of orders placed in the different discount intervals per supplier and per instance. Here, an order refers

402

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

Fig. 4. Gaps and runtimes of the Kernel search for varying minimum bucket sizes and accepted numbers of non-improving iterations.

Fig. 5. Transport flows of heuristic solutions for instances # 1 and # 8 (width of arcs proportional to logarithmized transport quantities).

to a non-zero quantity of one or more additive(s) supplied in a period. For SHQD each order is associated with selecting two discount intervals (for the total and the lsingle discount) which are both reflected in Fig. 7. As expected, most orders are placed with supplier STQD except in instances # 1 and # 8. Interestingly, the share of orders in high discount intervals becomes smaller with increasing time horizon. This means that saving purchasing cost by achieving higher discounts is exceeded by the required additional leasing costs. Similarly, the share of orders in discount interval 1 is higher for instances with 6 additives compared to their 9-additive counterparts. Again potential savings in purchasing cost discounts do not compensate leasing costs for storage items. To summarize, our case study shows that ordering under discounts and storage selection are highly interdependent decisions. Although an optimal supplier can be identified for each material, the identification of this supplier is not straightforward. Limited storage capacities and heterogeneous leasing costs for stor-

age items also influence supplier selection. In the presented results, savings due to purchasing price discounts are overcompensated by additional storage leasing cost particularly for instances with larger time horizons. Similarly, transportation costs are an important factor particularly for supplier and storage selection. Therefore, as is shown in the case study results, even storage facilities with no leasing costs can be unattractive if the associated transportation cost are too high. Stock-holding costs show only a marginal effect on operational total cost here and, thus, might be neglected too as the stock levels are mainly determined by product demands.

5.3. Semi-artificial instances To study the performance of the heuristic in more detail, we created a set of semi-artificial test instances that are based on the

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

403

Table 9 Usage of storage facilities. inst.

I1

#

∅l

∅y

s

∅l

∅y

s

∅l

∅y

s

∅l

∅y

s

1 2 3 4 5 6 7 8

43.5 63.7 62.5 61.2 30.6 60.7 54.1 38.2

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

3 5 4 22 1 4 12 52

18.5 17.2 14.4 13.0 17.5 18.2 13.5 15.5

1.22 1.29 1.06 1.12 1.15 1.08 1.01 1.04

10 19 22 35 17 17 25 51

— — 32.6 10.1 14.0 14.1 12.9 12.6

— — 1.00 1.00 1.00 1.00 1.00 1.00

— — 3 2 4 4 4 4

5.6 5.9 9.4 7.8 5.2 5.6 6.9 6.6

1.00 1.04 1.21 1.10 1.00 1.00 1.02 1.01

12 19 31 38 15 24 75 223

I2

I3

I0

Fig. 6. Shares of total order quantities per instance, supplier and additive.

Fig. 7. Shares of orders categorized by instance, supplier, and discount interval.

404

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406 Table 10 Average performance measures of the solutions from Gurobi and the heuristic for all artificial instances. Heur. (ext.)

Gurobi (best)

|T |

|S |

# solved

Heur. (def.) time [sec.]

gap (grb)

time [sec.]

gap (grb)

# solved

15

6 12 6 12 6 12 6 12

10 10 10 10 10 10 10 10

175 349 186 365 364 387 620 549

0.2% 0.4% 1.4% 0.6% 2.8% 1.9% 4.1% 3.5%

349 697 360 725 700 723 1038 963

0.2% 0.3% 1.1% 0.4% 2.2% 1.4% 2.0% 1.2%

10 10 10 10 10 10 6 3

31 62 181

Gurobi (first)

Gurobi (at heur.)

time [sec.]

gap (LB)

time [sec.]

gap (LB)

gap (heur)

# solved

gap (heur)

45744 55211 86400 86400 86400 86400 86400 86400

0.5% 0.5% 2.3% 1.2% 4.1% 3.5% 5.6% 3.9%

4 4 87 50 5089 6750 23005 3400

20.7% 12.2% 14.8% 16.2% 22.6% 17.7% 19.6% 35.9%

19.9% 11.2% 10.8% 14.1% 14.4% 11.5% 8.4% 26.7%

10 10 8 10 0 1 0 0

0.1% −0.3% 3.2% 1.0% — 11.9% — —

Table 11 ANOVA results for the effects of minimum bucket size and allowed number of non-improving iterations. Intercept

Estimate 1.144 Std. Error 0.032 ∗∗∗ Sig. level Model stats: R2 = 0.328; Note: ∗ p < 0.1;

∗∗

p < 0.05;

∗∗∗

nacc

min. bucket size

25%

50%

75%

100%

10%

20%

30%

−0.308 0.025

−0.482 0.025

−0.532 0.025

−0.568 0.025

0.147 0.022

0.293 0.022

0.266 0.022

∗∗∗

∗∗∗

∗∗∗

∗∗∗

∗∗∗

∗∗∗

∗∗∗

adj. R2 = 0.327;

RSE = 1.343;

F-stat = 1020∗∗∗

p < 0.01; setting effects estimated but omitted in table.

case study instances #5–#8. For each of the instances #5–#8, the set of suppliers and their corresponding parameters are now generated randomly. We consider an enlarged set of 6 and 12 suppliers for each instance with 2 and 4 suppliers per discount scheme, respectively. For each of the resulting eight settings, 10 instances are randomly generated, which leads to a total of 80 test instances. All instance files can be found on https://prodlog.wiwi.uni-halle. de/forschung/research_data/instanzen_mp_mc_lsp_ss/. In each instance, each supplier offers between 1 to 9 randomly drawn materials. lTransport relations between suppliers and storage facilities are randomly sampled, too, while the transport relations between the storage facilities are the same as in the case study (see Fig. 3). For each supplier, a random transport time is sampled from the interval [1, 6] along with the supplied storage facilities. The sampled transport time of a supplier applies to all of its transport relations. The transport cost rates per relation are normally distributed with μ = 50 and σ = 15. The number of offered discount intervals is sampled from the interval [2, 6]. The overall average unit prices R for all materials are drawn from a normal distribution with μ = 60 and σ = 20 and serve as a reference for all corresponding suppliers. Then, each supplier offers an individual lbase line unit price that is calculated by multiplying the overall average unit price with a normally distributed multiplier (B ∼ N(1, 0.1)). Hence, each supplier offers a base line price that varies by about 10% from the overall average unit price and it is equally likely that the base line price is below or above the overall average unit price. The discount levels are determined in a similar way by drawing multipliers from N ∼ N(0, 0.1) for SQD and TQD suppliers or N ∼ N(0, 0.05) for HQD suppliers. The price interval offered by a supplier is then calculated as [R · B, R · B · (1 + |N| )] if B < 1 and [R · B · (1 − |N| ), R · B] if B > 1. I.e., if the base line price R · B is below the overall average unit price R, the base line price serves as the minimal price offered while it is the maximum price otherwise. This way the price intervals offered by the suppliers are overlapping to a large extend. Finally, the discount intervals are constructed for each material-supplier combination by sampling a minimal ordering quantity from the interval [0,5] and an interval width from [5,20]. I.e., for each material-supplier combination, all

discount intervals have the same width (except for the last interval whose upper bound is at least 50). Table 10 shows the average performance measures over the 10 instances per setting. It reports results for the heuristic in the default setting and with the extended runtime limit as well as for Gurobi. The results confirm the preliminary conclusions drawn from the case study instances. Overall, the heuristic performs well in comparison to Gurobi with an average gap to the best Gurobi solutions of at most l4.1% (base setting). When the runtime limit is doubled in the extended setting, results for small instances with |T | ∈ {15, 31} improve only slightly, while for the remaining instances noteworthy improvements can be realized (compare columns “gap (grb)”). Furthermore, Gurobi delivers feasible solutions only for l29 out of the 40 largest instances (|T | = {62, 182}) whereas the heuristic can solve all instances within the given time limits (compare columns “# solved”). Considering Gurobi’s performance at the time when the heuristic finishes in the default setting, the performance is almost identical for |T | = 15 or even slightly worse for |T | = 31. For the larger instances, Gurobi solves only 1 out of 40 instances within this time span. lIt is noteworthy that for the instances with 12 suppliers average gaps to the lower bound are consistently smaller than for the instances with 6 suppliers. This indicates that larger numbers of suppliers allow the solver to calculate tighter lower bounds. On the other side, instances with more suppliers seem to be more complex to solve to optimality. This is indicated by smaller average gaps of the heuristic’s solutions for instances with 12 suppliers (compare columns “gap (grb)”). For instances with |T | = 15, moreover, average solution time is larger for the instances with 12 suppliers than for instances with 6 suppliers (compare column “time” in first two rows). The effects of bucket handling are also examined using the semi-artificial instances. For this purpose, all 80 instances are solved for all of the 20 configurations regarding minimum bucket size and allowed number of non-improving iterations as outlined above. To analyze the results across the tested configurations and instances, gaps to the best found solution are calculated. The best solutions are delivered either by Gurobi or the tested standard

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

405

Table 12 Average number of used suppliers and frequency of suppliers per material.

|T |

|S |

15

avg. number of

Frequency of suppliers per material

used suppliers

1

2

3

4

85.0% 71.7% 70.0% 65.4% 72.2% 48.9% 67.9% 33.3

15.0% 28.3% 27.5% 32.1% 25.6% 41.1% 27.2% 51.9%

0.0% 0.0% 2.5% 2.5% 2.2% 8.9% 4.9% 13.6%

0.0% 0.0% 0.0% 0.0% 0.0% 1.1% 0.0% 1.2%

6 12 6 12 6 12 6 12

31 62 181

2.7 3.2 3.6 4.1 3.6 5.2 2.7 4.1

configurations of the heuristic (as summarized in Table 10). Table 11 displays the results of an Analysis of Variance (ANOVA) of the gaps to the best solutions. Detailed results in terms of the medians over all instances in a certain setting and for each tested configuration can be found in the appendix in Table A.1. The reference setting (coded in the “Intercept” parameter) refers to the default setting with 5% minimum bucket size and only improving iterations allowed. Table 11 confirms the observations of the case study experiments. Increasing the allowed number of non-improving iterations leads to significantly reduced gaps. Thereby, the marginal effects are largest if this parameter is set to lthe maximum number of iterations. Similarly, minimum bucket sizes larger than lthe reference value of 5% lead to significantly higher gaps. Hence, we conclude that testing up to 50% of the available buckets contributes to an improvement of solution quality. Additionally, the minimum bucket size lshould be set to 5% of the total number of binary variables. To analyse the role of suppliers, Table 12 displays the average number of suppliers engaged in the heuristic’s solutions as well as the frequency of using 1, 2, 3, or l4 suppliers for sourcing a material over all instances and materials. The results indicate that a single supplier is used for most of the materials. Still, the frequency of multiple sourcing (i.e., at least two suppliers are engaged for a material) as well as the total number of engaged suppliers increases with an increasing number of available suppliers. Both effects intensify with increasing time horizons. To summarize, the proposed heuristic delivers consistently good results with very low optimality gaps. The method even

outperforms Gurobi in terms of the number of solved instances and solution time for larger problems. Furthermore, for large instances, extending the runtime of the heuristic can considerably improve solution quality. 6. Conclusions In this manuscript a multi-period lot-sizing problem with supplier selection, discounts and storage selection is described, which appears for example in the chemical industry. We proposed an optimization model, a kernel search based heuristic, and a comprehensive case study for this problem. It turns out that a standard solver has difficulty with solving the problem to optimality, whereas our heuristic can solve all instances in at most seven minutes with very low gaps. The case study revealed a strong interdependence between supplier selection and storage selection as logistics cost contribute noticeably to total cost. It also showed that the offered discount schemes have an impact on the supplier selection and ordering decisions. Therefore, integrating all these fields of decision making offers managers additional saving potentials in operational procurement management. For future research, incorporating stochasticity in demand and supply processes is an extension, which deserves attention as it is of high relevance in practice also in the process industries. Likewise, studying the dynamic behavior of supply systems when data change in time is a worthwhile field of future research. Appendix A. Detailed results of the semi-artificial instances

Table A1 Medians of the gaps to the best solution for all tested configurations of the heuristic.

|B |

5%

10%

20%

30%

nacc

0%

25%

50%

75%

100%

0%

25%

50%

75%

100%

0%

25%

50%

75%

100%

0%

25%

50%

75%

100%

15 6 15 12 31 6 31 12 62 6 62 12 181 6 181 12 aver.

0.4% 1.0% 1.6% 0.9% 1.9% 2.0% 2.7% 1.5% 1.5%

0.2% 0.1% 0.9% 0.4% 1.4% 1.1% 1.9% 1.5% 0.9%

0.0% 0.0% 0.7% 0.3% 1.1% 0.8% 1.7% 0.8% 0.7%

0.0% 0.0% 0.9% 0.2% 1.5% 0.8% 1.4% 0.9% 0.7%

0.0% 0.0% 0.7% 0.3% 1.7% 0.8% 1.3% 0.9% 0.7%

0.4% 0.5% 1.5% 0.8% 2.1% 1.8% 2.8% 1.5% 1.4%

0.2% 0.0% 1.3% 0.4% 1.8% 1.0% 2.7% 1.5% 1.1%

0.1% 0.0% 1.0% 0.3% 1.8% 0.9% 1.8% 1.5% 0.9%

0.1% 0.0% 0.9% 0.2% 1.6% 0.9% 1.8% 1.5% 0.9%

0.1% 0.0% 1.0% 0.2% 1.6% 0.9% 1.7% 1.2% 0.8%

0.2% 0.3% 1.2% 0.6% 2.2% 1.5% 2.1% 1.5% 1.2%

0.2% 0.0% 1.3% 0.5% 1.6% 1.3% 2.8% 1.5% 1.2%

0.1% 0.0% 1.2% 0.4% 1.8% 1.2% 2.5% 1.5% 1.1%

0.2% 0.0% 1.1% 0.4% 1.6% 1.1% 1.9% 1.5% 1.0%

0.2% 0.0% 1.2% 0.4% 1.7% 0.9% 2.0% 0.7% 0.9%

0.2% 0.3% 1.1% 0.6% 1.9% 1.3% 2.6% 2.2% 1.3%

0.1% 0.1% 1.1% 0.5% 1.8% 1.2% 2.4% 2.2% 1.2%

0.0% 0.0% 1.1% 0.4% 1.7% 1.1% 2.3% 2.1% 1.1%

0.0% 0.0% 1.3% 0.3% 1.8% 1.2% 2.1% 2.0% 1.1%

0.0% 0.0% 1.1% 0.4% 1.4% 1.1% 2.2% 1.9% 1.0%

406

T. Kirschstein and F. Meisel / European Journal of Operational Research 279 (2019) 393–406

References Aissaoui, N., Haouari, M., & Hassini, E. (2007). Supplier selection and order lot sizing modeling: A review. Computers & Operations Research, 34(12), 3516–3540. Angelelli, E., Mansini, R., & Speranza, M. (2007). Kernel search: A heuristic framework for MILP problems with binary variables. Technical report, Department of Electronics for Automation, University of Brescia. Angelelli, E., Mansini, R., & Speranza, M. G. (2010). Kernel search: A general heuristic for the multi-dimensional knapsack problem. Computers & Operations Research, 37(11), 2017–2026. Angelelli, E., Mansini, R., & Speranza, M. G. (2012). Kernel search: A new heuristic framework for portfolio selection. Computational Optimization and Applications, 51(1), 345–361. Balakrishnan, A., & Natarajan, H. P. (2014). Integrated procurement planning in multi-division firms. Production and Operations Management, 23(10), 1795–1810. Benton, W., & Park, S. (1996). A classification of literature on determining the lot size under quantity discounts. European Journal of Operational Research, 92(2), 219–238. Bohner, C., & Minner, S. (2017). Supplier selection under failure risk, quantity and business volume discounts. Computers & Industrial Engineering, 104, 145–155. Chaudhry, S. S., Forst, F. G., & Zydiak, J. L. (1993). Vendor selection with price breaks. European Journal of Operational Research, 70(1), 52–66. Coles, S., Bawa, J., Trenner, L., & Dorazio, P. (2001). An introduction to statistical modeling of extreme values, Springer, 208. Crama, Y., & Torres, A. (2004). Optimal procurement decisions in the presence of total quantity discounts and alternative product recipes. European Journal of Operational Research, 159(2), 364–378. Federgruen, A., & Lee, C. Y. (1990). The dynamic lot size model with quantity discount. Naval Research Logistics, 37(5), 707–713. Ghaniabadi, M., & Mazinani, A. (2017). Dynamic lot sizing with multiple suppliers, backlogging and quantity discounts. Computers & Industrial Engineering, 110, 67–74. Goossens, D. R., Maas, A., Spieksma, F. C., & Van de Klundert, J. (2007). Exact algorithms for procurement problems under a total quantity discount structure. European Journal of Operational Research, 178(2), 603–626. Guastaroba, G., Savelsbergh, M., & Speranza, M. (2017). Adaptive kernel search: A heuristic for solving mixed integer linear programs. European Journal of Operational Research, 263(3), 789–804. Hadley, H., & Whitin, T. M. (1963). Analysis of inventory systems. Prentice Hall. Jackson, J. E., & Munson, C. L. (2016). Shared resource capacity expansion decisions for multiple products with quantity discounts. European Journal of Operational Research, 253(3), 602–613.

Lee, A. H., Kang, H.-Y., Lai, C.-M., & Hong, W.-Y. (2013). An integrated model for lot sizing with supplier selection and quantity discounts. Applied Mathematical Modelling, 37(7), 4733–4746. Low, R., & Waddington, J. (1967). The determination of the optimum joint replenishment policy for a group of discount-connected stock lines. Journal of the Operational Research Society, 18(4), 443–462. Manerba, D., & Mansini, R. (2012). An exact algorithm for the capacitated total quantity discount problem. European Journal of Operational Research, 222(2), 287–300. Manerba, D., & Mansini, R. (2014). An effective matheuristic for the capacitated total quantity discount problem. Computers & Operations Research, 41, 1–11. Manerba, D., Mansini, R., & Perboli, G. (2018). The capacitated supplier selection problem with total quantity discount policy and activation costs under uncertainty. International Journal of Production Economics, 198, 119–132. Manerba, D., & Perboli, G. (2019). New solution approaches for the capacitated supplier selection problem with total quantity discount and activation costs under demand uncertainty. Computers & Operations Research, 101, 29–42. Mansini, R., Savelsbergh, M. W., & Tocchella, B. (2012). The supplier selection problem with quantity discounts and truckload shipping. Omega, 40(4), 445–455. Mazdeh, M. M., Emadikhiav, M., & Parsa, I. (2015). A heuristic to solve the dynamic lot sizing problem with supplier selection and quantity discounts. Computers & Industrial Engineering, 85, 33–43. Munson, C. L., & Jackson, J. (2015). Quantity discounts: An overview and practical guide for buyers and sellers. foundations and trends in technology. Information and Operations Management, 8(1–2), 1–130. Munson, C. L., & Rosenblatt, M. J. (1998). Theories and realities of quantity discounts: An exploratory study. Production and Operations Management, 7(4), 352–369. Parsa, I., Khiav, M., Mazdeh, M., & Mehrani, S. (2013). A multi supplier lot sizing strategy using dynamic programming. International Journal of Industrial Engineering Computations, 4(1), 61–70. Shahsavar, A., Zoraghi, N., & Abbasi, B. (2018). Integration of resource investment problem with quantity discount problem in material ordering for minimizing resource costs of projects. Operational Research, 18(2), 315–342. Stadtler, H. (2007). A general quantity discount and supplier selection mixed integer programming model. OR Spectrum, 29(4), 723–744. Tempelmeier, H. (2002). A simple heuristic for dynamic order sizing and supplier selection with time-varying data. Production and Operations Management, 11(4), 499–515. Xu, J., & Lu, L. L. (1998). The dynamic lot size model with quantity discount: Counterexamples and correction. Naval Research Logistics, 45(4), 419–422.