1. Introduction

Many companies strive to synchronize their production, transportation and replenishment planning by adopting supply chain management practices such as vendor-managed inventory, efficient consumer response, and collaborative planning, forecasting and replenishment. The main objective of these collaborative methods is to reduce inefficiencies and to eliminate redundancies between the different partners. Moreover, collaborative methods allow suppliers to have better visibility of their clients forecast and planned demand and thus be able to plan more effectively ahead of time. However, to take full advantage of these methodologies, companies are in need of decision-making tools to help them turn the large amount of shared data to smart actions.

In this paper, we tackle the problem of integrating and synchronizing production, distribution and inventory planning in a single manufacturing plant and multiple delivery points network with dynamic demands. In practice, a manufacturing plant usually supplies multiple locations. These locations can be distribution centres and/or final demand points. The nodes of the production–distribution network can belong to a single company or to several companies. Transportation costs between the plant and its client-locations often offer economies of scale that can be represented by piecewise linear functions. These transportation economies of scale may have a major impact on inventory planning and replenishment strategies at both the plant and its clients. Major cost savings can be achieved by integrating production, inventory and transportation planning.

Integrating lot-sizing decisions in a single-origin multi-destination distribution network has attracted the attention of researchers. In these problems, the origin is responsible for replenishing multiple locations (clients) with finished goods to satisfy planned and forecasted orders. Anily (1994), Gallego and Simchi-Levi (1990), Anily and Federgruen (1990, 1993) and Herer and Roundy (1997) tackled the problem in the case of a single product with deterministic static demand. In their work, the origin is assumed to be a distribution centre and the destinations represent multiple retailers. Retailer orders are assumed to be small and hence multiple orders should be consolidated to take advantage of transportation economies of scale. Transportation costs are made up of a cost per mile plus a fixed charge for using a truck. The objective is to determine replenishment policies that specify the delivery quantities and the vehicle routes so as to minimize long-run average inventory and transportation costs. Viswanathan and Mathur (1997) generalized Anily and Federgruen (1990) to the multi-item version of the problem. Diaby and Martel (1993) and Chan et al (2002) assumed deterministic dynamic demand and general piecewise linear transportation costs for the single product case.

Most of the above contributions tackled the flow integration problem for a single item, in a single-origin multi-destination network. Also, they approach the problem from the retailer's perspective since the origin is a distribution centre. However, often in practice the origin is a plant where multiple products share the same production and transportation resources.

In this paper, we consider the integration of production, distribution and inventory planning of multiple products in a single-plant multi-destination network of the type illustrated in Figure 1. This production–distribution context is found in the pulp and paper industry and in several other process industries. Production at the manufacturing location is assumed to involve multiple stages but with one of them creating a bottleneck. In the bottleneck stage, a number of intermediate products (IPs) are manufactured by parallel machines, each machine having a predetermined fixed production sequence. In the succeeding stages, the IPs are transformed into a large number of finished products (FP). No inventory of IPs is kept but the FPs can be stocked at the plant before they are shipped. An inventory of FPs is also kept in the distribution centres.

Figure 1.

Single-plant multi-destination network of a paper manufacturer.

Full figure and legend (74K)

Multiple transportation modes can be used to replenish the different destinations with finished goods. These transportation modes may have different transportation lead times from the plant to each destination and their cost structure can be represented by a general piecewise linear function to reflect economies of scale. Transit inventory costs are explicitly taken into consideration and are embedded in the cost structure of each transportation mode as shown in Figure 2. It is important to note that when transportation modes have the same lead time to a given destination, their cost structures can be represented by a single piecewise linear function. However, if the transportation modes have different transportation lead times then their cost structures must be represented by different piecewise linear functions as shown in Figure 2.

Figure 2.

An example of two different transportation modes cost functions.

Full figure and legend (14K)

The rest of the paper is structured as follows. In Section 2, a mixed integer programming formulation of the production–distribution planning problem considered is proposed. In Section 3, valid inequalities to strengthen the proposed formulation are derived and the strategy of adding extra 0-1 variables to strengthen the branching process is discussed. A sequential solution approach based on the independent, but synchronized, solution of the production and distribution problem is also introduced. Finally, in Section 4, using data from the production–distribution network of a large international paper company, the computational performances of the different solution strategies proposed are tested and the advantage of integrated planning over independent but synchronized planning is assessed.

2. Problem formulation

In this section, the problem studied is defined more precisely and a mathematical programming formulation of the production–distribution process is proposed.

2.1. Production problem

In the manufacturing plants of the type illustrated in Figure 1, there is typically a significant IP changeover effort (cost and time) required and the intermediate production context is very similar to the well-known discrete lot-sizing problem (Salmon, 1991) where short time periods are considered and at most one production setup/changeover is allowed per machine per period of time. On the other hand, the production of a larger number of FPs does not involve significant setup/changeover costs and their production resources do not represent bottlenecks.

Let g_ii' be the number of units of IP product i required to produce one unit of FP product i', taking into account any waste incurred in the transformation process. Since each FP is made from a single IP, the set of FP products can be partitioned according to the IP product it is made of. In addition, it is assumed that a standard production sequence of IP products must be maintained for each machine m=1,...,M, and that at most one product type can be produced in a given time period. Let e_m denote the index of the IP product in the eth position in machine m production sequence, e_m=1_m,...,f_m, where f_m represents the product in the final position of machine m production sequence. Thus, when l<f product (e+1)_m can be produced on machine m, only after product e_m has finished its production batch. The production resource consumption for IPs is assumed to be concave, that is, a fixed resource capacity consumption is incurred whenever production switches from one IP product to another (changeover time), and linear resource consumption is incurred during the production of a batch of IP products (see Figure 3). Inventory holding costs are assumed to be linear. The notation required to formulate the production planning model is summarized in Tables 1 and 2.

Figure 3.

An example of a production sequence at machine m.

Full figure and legend (52K)

Table 1 - Production indexes, parameters and sets.

Full table

Table 2 - Production decision variables.

Full table

Using this notation, the production planning problem of the manufacturing location can be represented by the following optimization model:

In model (P), (P.1) and (P.2) are the flow conservation constraints of IP and FP products at the manufacturing location. Constraints (P.3) ensure that production capacity is respected. Constraints (P.4) and (P.5) make sure that the production sequence is respected for each machine. For a given machine m, when e<f, constraint (P.4) enforces the number of product (e + 1)_m changeovers to be less than or equal to the number of product e_m changeovers for any given period of time. Hence, it forces product (e + 1)_m production to start only after the production batch of product e_m is completed. Constraints (P.5) do the same job for product f_m which has the particularity of being the last of the production sequence on machine m. Thus, after its production batch, machine m has to switch production to product 1_m and start another sequence. Constraints (P.6) makes sure that at most one product is manufactured per period of time for each machine. Finally, constraints (P.7) restrict the changeovers on a machine to the periods in which there is some production.

2.2. Transportation costs

For most practical situations, joint shipping cost structures can be represented by a general (ie, not necessarily concave, nor continuous) piecewise linear function z(S) of the load S (measured in cwt, tons, pallets, etc) to be shipped. These general piecewise linear functions can be represented as a series of linear functions, as shown in Figure 4. Let S_j, j=0,...,,(S₀=0) denote the break points of the piecewise linear function and let b_j=S_j-S_j-1,j=1,..., denote the length of the jth interval on the S-axis defined by the break points (S₀,...,S). Finally, for interval j, let v_j be the slope of its straight line (variable cost), A_j be the discontinuity gap at the beginning of the interval and E_j be the value of the function at the end of the interval, that is, E_j=z(S_j). Then, it is seen that for S_j-1<SS_j, we have z(S)=(E_j-1+A_j)+v_js_j,s_j=(S-S_j-1).

Figure 4.

General piecewise linear function.

Full figure and legend (48K)

In Rizk et al (2006b), three different formulations to model this type of piecewise linear functions were compared. Based on extended experimental tests, it was shown that one of these formulations performs more efficiently than the others. In what follows, this formulation is used to model our piecewise linear transportation costs. For an amount S to be shipped to a given destination, using a given transportation mode, in a given period of time, let j be the interval for which S_j-1<SS_j,j1,S₀=0. The amount shipped can be expressed as S=_jS_j, where _j=S/S_j for j1. Based on the above, S can be written in general as S=_{j=0 j}S_j, where (S_j-1/S_j)<_j1 if S_j-1<SS_j and _j=0 otherwise, for j=1,...,. The last two conditions can be represented by a binary variable _j, where , j=1,...,, and . Using the above observation, S can be expressed in an LP model by the following set of constraints:

From the definition of E_j and the above set of constraints, it is seen that z(S) can be expressed as a linear function of variables _j and _j,j=1,...,:

The additional notation required to formulate the transportation and replenishment planning model is summarized in Tables 3 and 4.

Table 3 - Transportation and replenishment indexes, parameters and sets.

Full table

Table 4 - Transportation and replenishment decision variables.

Full table

2.3. Production–distribution model

For a given destination wW, let ^w=Min_uU^w(_u^w)^w is the shortest transportation lead time to destination w. Let us assume that in period 1, a quantity of product iFP is manufactured at the plant and that at the end of period 1, we decided to ship an amount of product i to destination w. Because of the production and transportation lead times, the quantity of product i shipped cannot get to destination w earlier than time period +^w+1. Thus, destination w replenishment planning should at the earliest start at period +^w+1. Figure 5 illustrates an example of different transportation modes shipments (R_uit^w) to satisfy product i demand at destination w. In this example, the planning horizon has five periods (T=5). There are three available transportation modes to ship FPs from the plant to destination w(U^w={1,2,3}). The transportation modes lead times from the plant to destination w are ₁^w=1,₂^w=2 and ₃^w=3. Production planning in the plant starts at period 1 and ends at period T=5. On the other hand, because of production and transportation lead times, as stated above, replenishment planning in destination w starts at period +^w+1 and ends at period T++^w. Note that for a given transportation mode uU^w, only shipments that are made before time period T++^w-_u^w can get to destination w within its replenishment planning horizon ([+^w+1,T++^w]).

Figure 5.

Example of multiple transportation modes shipments.

Full figure and legend (64K)

For a given destination wW and transportation mode uU^w, let _u^w=_u^w-_u^w. Using the notation in Tables 1, 2, 3 and 4, the synchronized production–distribution planning problem in a single-manufacturer multi-destination network with multiple transportation modes can be formulated as follows (F): (Figure 5):

subject to (P.1), (P.3), (P.4), (P.5), (P.6), (P.7), (P.8) and to constraints:

The objective function in model (F) is to reduce operations costs across the whole network. Operations costs include changeover costs at the plant, inventory costs at both the plant and the destinations and transportation costs. The plant flow balance constraint (F.1) is similar to (P.2) but it includes the different destinations shipment quantities. (F.2) is the flow balance constraint at each destination. Following relations (1), (2) and (3), constraints (F.3), (F.4) and (F.5) define the transportation capacity and cost structure between the plant and the different destinations wW, for each transportation mode uU^w.

3. Model solution

In this section, valid inequalities to strengthen the formulation (F) are proposed. The use of cuts to improve the solution of lot-sizing problems was first introduced by Barany et al (1984). The cuts proposed here, however, are derived from the specific properties of our problem. We also show how a reformulation technique proposed by Orgyczak (1996) for multiple choice requirement constraints can be applied to our problem. Orgyczak's reformulation technique is known in the literature as special ordered inequalities (SOI). The strategy of adding extra 0-1 variables to improve the branching process is also discussed. Lastly, we propose a heuristic sequential solution approach to solve large versions of program (F).

3.1. Valid cuts

Rizk et al (2006b) presented five valid cuts for a synchronized production–distribution planning problem with a single destination and where all transportation modes are assumed to have the same lead time. Three of these valid cuts are 'transportation' cuts involving the binary variables associated with a piecewise linear transportation cost function. In this section, we generalize two of these valid cuts to the case of multiple destinations and multiple transportation lead times. These generalizations are presented in proposition 1 and proposition 2 below. Two 'production' cuts, involving only binary variables associated to the production set-up, are also presented in Rizk et al (2006b). These cuts are also valid for (F) and they are presented below as proposition 3 and proposition 4.

Proposition 1

are valid cuts for (F)

Proof

Let _u(t-u^w)0^w=1 for a given location wW, a transportation mode uU^w and a period of time +_u^w+1tT++_u^w. In this case, _u(t-u^w)j^w=0, 1j_u(t-u^w)^w (constraints (F.5)) and hence R_ui(t-u^w)^w=0 for iFP given (F.3) and (F.4). Thus, based on (F.2) and (F.6), .

(Cut.1) is based on the observation that, if no shipments using transportation mode u are made to destination w in period (t-_u^w), there must be enough quantity of each product shipped via other transportation modes and/or there is enough inventory in the destination at the end of period (t-1) to satisfy the demand of period t.

The next valid cut (Cut.2) forces an explicit lower bound on the amount shipped during a time interval and hence tightens the constraints on the interval selection binary variables.

Proposition 2

The following inequalities are valid cuts for (F)

where wW, t=+^w+1,...,T++^w and tLT++^w.

Proof

For a location wW and a finished product iFP, we know that I_i(t-1)^w+_t'=t^L[_uU^wR_ui(t'-u^w)^w]_t'=t^Ld_it'^w (since this is an aggregate sum of constraints (F.2) and since I_iL^w0) and hence r_iI_i(t-1)^w+_t'=t^L[_uU^wr_iR_ui(t'-u^w)^w]_t'=t^L r_id_it'^w. If we sum the previous inequalities across all finished products iFP, we obtain _iFP[r_iI_i(t-1)^w+_t'=t^L[_uU^wr_iR_ui(t'-u^w)^w]]_iFP[_t'=t^Lr_id_it'^w]. Furthermore, we can observe from constraints (F.3) and (F.4) in (F) that _iFPr_iR_ui(t'-u^w)^w _j=1^w_u(t'-uw)S_u(t'-u^w)j^w_u(t'-u^w)j^w. Thus: (Cut.2) is valid.

An important special instance of the above cut is the case where for a given location w, there is only one mode of transportation u in U^w, and for all time period t the interval lengths b_utj^w are the same, that is, b_utj^w=b_u^w, S_utj^w=jb_u^w. In this case, since I_i(+^w)^w=0, (Cut.2) for destination w, for time period t=+^w+1, and for tLT++^w can be reduced to:

and since the left side is integer, the above inequalities can be replaced by:

Proposition 3

For an intermediate product iIP and a period 1tT, the following inequalities:

where CAP_it=Max_(mMi),(1t't) (C_it'^m/a_it'^m), are valid cuts for (F). See proof in Rizk et al (2006b).

Based on the total demand of a given product during a time interval and its available production capacity, (Cut.3) forces a tight lower bound on the number of product setups that need to occur during the same time interval:

Proposition 4

For each machine mM, each product iIP, and each period 1tT, the following inequalities are valid cuts for (F):

See proof in Rizk et al (2006b).

(Cut.4) simply states that an intermediate product i changeover must occur (_it^m=1) at the beginning of period t if a new production batch for this product is initiated (_i(t-1)=0 and _it=1) at the beginning of the same period.

3.2. SOI

In formulation (F), for a given destination wW, a transportation mode uU^w and a time period t, variables _utj^w, j=0,...,_ut^w form a Special Order Set (SOS) of type 1 as defined by Beale and Tomlin (1970). To model SOS1 more efficiently in a MIP program, Orgyczak (1996) introduced a set of new variables _utj^w where _ut0^w=_ut0^w, _utj^w=_ut(j-1)^w+_utj^w for j=1,...,_ut^w, and he proposed the following variable upper bounding (VUB) constraints to replace (F.5):

This formulation of multiple choice requirement constraints is known as SOI. Orgyczak (1996) reported, based on experimental tests, that SOI reduces the CPU time of simple MIP solvers. Branching on the new variables introduced by the SOI is expected to result in a smaller and more balanced branching tree and thus accelerate the convergence of the branch-and-bound procedure.

Note that SOI can be used in formulation (F) simply by replacing _utj^w by (_utj^w-_ut(j-1)^w and constraints (F.5) by (F.9). Let us denote this new formulation as (F-SOI). The impact of using the SOI will be tested in our experimental test section.

3.3. Special 0-1 variables

The number of 0-1 variables is usually regarded as an indicator of the computational difficulty of models of the type proposed in this paper. This is because the potential size of the branch-and-bound tree increases exponentially with the number of 0-1 variables. However, the strongest virtue of the branch-and-bound method is its ability to eliminate large sections of the branch-and-bound tree as was shown in practice. There is even sometimes benefit in adding extra 0-1 variables that represent 'dichotomies' in the system being modelled to facilitate the branching process. Based on this observation, for a destination wW and a transportation mode uU^w, we introduce new binary variables

These new variables can be used to represent the dichotomy 'either the quantity shipped is higher than S_utj^w or not'. For formulation (F), X_utj^w can be expressed as:

For formulation (F-SOI), X_utj^w can be expressed as:

Branching on these new 0-1 variables has the effect of a cutting plane. In our test section, in addition to testing the performance of the valid cuts, the impact of adding the above 0-1 variables on the computation time is also tested.

3.4. Sequential solution approach

The valid cuts proposed above can help enormously in solving (F) with commercial solvers such as Cplex as shown in the experimental section at the end of the paper. However, for some practical cases, solving (F) with available commercial solvers can still be difficult. Hence, it is worth developing a heuristic approach to solve larger versions of (F). The approach proposed is based on a sequential solution strategy often used in practice (Chandra and Fisher, 1994). A production plan is first prepared. The objective of this plan is to minimize the changeover and plant inventory holding costs, subject to meeting the network total effective demand per period. Replenishment plans for all destinations are then elaborated, one by one, so as to minimize each destinations total inventory holding and transportation costs, subject to the inventory availability specified in the production plan and to the time-varying demands of each destination. The following sub-sections describe further the production planning and the destinations replenishment sub-models.

3.4.1. Plant production plan

Let D_it=d_it+_w=1^W d_i(t+^w)^w, t=+1,...,T+ be the network total effective demand for product iFP in period t. The plant production plan can then be determined using the following model:

Model (P') is similar to model (P) described in Section 2, except that in constraint (P'.2), the plant demand for each product i, d_it, is replaced by D_it to take the effective demand of the whole network into consideration. (P') is then solved using Cplex 9.0 libraries, in combination with the valid cuts (Cut.3) and (Cut.4) proposed in Section 3.1.

3.4.2. Destinations replenishment plan

Based on the plant production plan, replenishment to the different destinations wW must be planned in a way that minimizes their inventory holding and transportation costs subject to meeting the products demand in each time period. The replenishment plan of destination w specifies the shipping quantities R_uit^w. In our sequential solution approach, the different destination replenishment problems are solved one by one as they are listed: the problem for destination 1 is solved first, then the problem for destination 2, etc. While solving a replenishment problem for a given destination w, one must take into account the replenishments already planned for all destinations w'<w, as well as make sure that enough inventory will still be available at the plant to meet demand for all destinations w'>w. This way, the different destination replenishment plans are synchronized and the inventory levels planned at the plant are respected.

Let P_it^w represent the quantity of product i demanded, in period t, by sites other than destination w. P_it^w can be expressed as P_it^w=d_it+_y<wuU^yR_uit^y+_y>wd_i(t+^y). For destination w, the replenishment plan can be determined using the following model:

Model (D-w) objective is to reduce inventory holding and transportation costs of destination w. Constraints (D-w.1) and (D-w.3) to (D-w.8) are similar to constraints (F.2) to (F.8). Constraints (D-w.2) make sure that inventory availability as implied by the production plan of model (P') is respected.

To strengthen the mathematical formulation of the different destination replenishment models, the valid cuts (Cut.1) and (Cut.2) proposed in Section 3.1 can be used. Note that when multiple transportation modes are considered, model (D-w) can be very difficult to solve even if it applies to a single destination. In that case, the valid inequalities (Cut.2) are not as efficient as for the case of a single mode. However, one can observe that in that case, only the fastest transportation mode can be used to deliver shipments and meet demand of destination wW for periods t=+^w+1,...,+₂^w. Let, without loss of generality, mode 1 and mode 2 designate the fastest (₁^w=^w) and the second fastest (₂^w) transportation modes to destination wW. Then, if for all time periods t=+^w+1,...,+₂^w, the interval lengths b_1tj^w are the same, the strengthened cuts (Cut.2'), can be used instead of (Cut.2) for periods t=+^w+1,...,+₂^w and for tL+₂^w. Models (D-w), wW, are then solved using Cplex 9.0 Callable Libraries.

The sequential solution approach can be summarized as follows:

Step 1:

List the different destinations according to a specific criterion.

Step 2:

Calculate D_it=d_it+_wWd_i(t+^w)^w, iFP; +1tT+.

Step 3:

Determine the production plan at the plant by solving model (P') (using cuts (Cut.3) and (Cut.4) and Cplex 9.0 libraries).

Step 4:

Following the sequence defined in Step 1, determine the replenishment plan for each destination wW by solving model (D-w) (using cuts (Cut.1) and cuts (Cut.2) or (Cut.2') and Cplex 9.0 libraries).

The sequence in which the different destinations are listed, in Step 1 of the sequential solution approach proposed, may have an impact on the overall solution quality. A destination in a given position of the sequence may have more (less) inventory available at the plant than its successor (predecessor) to take advantage of transportation economies of scale. In our experimental test section, the quality of the solutions obtained is assessed when the different destinations are listed in Step 1 based on the following criteria:

• Destinations total volume (DTV): The destinations are listed in decreasing order of their total demand for the whole planning horizon:
  
  
  
  
• Total transportation cost (TTC): The destinations are listed in decreasing order of their lowest TTC. For a given destination, the TTC is given by:
  
  
  
  
• Transportation cost per unit (TCPU): The destinations are listed in decreasing order of their unit transportation cost. For a given destination, the TCPU is given by:
  
  
  
  Also, we noted that when solving model (D-w) in Step 4, Cplex finds optimal and/or quality solutions relatively quickly but spends the majority of its computation time proving the optimality of these solutions. Moreover, the objective of our sequential approach is not to find optimal solutions for models (D-w) but rather find a good solution for model (F) within an acceptable time. To reduce computation times, this suggest that the MIP gap to use while solving models (D-w) should not be too small. This observation is further demonstrated in our experimental results section.

In addition to providing a heuristic to solve problem (F), the solution approach proposed in this section gives the best solution that can be obtained when the production and distribution problems are solved sequentially, as is usually done in practice. By comparing the optimal solution of (F) with the solution provided by this approach, it is therefore possible to asses the gain of integrated production–distribution planning over independent and synchronized planning.

4. Experimental testing

4.1. Test application

The models proposed in this paper were developed in the context of a project with a large pulp and paper company (the second largest producer of uncoated free-sheet in North America and the third in the world), and our numerical tests and comparisons are based on real data. Paper production starts with pulp mills transforming wood logs into a mix of pulp known as furnish (cellulose fibre). Paper machines then convert the furnish to a paper sheet by spraying it into a wide moving mesh. From this point, water is removed from the paper sheet in the press section and high-pressure rollers smooths its texture. The paper sheet is rolled around solid cores and then cut to the width and diameter required by the customers or by the converting operations required to produce the FPs, such as offset and photocopying paper sizes, sold on the market. Transportation to different warehouses and customers is performed using full truck loads (TL).

The mill considered in the study has two paper machines. The two machines use different types of pulp to produce three different grades of paper. Business grade paper, as its name indicates, is used in business offices and for personal computers. Printing and publishing grade papers are used for publications and glossy publications. Specialty and technical grade papers are highly specialized papers mainly used in flexible packaging and in industrial applications. Every time a paper machine switches its production to a different grade, there are losses of paper sheets that vary depending on which grade was produced before (ie changeover costs are sequence dependent). Changeover costs are estimated based on the paper sheet losses. The plant ships FPs to distribution centres, as well as, directly to customer locations. It is assumed that inventory costs are the same for all FPs and are time independent in both the plant and the warehouse locations.

4.2. Experimental design

The aim of our experiments is to measure the quality of the different strategies proposed in Section 4 to solve model (F). Using the pulp and paper company demand, capacity and cost data, different problem instances were generated as listed in Table 5. Column headings are: the name of the problem instance (Instance), the number of planning periods (T), the number of potential transportation modes to each destination (Modes Number), the number of destinations (Destinations Number), the number of FPs considered, the number of continuous variables (Cont. Var.), the number of 0-1 variables (0-1 Var.) and the number of constraints (Constraints) in the model (without additional cuts). In problem instances Pb1 and Pb2, products are shipped from the plant to the different destinations using TL only. In addition to truck, shipments to the different destinations can also be made via rail in problem instances Pb3 and Pb4. To measure the impact of including, in our model, a slower but cheaper transportation mode on the different cost components, costs as well as computation times of Pb3 and Pb4 are compared to Pb1 and Pb2, respectively. The callable libraries of CPLEX version 9.0 were used to solve the different instances of our test problems. A time limit of 2 h (7200 s) was imposed on the solution time of all runs. All experiments were performed on a 1.66 GHz Pentium III workstation with 1 GB main memory running under Windows NT.

Table 5 - Test instances.

Full table

4.3. Experimental results

To measure the impact of using SOI (Section 3.2), as well as, the impact of adding special 0-1 variables (Section 3.3), we run tests to solve problem instance Pb1 using the four different formulations (F), (F.Bin), (F-SOI) and (F-SOI.Bin). Adding valid inequalities to a problem formulation may have a positive or a negative impact on its solution time. Therefore, we tested the impact of different combinations of the valid inequalities proposed in Section 3.1 on the computation time of solving the different formulations (F), (F.Bin), (F-SOI) and (F-SOI.Bin). CPLEX MIP solver also adds general cuts to a problem formulation while solving it. To isolate the impact of the valid cuts proposed in this paper from those of CPLEX, two sets of experiments were performed. In the first set, CPLEX cuts were allowed to be generated, while in the second set they were not allowed to be generated.

Results showing the impact on computation times of adding different combinations of valid cuts to models (F), (F.Bin), (F-SOI) and (F-SOI.Bin) are presented in Tables 6 and 7. These tables report the average elapsed CPU times, in seconds, without CPLEX cuts (w/o-CPLEX cuts) and with CPLEX cuts (w-CPLEX cuts). These results show that major improvements can be achieved by adding some valid cuts to the different formulations. In addition, it is clear that, out of the five valid cuts proposed in Section 4.1, (Cut.2) are the most effective in improving the different formulations solution time. Adding (Cut.2), (Cut.3) and (Cut.4) yields the best results for (F), (F.Bin), (F-SOI) and (F-SOI.Bin). The results also suggest that allowing or not allowing CPLEX to generate its own cuts does not have a major impact on the performance of our valid cuts. When used with valid inequalities (Cut.2), (Cut.3) and (Cut.4), formulation (F.Bin) proved to be the most efficient. Formulation (F.Bin) along with valid inequalities (Cut.2), (Cut.3) and (Cut.4) are then used in what follows to solve problem instances Pb2, Pb3 and Pb4.

Table 6 - Impact of valid cuts and special 0-1 variables on (F) with and without CPLEX cuts.

Full table

Table 7 - Impact of valid cuts and special 0-1 variables on (F-SOI) with and without CPLEX cuts.

Full table

Table 8 compares the computation time, the different cost components, as well as, the total cost of problem instances Pb3 and Pb4 solutions with those of problem instances Pb1 and Pb2, respectively. Column headings of Table 8 are: the name of the problem instance (Instance), production cost deviation percentage (% ProdC dev.), plant inventory carrying cost deviation percentage (% PinvcC dev.), destinations inventory carrying cost deviation percentage (% DinvcC dev.), transportation cost deviation percentage (% TransC dev.) and total cost deviation percentage (% TotalC dev.). Table 8 shows that by including rail in the model as a competing transportation mode to TL, our problem instances total costs decreased, in average, by 17%. It is interesting to note that even though the costs of production and the inventory carrying costs at the destinations slightly increased, major gains in transportation costs were achieved. Reduction of the plant inventory carrying costs were also observed. As should be expected, considering a slower but cheaper transportation mode (rail), the model solution increased its production changeovers so as to increase its utilization of the cheaper transportation mode. This result is particularly interesting because for the real case considered, all transportation is currently made by truck.

Table 8 - Impact of adding rail as an additional transportation mode.

Full table

To test the quality and the efficiency of our sequential solution approach, we first performed a comparison of the computation times and the solutions quality of the approach using the different destinations sorting criteria proposed in Section 3.2.2. The three different sorting criteria DTV, TTC and TCPU were respectively used in Step 1 of the sequential solution approach to solve problem instances Pb1 and Pb2. The average computation time as well as the total cost comparison of the solutions obtained are presented in Table 9. The results show that different sorting of the destinations did not have a major impact on the solution time and quality of our heuristic approach. The DTV criteria yield a slightly better solution quality than the others. The DTV criteria was consequently used in Step 1 of our approach to solve problem instances Pb3 and Pb4.

Table 9 - Comparison of the sequential solution approach results using different criteria to prioritize the destinations in Step 1.

Full table

The computation time and quality of the solutions obtained with our sequential solution approach are compared, in Table 10, to those of optimal and/or best solutions obtained using model (F.Bin) along with valid inequalities (Cut.2), (Cut.3) and (Cut.4). Note that the 'MIP gap' column in Table 10 indicates Cplex 9.0 MIP gap used to solve models (D-w) in Step 4 of the sequential solution approach. Problem instances Pb1 and Pb2 are solved using a MIP gap of 1 and 0.01% to show the impact of using a relatively higher MIP gas on both the computation time and the solution quality. For each problem instance, Table 10 shows the percentage deviation in computation time (% Time dev.), in total cost (% TotalC dev.), in production cost (% ProdC dev.), inventory cost at the plant (% PinvC dev.), in inventory cost at the destinations (% DinvC dev.) and in transportation cost (% TransC dev.) of the heuristic approach solution relatively to the one obtained by solving (F.Bin) with valid inequalities (Cut.2), (Cut.3) and (Cut.4). Table 10 shows that instead of a MIP gap of 0.01%, using a MIP gap of 1% in Step 4 of the sequential solution approach, resulted in a computation time reduction of 14, 39, 88 and 97% for problem instances Pb1, Pb2, Pb3 and Pb4 respectively. On the other side, the solution total cost increased only by 0.23% for Pb1, 0.08% for PB2, 0.20% for Pb3 and 0.26% for Pb6. Table 10 also shows that the sequential solution approach is very effective in solving problem instances with a single transportation mode (Pb1 and Pb2). In this case, it was able to find quality solutions (within 1% of the optimal solutions) in a computation time that represents less than 15% of the time required to find optimal solutions of model (F.Bin) with Cplex 9.0 and valid inequalities (Cut.2), (Cut.3) and (Cut.4).

Table 10 - Sequential solution approach computation time and solution quality.

Full table

As can be observed, however, for problem instances Pb3 and Pb4, the sequential solution approach only managed to find solutions with total costs that are within 15% of those of the integrated approach solutions obtained by solving model (F.Bin) with Cplex 9.0 and valid inequalities (Cut.2), (Cut.3) and (Cut.4). For these problems, the approach fails to exploit the potential savings offered by the multiple transportation modes available. One of the factors that contributes to this shortcoming is the way the network total effective demand (D_it=d_it+_w=1^Wd_i(t+^w)^w) is computed for the production problem (P'). The total effective demand (D_it) is computed using the shortest transportation lead time to each destination (^w) to ensure the production problem (P') feasibility. However, it falls short of ensuring that large quantities of products can be available ahead of time to be able to take advantage of cheaper transportation modes, that usually have longer transportation lead times. On the other side, using a longer transportation lead time to compute (D_it) will force the plant to produce products much earlier than they are required and thus may result in an infeasible production problem or the plant will be forced to carry high levels of inventory as is the case usually in practice. This shows that, it may be very beneficial, in some practical situations, to try to solve the production–distribution planning problem with an integrated approach, instead of solving the production and distribution problems independently, as is often done. Note also, that theapproach proposed here to solve the production and distribution problems sequentially is designed to coordinate the solution of the sub-problems. This is why the solutions obtained for Pb1 and Pb2 are very good. If such a coordination and synchronization scheme had not been implemented, the difference between the quality of the integrated approach solution and the independent approach solution would have been even greater.

5. Conclusion

In this paper, we examined the production–distribution planning problem of multiple items between a manufacturing location and multiple destinations. The manufacturing location is assumed to operate a serial production process with a bottleneck stage, subject to a predetermined production sequence. The transportation costs between the plant and the different destinations are characterized by piecewise linear economies of scale. In addition, we consider that competing transportation modes can be available to make shipments between the plant and the different destinations. As is the case in practice between rail and truck for example, these transportation modes may have different transportation lead times and different cost structures. To model this problem, a tight mixed integer programming formulation is proposed, and valid inequalities to 'strengthen' the formulation are derived. Furthermore, we proposed the addition of special binary variables that represent 'dichotomies' in the mathematical formulation. The purpose of adding these special binary variables is to increase the branching efficiency of the branch-and-bound algorithm used by most commercial solvers. We also show how our piecewise linear transportation cost functions can be modelled using an alternative formulation known in the literature as SOI. Finally, we proposed a sequential solution approach to solve large versions of the problem and to evaluate the potential impact of integrated planning.

Our experimental results showed that in the case of a single transportation mode, adding our valid inequalities to (F) improved Cplex 9.0 solving times by a factor of 6. Moreover, adding special 0-1 variables that represent 'dichotomies' reduced computation times even further. However, using the SOI proposed did not seem to have a significant impact. Also, in the case of a single transportation mode, our sequential solution approach was able to find solutions that are within 1% of the optimal solution with computation times that are, in average, more than seven times faster than solving (F.Bin) using Cplex 9.0 and valid inequalities (Cut.2), (Cut.3) and (Cut.4).

On the other hand, it was shown that a sequential solution approach fails to capture the potential savings offered by the use of multiple transportation modes. Given the magnitude of these potential savings, it is clear that it would be rewarding to invest in research to find heuristic methods that are capable of solving model (F) more efficiently. We believe in particular that heuristics based on Lagrangean relaxation (see Rizk et al 2006a) and metaheuristics are worth investigating.

References

1. Anily S (1994). The general multi-retailer EOQ problem with vehicle routing costs. Eur J Opl Res 79: 451–473. | Article |
2. Anily S and Federgruen A (1990). One warehouse multiple retailers systems with vehicle routing costs. Mngt Sci 36: 92–114.
3. Anily S and Federgruen A (1993). Two-echelon distribution systems with vehicle routing costs and central inventories. Opns Res 41: 37–47.
4. Barany I, Van Roy TJ and Wolsey LA (1984). Strong formulations for multi-items capacitated lotsizing. Mngt Sci 30: 1255–1261.
5. Beale EML and Tomlin JA (1970). Special facilities in general mathematical programming system for non-convex problems using ordered sets of variables. In: Lawrence J. (Ed). Proceedings of the Fifth International Conference on Operational Research. Tavistock Publications: London, pp. 447–454.
6. Chan LMA, Muriel A, Shen ZJ, Simchi-Levi D and Teo C (2002). Effective zero-inventory-ordering policies for the single-warehouse multiretailer problem with piecewise linear cost structures. Mngt Sci 48: 1446–1460. | Article |
7. Chandra P and Fisher ML (1994). Coordination of production and distribution planning. J Opl Res Soc 72: 503–517.
8. Diaby M and Martel A (1993). Dynamic lot sizing for multi-echelon distribution systems with purchasing and transportation price discounts. Opns Res 41: 48–59.
9. Gallego G and Simchi-Levi D (1990). On the effectiveness of direct shipping strategy for the one warehouse multi-retailer R-systems. Mngt Sci 36: 240–243.
10. Herer Y and Roundy R (1997). Heuristics for a one-warehouse multiretailer distribution problem with performance bounds. Opns Res 45: 102–115.
11. Orgyczak W (1996). A note on modeling multiple choice requirements for simple mixed integer programming solvers. Comput Opns Res 23: 199–205. | Article |
12. Rizk N, Martel A and Ramudhin A (2006a). A Lagrangean relaxation algorithm for multi-item lot sizing problems with joint piecewise linear resource costs. Int J Prod Econom 106: 344–357. | Article |
13. Rizk N, Martel A and D'Amours S (2006b). Multi-item dynamic production distribution planning in process industries with divergent finishing stages. Comput Opns Res 33: 3600–3623. | Article |
14. Salmon M (eds) (1991). Multi-stage Production Planning and Inventory Control, Lectures Notes in Economics and Mathematical Systems, vol. 355. Springer-Verlag: Berlin, pp 92–108.
15. Viswanathan S and Mathur K (1997). Integrating routing and inventory decisions in one-warehouse multiretailer multiproduct distribution systems. Mngt Sci 43: 294–312.

Acknowledgements

This project would not have been possible without the collaboration of FOR@C's partners especially of NSERC and Domtar.