A primal decomposition method for the integrated design of multi-period...

KORAY DOGAN [1]

MARC GOETSCHALCKX [2]

We study the integrated design of strategic supply chain networks and the determination of tactical production-distribution allocations in the case of customer demands with seasonal variations. Given a set of potential suppliers, potential

manufacturing facilities and distribution centers with multiple possible configurations, and customers with seasonal demands, the goal is to determine the configuration of the production-distribution system with the lowest sum of supply, production, transportation, inventory, and facility costs such that seasonal customer demands are met. We develop a mixed integer programming formulation and an integrated design methodology based on primal (Benders) decomposition. For a case study in the packaging industry, specialized acceleration techniques reduced the running times by a factor of 480. The company projects savings of 2% or $8.3 million by using the integrated rather than the optimal hierarchical configuration.

1. Integrated production--distribution systems

1.1. Problem definition

A supply chain is a network of suppliers, production facilities, and distribution facilities that transforms raw materials into finished products through multiple manufacturing stages and delivers them to the customers through one or more distribution stages. Strategic supply chain decisions typically involve a time horizon of more than 1 year and determine the location of facilities, configuration of the manufacturing and distribution facilities, and selection of distribution channels. We also consider the establishment and size of individual production lines part of the strategic decisions, because of the cost and time involved in acquiring or moving such a production line. Tactical decisions typically deal with a time horizon of 3 months to 1 year and determine the production mixes and levels at the various facilities, the inventory levels in the facilities, the allocation of customers, and the transportation activities.

Without loss of generality, we considered a manufacturing process that consists of a primary and secondary production phase. For each phase one or more alternative manufacturing lines exist that differ by their technology or capacity. The manufacturing lines have maximum production capacities and the resource requirement and the constant marginal cost for manufacturing a particular product on a particular production line are known. The two-stage production--distribution network is illustrated in Fig. 1.

The different manufacturing and distribution facilities are connected geographically by transportation channels and temporally by seasonal inventory. At all facility sites cycle inventories are created by the incoming and outgoing shipment frequencies. Selection of the transportation channels between the facilities is driven by the transportation costs and the cycle and pipeline inventory costs. The transportation channels have a throughput capacity.

The only customer service constraints considered are the full satisfaction of the seasonal demands of the customers. We assumed that customer demand is deterministic and ignored the impact of safety stocks.

1.2. Literature review

The design of production and distribution systems has been active area of research during the last 30 years. There exist a large body of research on the modeling and design of various components of the integrated multi-period production--distribution problem. Inclusion of all the papers that focus on a single component of the system is impossible. Several comprehensive reviews have been published in the last few years. Geoffrion and Powers [1] discuss the evolution of strategic distribution systems design over the last 20 years. They observe that the associated large-scale models have proven to be extremely difficult to solve to optimality without the application of Benders decomposition or factorization methods. They report however that Benders decomposition is now being removed from the optimizers because of the level of technical resources and implementation effort that it requires to perform well. Thomas and Griffin [2] present a comprehensive review of the literature related to the coordination of two or more of the main stages of the supply chain, that is, procurement, production, and distribution. Vidal and Goetschalckx [3] provide a comprehensive review of models in the literature that have a special emphasis on global logistics systems.

A variety of solution algorithms for the capacitated facility location problem has been proposed. An early review is given in Aikens [4]. The capacitated facility location problem exhibits a related, but significantly simpler, structure to the problem discussed in this paper. The most efficient solution reported to date is based on cross decomposition integrating primal (Benders) decomposition with dual decomposition (Langrangean relaxation), Van Roy [5,6].

In an early seminal paper, Geoffrion and Graves [7] present an algorithm based on Benders decomposition to solve a multicommodity single-period production-distribution problem. The model uses a path-based formulation tracing products from suppliers to customers through one intermediate facility. The model also incorporates single sourcing constraints, a feature required by the solution procedure. Cohen and Lee [8,9] propose a modeling framework for general production-distribution systems. The resulting model has a unified, hierarchical and stochastic network structure. This overall model could not be solved to optimality. The authors develop a hierarchical approach and solve several of the subproblems. Arntzen et al. [10] develop a model for a global multi-period supply chain system, incorporating bill of material constraints. The model contains several international features such as duty relief and drawback. The authors report the solution of several real life cases by employing non-traditional solution met hods, such as elastic constraints, row factorization, cascaded problem solution, and constraint-branching enumeration. In a series of articles Lee [11,121 implements primal and cross decomposition algorithms for modeling multiple facility types. The models impose capacity constraints by commodity and lack joint capacity constraints. Goetschalckx et al. [13] develop a comprehensive model for a multicommodity, multi-facility, multi-facility-type, and multiechelon distribution system. The arc-based model is solved by a general purpose mixed integer programming solver using preprocessing and partial implementation of strong formulation constraints.

The facility location problem has been a popular research area due to its mathematically challenging nature and relation to the real life strategic distribution design problem. The most efficient methods for solving it are based on cross decomposition, which exploits the primal and/or dual structure of the model. However, the research on integrated multi-period production-distribution problems is much more limited. A few integrated models have been proposed. Their computational complexity and size have forced researchers to use hierarchical decomposition and/or heuristic algorithms to obtain solutions. We will first develop a model for the integrated strategic-tactical design and then present a primal decomposition algorithm with specialized acceleration techniques that yields the optimal solution in a reasonable amount of time. We will then illustrate the value of the integrated versus the hierarchical approach based on a case study in the packaging industry.

2. Model development

2.1. Network flow representation

The integrated production-distribution network design problem can be represented by a fixed charge and multicommodity network flow problem. The network has four regions. The first region corresponds to the suppliers or mills in our case study. The second region corresponds to the first phase of the manufacturing process and the third region corresponds to the second phase of the manufacturing process. Finally, the fourth region corresponds to the customers. Transportation arcs connect the regions. A schematic of the network structure for a single period is shown in Fig. 2.

The arcs between a pair of nodes in the supplier region allows us to model joint product capacity constraints in each of the mills. Each facility in the primary and finishing region is represented by a sequence of four nodes. The arcs between the first and second node allow us to model fixed and variable cost and capacity constraints for the production section of the facility. The arcs between the second and third node allow us to model fixed and variable costs and capacity constraints for individual production lines in the facility. Finally, the arcs between the third and fourth node allow us to model fixed and variable costs and capacity constraints for the warehousing section of the facility.

Since this is a multi-period problem, the complete network consists of a number of identical planes, one dedicated to each period. The only connections between the parallel planes of the network are the arcs between the nodes associated with the same warehouse in subsequent periods.

The overall production-distribution network design problem can be decomposed into two subproblems. The first subproblem determines the resource location and sizing and the production allocation. The second subproblem determines the tactical transportation flows.

2.2. Strategic resource sizing and production allocation problem

The arcs in this problem correspond to the production and warehousing facilities and the production lines in the manufacturing facilities. The arcs have a fixed and variable costs and a maximum throughput capacity. The binary variables correspond to the open or closed status of each manufacturing and warehouse facility and of the manufacturing lines in the facilities. The continuous variables correspond to the amount of each product purchased from the various suppliers, the amount manufactured on the various production lines, and the seasonal inventory levels at the end of every period.

To ensure that the solution of any resource sizing and allocation formulation has a feasible transportation solution three constraints are added. They ensure that for every combination of product and time period the total quantity of raw materials, the total throughput of the first stage manufacturing process, and the total throughput of the second stage manufacturing process, respectively, are at least equal to the total demand for that product during that period. These constraints will later eliminate the need to consider extreme rays in the primal decomposition method. The resource sizing and allocation problem is solved with the MIP module of CPLEX, using its linear programming relaxation to provide bounds.

2.3. Multicommodity network flow problem

The arcs in this problem correspond to the various transportation channels. The arcs have variable costs and joint throughput capacities. The continuous variables correspond to transportation quantities. Parallel arcs between the same origin and destination facility can exist to model transportation channels with different shipment frequencies and/or transportation times. The frequency then determines the amount of cycle inventory at the origin and destination facility.

All products must undergo the same two manufacturing processes, so arcs from suppliers to the second manufacturing stage or from the first manufacturing stage to the customers are not allowed. This implies that the single period network decomposes into three networks: from suppliers to first stage, from first stage to second stage, and from second stage to customers, respectively.

The multicommodity network flow problems are solved with the network-based linear programming module of CPLEX.

3. Formulation

3.1. Verbal model

Minimize:

Total Cost = Supply cost + Fixed manufacturing cost

+ Variable manufacturing cost

+ Fixed facility operating cost

+ Variable facility operating cost

+ Warehousing cost

+ Cycle inventory cost at the facilities

+ Pipeline inventory cost

+ Inventory carry-over cost

+ Transportation cost.

Subject to: Customer demand satisfaction,

Conservation of flow at facilities,

Conservation of flow at suppliers,

Conservation of flow at machines,

Supplier capacity, Facility capacity,

Machine capacity, Single facility type at a site,

Linkage constraints between machines and facilities.

3.2 Notation

Sets and indices

T = set of time periods;

t = index for time periods;

P = set of products;

p = index for products;

S = set of suppliers;

s = index for suppliers;

I = set of production facilities;

I = index for production facilities;

J = set of finishing facilities;

W = set of production-warehousing facilities;

X = set of finishing-warehousing facilities;

j = index for finishing-warehousing facilities;

[L.sub.i] = set of available facility types for facility i;

l = index for facility types;

K = set of customers;

k = index for customers;

[M.sub.i] = set of available machines at facility i;

[M.sub.I] = set of available machines at production facilities;

[M.sub.J] = set of available machines at finishing facilities;

m = index for machines.

Parameters

[SCap.sub.ts] = resource capacity of supplier s during period t;

[SReq.sub.tps] = supply resource requirement for product p by supplier s during period t;

[SCost.sub.tps] = supply cost of product p by supplier s during period t;

[Fcap.sub.til] = resource capacity for facility i of type l during period t;

[FReq.sub.tpil] = resource requirement of product p for facility type l of facility i during period t;

[FCost.sub.tpil] = processing cost of product p at facility type l of facility i during period t;

[FFxCost.sub.il] = fixed cost of facility type l of facility j;

[WCap.sub.til] = resource capacity for warehouse w of type l during period t;

[WReq.sub.tpil] = resource requirement of product p for facility type l of warehouse w during period t;

[WCost.sub.tpil] = handling cost of product p at facility type l of warehouse w during period t;

[WFxCost.sub.til] = fixed cost of facility type l of warehouse m during period t;

[MCap.sub.tim] = resource capacity for machine m at facility i during period t;

[MReq.sub.tpim] = resource requirement of product p for machine m at facility i during period t;

[MCost.sub.tpim] = processing cost of product p at machine m at facility i during period t;

[MFxCost.sub.im] = fixed cost of machine m at facility i;

[V.sub.p] = value of product p;

[D.sub.tpk] = demand of customer k for product p during period t;

[RF.sub.tij] = replenishment frequency for transportation channel ij during period t;

[TT.sub.tij] = transit time of channel ij during period t;

[TCost.sub.tpij] = transportation cost for product p through channel ij during period t;

r = inventory carrying cost rate per time period;

R = inventory carrying cost rate for planning horizon;

[MaxN.sub.jm] = maximum number of machines of type m at facility j.

Decision variables

[X.sub.tpij] = flow of product p shipped through transportation channel ij during period t;

[Z.sub.tpjm] = production of product p at machine m of facility j during time period t;

[S.sub.tps] = supply of raw material for product p from supplier s during time period t;

[n.sub.jm] = the number of machines of type m operating at facility j;

[q.sub.tpj] = the inventory level of product p at warehouse j at the end of period t;

[v.sub.tpjl] = flow of product p through facility j of facility type l during period t;

[y.sub.jl] = {1 if facility type l is selected open at facility l, 0 otherwise.

3.3. Model development

We will focus here on the features and assumptions of the model that are different from the model in Goetschalckx et al. [13]. Further details can be found in Dogan and Goetschalckx [14] and a complete description of the model can be found in Dogan [15].

3.3.1. Objective function cost components

Raw materials are assumed to be introduced into the supply chain system at the suppliers. The supply cost may vary at different suppliers at different periods. The total manufacturing cost of the system has two parts: (i) the fixed cost of having the machines operational; and (ii) the variable cost of the production itself. The variable manufacturing cost is modeled at the strategic level and operational issues, such as machine scheduling and setup costs, are estimated in in the cost and resource parameters. The cost of manufacturing product p at machine m of facility i is assumed to be proportional to the amount of production.

The inventory cost includes the cost of holding the work-in-process inventory and the finished goods at the facilities and in the transportation channels. The inventory costs at the suppliers are not considered in this model. The actual level of inventory at the facilities depends on the order quantities of the incoming and the outgoing shipments, phasing of these shipments and the operational issues. This makes modeling them at the strategic level a difficult task. To keep the model simple, the replenishment frequencies at the incoming transportation channels are assumed to be constant over the planning horizon and outgoing shipments are assumed to be continuous. In addition, the effect of inventory phasing is ignored. The cumulative inventory level is obtained by adding the inventory levels of individual shipments, where

[Q.sub.s] = [x.sub.tpsj]/[RF.sub.tsj].

The average inventory associated with regular shipments from facility is [Q.sub.s]/2 and the cost of holding inventory at the facility over a single period is calculated by (r[V.sub.p][Q.sub.s])/2. The total inventory holding cost at facilities is given by:

[[sigma].sub.t[epsilon]T] [[sigma].sub.p[epsilon]P [[sigma].sub.i[epsilon]{[S.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][W.sup.t]}] [[sigma].sub.j[epsilon]{[I.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][J.sup.t]}] r[V.sub.p]/2[RF.sub.tij] [x.sub.tpij].

Another major component of the inventory cost is the pipeline inventory cost, especially on the long trunking transportation channels.

Inventory cost is also charged for the inventory on hand at the end of each period. For ease of notation, the inventory carry-over cost is only considered at the warehouses between the production facilities. The inventory level of product p at warehouse j at the end of time period t is equal to the difference between the total outflow and the total inflow of the warehouse during time period t. The multi-season model is considered as cyclic, i.e., the inventory levels at the end of last time period wrap around as the beginning inventory levels of the first time period. Figure 3 illustrates the ending inventory levels for product p at warehouse j. Let [q.sub.1], [q.sub.2], and [q.sub.3] be the ending inventory levels of time period 1, 2, and 3, respectively. Observe that for the optimal solution, at least one of the inventory levels at the end of a period must be equal to zero. Using a linear approximation of the inventory levels, the area below the inventory line for each period is [q.sub.1] + [q.sub.2] + [q.sub.3].

The total inventory cost of carry-over items is given by:

[[sigma].sub.t[epsilon]T] [[sigma].sub.p[epsilon]P] [[sigma].sub.j[epsilon][X.sup.t]] [RV.sub.p][q.sub.tpj].

Each production facility has a location and can have multiple types with a different costs and capacity associated with each type. Different types can be used to model different sizes or different technologies. The operating cost of a facility-type is modeled with fixed and variable costs. The fixed portion of a facility-type operating cost is the cost of having the production or warehousing function of the facility-type in operation and is independent of the throughput of the facility. The system incurs variable cost for the throughput at the manufacturing facilities. These costs represent the expenses that are directly related to the throughput of the facility, such as the material handling, storage, and hourly labor costs. The operating cost rate depends on the type of the facility operational at the site. The transportation cost is assumed to be linearly dependent on the magnitude of the flow.

3.3.2. Constraints

The model has the standard customer demand satisfaction constraints and conservation of flow constraints. Figure 4 illustrates the conservation of flow constraints at the facilities.

The production and warehousing facility nodes are disaggregated into two nodes each to represent the flow variables as arcs. Parallel arcs between nodes (l)-(2) and (3)-(4) represent the flow through alternative facility types, and the parallel arcs between (2)-(3) represent the flow on different machines.

At the outflow nodes of production facilities, the throughput of the facility must be equal to the production at the machines (Fig. 4, nodes (2)):

[[sigma].sub.m[epsilon][[M.sup.t].sub.j]] [Z.sub.tpjm] - [[sigma].sub.l[epsilon][L.sub.j]] [V.sub.tpjl] = 0 t [epsilon] T, j [epsilon] {[I.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][J.sup.t]}, p [epsilon] P.

At the inflow nodes of warehousing facilities, the throughput of the facility plus the ending inventory must be equal to the production at the machines plus the beginning inventory (Fig. 4, nodes (3)):

[q.sub.tpj] + [[sigma].sub.l[epsilon][L.sub.j]] [V.sub.tpjl] - [[sigma].sub.m[epsilon][[M.sup.t].sub.j]] [Z.sub.tpjm] -[q.sub.t]-[1.sub.pj] = 0

t [epsilon] T, i [epsilon] {[W.sup.t] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sup.t]}, P [epsilon] P.

Capacity is the upper bound on the total resources that can be consumed by the throughput of the products. The capacity of a production facility is expressed at two levels: machine level and facility level. Total resource consumption is bounded by the available uptime of the machine. At the facility level, the resource consumption of the total production of the facility type is bounded by its capacity. Facility types are mutually exclusive design alternatives for the facility; hence at most one facility type is in use at a potential site. The model includes the standard singular facility type constraints and linkage constraints. Special linear side constraints can be added to represent joint decisions.

4. Primal decomposition method

4.1. Implementation of primal decomposition

The structure of the tactical production--distribution problem leads naturally to a primal decomposition scheme. At each iteration, the mixed integer master problem is solved to obtain a solution that specifies the status of the facilities, the production lines, and the production and inventory quantities. The transportation subproblem is then solved to obtain the transportation flows with the lowest costs. This transportation solution is added to the master problem as an additional constraint, and the mixed integer master problem is forced to search for a solution with a lower total cost. The algorithm terminates when the master problem cannot find a solution with a lower total cost; i.e., the master problem has become infeasible.

The decomposition principle is implemented by selecting [nu], z, y, and n variables as the complicating variables. When these variables are fixed, the constraints containing only these variables drop from the formulation, and the resulting mathematical program is expressed by the following form, LP(v,z,y,n).

Min [[sigma].sub.t[epsilon]T] [[sigma].sub.p[epsilon]P] [[sigma].sub.i[epsilon]{[S.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][W.sup.t]}] [[sigma].sub.j[epsilon]{[I.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][J.sup.t]}] r[V.su.p]/2[RF.sub.tij] [x.sub.tpij]

+ [[sigma].sub.t[epsilon]T] [[sigma].sub.p[epsilon]P] [[sigma].sub.i[epsilon]{[S.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][W.sup.t]}] [[sigma].sub.j[epsilon]{[I.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][J.sup.t]}] r[V.sub.p][T.sub.tij][x.sub.tpij]

+ [[sigma].sub.t[epsilon]T] [[sigma].sub.p[epsilon]P] [[sigma].sub.i[epsilon]{[S.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][W.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][W.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][X.sup.t]}] [[sigma].sub.j[epsilon][{[K.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][I.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][J.sup.t]}.sub.j]] [TCost.sub.tpij][v.sub.tpij].

s.t.

[[sigma].sub.j[epsilon][X.sup.t]] [x.sub.tpjk] = [D.sub.tpk] t [epsilon] T, p [epsilon] P, k [epsilon] K [[zeta]],

- [[sigma].sub.j[epsilon][K.sup.t]] [x.sub.tpij] = - [[sigma].sub.l[epsilon][L.sub.j]] [v.sub.tpjl] t [epsilon] T, j [epsilon] {[I.sup.t] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [J.sup.t]}, p [epsilon] P [[beta],[delta]],

[[sigma].sub.j[epsilon]{[S.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][W.sup.t]}] [x.sub.tpij] = [[sigma].sub.l[epsilon][L.sub.j]] [v.sub.tpjl] t [epsilon] T, i [epsilon] {[W.sup.t] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sup.t]}, p [epsilon] P [[gamma],[epsilon]],

[[sigma].sub.i[epsilon][I.sup.t] [x.sub.tpsi] = [S.sub.tps] t [epsilon] T, s [epsilon] [S.sup.t], p [epsilon] P [[alpha]],

This model is a multicommodity network flow problem and LP(v,z,y,n) can be decomposed into \P\ separate network flow problems. Let the dual variables associated with the above constraints be [alpha], [beta], [delta], [gamma], [epsilon] and [zeta] for the nodes in the sets, S, I, W, J, X, and K, respectively. These dual variables represent the relative unit transportation cost of supplying the demand at the associated nodes.

Solving the integrated production--distribution model is equivalent to solving the following mixed integer program, (IPDM).

Min [[sigma].sub.t[epsilon]T] [[sigma].sub.p[epsilon]P] [[sigma].sub.s[epsilon]S] [SCost.sub.tps][S.sub.tps] + [[sigma].sub.j[epsilon]{I[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]J}] [[sigma].sub.m[epsilon][M.sub.j] [MF.sub.x][Cost.sub.jm][n.sub.jm]

+ [[sigma].sub.t[epsilon]T] [[sigma].sub.p[epsilon]P] [[sigma].sub.j[epsilon]{[I.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][J.sup.t]}] [[sigma].sub.m[epsilon][[M.sup.t].sub.j]] [MCost.sub.tpjm][Z.sub.tpjm] + [[sigma].sub.t[epsilon]T] [[sigma].sub.p[epsilon]P] [[sigma].sub.j[epsilon][X.sup.t]] RV p[q.sub.tpj]

+ [[sigma].sub.j[epsilon]{[I.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][J.sup.t]}] [[sigma].sub.l[epsilon][L.sub.j]] [FF.sub.x][Cost.sub.jl][y.sub.jl] + [[sigma].sub.j[epsilon]{[W.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][X.sup.t]}] [[sigma].sub.l[epsilon][L.sub.j]] [WF.sub.x][Cost.sub.jl][y.sub.jl]

+ [[sigma].sub.t[epsilon]T] [[sigma].sub.p[epsilon]P] [[sigma].sub.j[epsilon]{[I.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][J.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][W.sup.t][MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII][X.sup.t]}] [[sigma].sub.l[epsilon][L.sub.j]] [FCost.sub.tpjl][v.sub.tpjl] + [Z.sub.0].

s.t.

[q.sub.tpj] + [[sigma].sub.l[epsilon][L.sub.j]] [v.sub.tpjl] - [[sigma].sub.m[epsilon][[M.sup.t].sub.j]] [Z.sub.tpjm] - [q.sub.t]-[1.sub.pj] = 0

t [epsilon] T, i [epsilon] {[W.sup.t] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sup.t]}, p [epsilon] P,

[[sigma].sub.m[epsilon][[M.sup.t].sub.j]] [Z.sub.tpjm] - [[sigma].sub.l[epsilon][L.sub.j]] [V.sub.tpjl] = 0 t [epsilon] T, j [epsilon] {[I.sup.t] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [J.sup.t]}, p [epsilon] P,

[[sigma].sub.p[epsilon]P] [SRep.sub.tps][S.sub.tps] - [SCap.sub.ts] [less than or equal to] 0 s [epsilon] [S.sup.t], t [epsilon] T,

[[sigma].sub.p[epsilon]P] [MReq.sub.tpi][Z.sub.tpim] - [MCap.sub.tim][n.sub.im] [less than or equal to] 0

t [epsilon] T, i [epsilon] {[I.sup.t] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [J.sup.t]}, m [epsilon] [[M.sup.t].sub.i],

[[sigma].sub.[p.sup.[epsilon]P]] [FReq.sub.tpjl][V.sub.tpjl] - [FCaP.sub.tjl][Y.sub.jl] [less than or equal to] 0

t [epsilon] T, j [[epsilon] {[I.sup.t] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [J.sup.t]}, l [epsilon] [L.sub.j],

[[sigma].sub.l[epsilon][L.sub.j]] [y.sub.jl] [less than or equal to] 1 j [epsilon] {I [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] J [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] W [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] X},

[n.sub.jm] - Max[N.sub.jm] [[sigma].sub.l[epsilon][L.sub.j]] [y.sub.jl] [less than or equal to] 0,

[[sigma].sub.[s.sup.[epsilon]][S.sup.t]] [S.sub.tps] - [[sigma].sub.j[epsilon][I.sup.t]] [[sigma].sub.l[epsilon][L.sub.j] [V.sub.tpjl] = 0 t [epsilon] T, p [epsilon] P,

[[sigma].sub.j[epsilon][W.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [V.sub.tpjl] - [[sigma].sub.j[epsilon][J.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [V.sub.tpjl] = 0 t [epsilon] T, p [epsilon] P,

[[sigma].sub.j[epsilon][I.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [V.sub.tpjl] = [TD.sub.tp] t [epsilon] T, p [epsilon] P,

[[sigma].sub.t[epsilon]T][[[sigma].sub.p[epsilon]P] [[[sigma].sub.s[epsilon][S.sup.t]] [[alpha].sub.tps][S.sub.tps] + [[sigma].sub.j[epsilon][W.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [[gamma].sub.tpj][v.sub.tpjl]

+ [[sigma].sub.j[epsilon][X.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [[[epsilon].sup.d].sub.tpj][V.sub.tpjl] + [[sigma].sub.k[epsilon][K.sup.t]] [[z.sup.d].sub.tpk][D.sub.tpk]]]]

- [[sigma].sub.t[epsilon]T][[[sigma].sub.[p.sup.[epsilon]P]][[[sigma].s ub.j[epsilon][I.sup.t]] [[sigma].sub.l[epsilon][L.sub.j] [[[beta].sup.d].sub.tpj][V.sub.tpjl] + [[sigma].sub.j[epsilon][J.sup.t] [[sigma].sub.l[epsilon][L.sub.j] [[[delta].sup.d].sub.tpj][V.sub.tpjl]]] [less than or equal to] [Z.sub.0] d [epsilon] D,

where D is the complete set of dual extreme points of LP(v, z, y, n). ([IPDM.sup.*]) is obtained by relaxing most of the primal cuts, and is the master problem of the integrated production distribution in the primal decomposition scheme. More importantly, it is equivalent to the resource allocation problem with the additional variable [Z.sub.0] and the additional constraints of primal cuts generated by a small number of the extreme points of LP(v, z, y, n). The last three constraint sets before the cuts force ([IPDM.sup.*]) to generate solutions that guarantee the feasibility of LP(v, z, y, n) by ensuring there is enough capacity at every stage.

At each iteration, the primal cut restricts ([IPDM.sup.*]) from the solution of the previous iteration by imposing the transportation costs associated with the configuration, and forces a configuration that has a better total cost than the best overall solution attained so far. The iterations stops when ([IPDM.sup.*]) fails to find a better potential resource allocation configuration than the best feasible solution so far. The solution algorithm presented so far is a standard implementation of the primal or Benders decomposition algorithm. The algorithm has the advantage that it significantly reduces the memory requirements compared to a monolithic model and solution algorithm. However, significant reductions in solution times can be obtained by further exploiting the problem structure.

4.2. Acceleration techniques

The transportation subproblem exhibits significant primal degeneracy; i.e., there are many dual solutions associated with the primal optimal solution. The selected dual optimal solution has an effect on the strength of the generated cut and hence on the number of cuts needed to prove optimality. Let the cost on the transportation channel between the nodes i and j for product p in season t, denoted by [C.sub.tpij], be the sum of the cycle and pipeline inventory and the transportation costs.

The dual of the transportation problem is expressed by:

Max [[sigma].sub.t[epsilon]T][[[sigma].sub.[p.sup.[epsilon]P]] [[[sigma].sub.[s.sup.[epsilon]][S.sup.t]] [[alpha].sub.tps][S.sub.tps] + [[sigma].sub.j[epsilon][W.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [[gamma].sub.tpj][V.sub.tpjl]

+ [[sigma].sub.j[epsilon][X.sup.t]] [[sigma].sub.l[epsilon][L.sub.j] [[epsilon].sub.tpj][V.sub.tpjl] + [[sigma].sub.k[epsilon][K.sup.t]] [[xi].sub.tpk][D.sub.tpk]]]

- [[sigma].sub.[t[epsilon].sup.T]] [[[sigma].sub.[p.sup.[epsilon]P]] [[[sigma].sub.j[epsilon][I.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [[beta].sub.[tpjv.sub.tpjl]] + [[sigma].sub.j[epsilon][J.sup.t]] [[sigma].sub.l][epsilon][L.sub.j]] [[delta].sub.tpj][V.sub.tpjl]]],

s.t.

[[alpha].sub.tps] - [[beta].sub.tpj] [less than or equal to] [C.sub.tpsj] t [epsilon] T, p [epsilon] P, s [epsilon] [S.sup.t], j [epsilon] [I.sup.t],

[[gamma].sub.tpi] - [[delta].sub.tpj] [less than or equal to] [C.sub.tpij] t [epsilon] T, p [epsilon] P, i [epsilon] [W.sup.t], j [epsilon] [J.sup.t],

[[epsilon].sub.tpi] - [z.sub.tpk] [less than or equal to] [C.sub.tpik] t [epsilon] T, p [epsilon] P, i [epsilon] [X.sup.t], j [epsilon] [K.sup.t].

The objective function of the dual problem contains the solution of the resource allocation problem as the coefficients. Any change on the optimal dual variables with zero coefficients will not change the objective function value. Using its structure, the primal cut is strengthened by revising the non-basic dual variables within the dual feasibility constraints above, i.e., increasing [alpha], [gamma], [epsilon] and decreasing [beta], [delta].

4.3. Disaggregation of the primal cut

A disaggregation of the primal cut is obtained by exploring the subproblem structure. The subproblem consists of \P\ disconnected networks for each product, which can also be separated into \Tdisconnected problems for each season. As it is presented in the dual feasibility constraints, the transportation channels are grouped in three disconnected networks, i.e., from suppliers to first stage production facilities, from first stage warehousing facilities to second stage production facilities, and from second stage warehousing facilities to customers. Using this disconnected networks structure of the sub-problem, the primal cut is strengthened by being replaced by 3 x \P\ x \T\ constraints below.

[[sigma].sub.s[epsilon][S.sup.t]] [[[alpha].sup.d].sub.tps][S.sub.tps] - [[sigma].sub.j[epsilon][I.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [[[beta].sup.d].sub.tpj][V.sub.tpjl] [less than or equal to] [Z.sub.1tp] t [epsilon] T, p [epsilon] P,

[[sigma].sub.j[epsilon][W.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [[[gamma].sup.d].sub.tpj] [V.sub.tpjl] - [[sigma].sub.j[epsilon][J.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [[[delta].sup.d].sub.tpj] [V.sub.tpjl] [less than or equal to] [Z.sub.2tp] t [epsilon] T, p [epsilon] P,

[[sigma].sub.j[epsilon][X.sup.t]] [[sigma].sub.l[epsilon][L.sub.j]] [[[epsilon].sup.d].sub.tpj] [V.sub.tpjl] - [[sigma].sub.k[epsilon]K] [[z.sup.d].sub.tpk] [D.sub.tpk] [less than or equal to] [Z.sub.3tp] t [epsilon] T, p [epsilon] P.

To incorporate these constraints in ([IPDM.sup.*]), [Z.sub.0] variable in the objective function is replaced by

[[sigma].sub.t[epsilon]T] [[sigma].sub.[p.sup.[epsilon]P] ([Z.sub.1tp] + [Z.sub.2tp] + [Z.sub.3tp]).

4.4. Initial set of primal cuts

One of the major shortcomings of the primal decomposition approach is that the stability and the performance of the algorithm significantly depend on the initial values of dual variables used. The feasibility conditions of the dual variables are independent of the variables of the master problem. Hence relaxing the integer restrictions on the y and n variables of the master problem would yield feasible dual variables of the sub-problem to the unrelaxed problem. At the initial stages of the algorithm, these dual variables would provide the initial primal cuts.

4.5. Optimality checking mechanism

Solving the master problem to optimality at each iteration takes a significant amount of computing time. Especially at the initial stages of the algorithm, the primal cuts do not contain good information on the transportation costs. The alternative approach is to use the master problem as an optimality checking mechanism. A constraint on the objective function would force the master problem to search for a solution better than the best available total cost, which provides an upper bound to the optimal solution value. At each iteration the right hand side of the constraint is updated with the current upper bound, and continued until the master problem fails to find a better feasible solution. In this case, the master problem ceases to generate a lower bound to the original problem. The flow diagram of the overall solution procedure is shown in Fig. 5.

5. Experimental evaluation

5.1. Test case description

The seasonal production-distribution model and the associated primal decomposition method described above were tested on a real-life logistics reorganization project of a company that supplies cardboard packages to breweries and soft drink manufacturers. The company ships 12 types of paper products from paper mills, through a two-stage manufacturing process, to more than 200 customers around the nation. Both stage one and stage two manufacturing lines are constrained by the available production hours on the machines. Different manufacturing lines can be installed with different processing capabilities and costs for each product. The demand exhibits a strong seasonal pattern. The production resources required differ significantly from cardboard type to cardboard type and from manufacturing line to manufacturing line. In order to satisfy the peak period's demand, additional manufacturing capacity must be purchased, inventory needs to be built up during the off period, or the company must use sub-contractors.

The objective of the study was to design a strategic plan for the allocation of the existing and new manufacturing lines to facilities and to products by considering the location of customers and the corresponding production and transportation costs. Because the system throughput was constrained by the production capacities, the allocation of customers to warehouses had previously not received sufficient attention and inefficient transportation patterns existed. All transportation was executed in full truck loads or rail boxcar loads from the paper mills to the first manufacturing phase.

5.2. Computational experience

Initially, the system was configured using the data from a single season. The demands for this season were the peak demands. The season ran from May through July and all costs and parameters were adjusted based on the length of this planning period. The size parameters of the test case and of the resulting mixed integer programming formulation are given in Table 1.

The total running time of a mixed integer programming algorithm depends on the interaction between a variety of factors, such as the strength of the integer cuts and linear and integer programming tolerances. The disaggregation of the cuts in the master problem and the strengthening of the dual variables reduced the running times from 865 seconds for a standard implementation of the primal decomposition method to 21 seconds, a reduction by a factor of more than 40.

The multi-season model divided the year into three unequal seasons. The seasons run from January through April, May through July, and August through December, respectively. All costs are adjusted based on the length of the corresponding season. The number of customers was three times as large to represent the three seasons. For the seasonal model, the cuts are also disaggregated by season. This reduced the solution times further from 1958 seconds to 154 seconds, a reduction by a factor of more than 12. Detailed statistics on the data and the computational experiment can be found in Dogan [15]. Comparing the solution configurations of the one-season and seasonal model showed that inventory was built up during the off-peak periods, which in turn reduced the number of required machines from 19 to 16. The average cost of paper product was reduced by $20 per ton when inventory buildup was allowed, which yielded approximate savings of $8.3 million dollars per year. The total cost was approximately $401 million dol lars per year, for savings of 2% versus the best case without inventory. These savings clearly show the importance of integrating tactical production and inventory decisions, when making strategic decisions such as the location of manufacturing plants and manufacturing lines for logistics systems with seasonal demand variances. It should be noted that both configurations are optimal with respect to their objective functions and the savings differential is possible only because the multiperiod system has more flexibility. It is expected that the savings of the integrated approach will increase with the size of the variance in the seasonal demand.

The effect of the ratio of transportation cost to total system cost was also investigated for a reduced case with only the four most popular products. In the original data, the transportation cost accounted for about 15% of the total cost. Three more cases were generated by multiplying the transportation rates by two, five, and 10, respectively. The required running times increased from 19 seconds, to 26, 31, and 211 seconds, respectively, while the number of MIP iterations increased from two, to two, three, and five, respectively. This is consistent with the principle behind the primal decomposition, where zero MIP iterations are required if no transportation costs are present and a large number of MIP iterations is required if there are no production costs.

6. Conclusions

A comprehensive model for the integrated design of multi-period production-distribution systems has been developed for the case of seasonal customer demands. The model allows the simultaneous determination of vendor sourcing; facilities location, types, and sizes; production lines location and sizing; production allocation; inventory location and sizing; customer allocation; and transportation flows. It integrates the strategic decisions concerning facilities and production lines with the tactical decisions concerning production, inventory, and customer allocation. In the real life case study presented, this integrated approach saved the company an additional 2% over the hierarchical approach where optimal strategic and tactical decisions were made sequentially.

Such integrated and detailed models tend to be very large for industrial logistics systems. Standard branch-and-bound algorithms based on linear programming relaxations require unacceptably long running times and very large computer memories. A primal decomposition solution methodology has been developed to obtain optimal solutions for this model within acceptable limits of running times and memory requirements. The model and the problem structure allowed the additional implementation of specialized acceleration techniques. The computational experiment showed significant reduction in computation times caused by using dissagregated cuts over products, periods, and stages. In addition, the strengthening of the dual variables and using the LP relaxation in the initial iterations of the decomposition also significantly improved the efficiency of the solution procedure. The combined effect of all these acceleration techniques reduced the running time by a factor of 480 versus a standard Bender's decomposition alg orithm. This decomposition and acceleration methodology made it possible to solve industrial sized problems in a reasonable amount of time. Further runs made during the sensitivity analysis required even less time.

The developed decomposition and acceleration methodology can be applied to a large variety of practical integrated multi-period production-distribution problems and makes it feasible to obtain optimal configurations for these problems. This in turn yields significant cost reductions for manufacturing corporations.

A fertile area for future research would be the incorporation of linear product transformations in the production phases based on the bill of materials. The current model traces each product from supplier, through the various production stages, to the final customer. However in most industries, product transformations and assembly operations are common. The authors are currently incorporating bill of materials into the design of global distribution systems. The current model also focuses on a cyclic multi-period system where each facility and production line is either open or closed for all the periods. A second extension to the current model would allow the facilities to be opened and closed in different time periods. This would allow the design of logistics systems with growing or declining demand.

Biographies

Marc Goetschalckx is an Associate Professor in the School of Industrial and Systems Engineering and a faculty member of the Logistics Institute of the Georgia Institute of Technology. He received a M.E. in Mechanical Engineering from the Catholic University of Leuven, Belgium, an M.E in Industrial Engineering of the University of Florida and a Ph.D. Degree from the Georgia Institute of Technology. He is a Senior Member of the Institute of Industrial Engineers and a member of INFORMS. He is on the Editorial Board of IIE Transactions. He has functioned as Research Director and Director of the Facilities Planning and Design division of IIE. His research areas are in the design and operation of industrial logistics systems such as supply chains, global logistics systems, transportation systems, facilities design, material handling and network design. He has authored numerous papers on these topics in journals such as IIE Transactions, Management Science, European Journal of Operational Research, and Journal of Bu siness Logistics and is a frequent speaker on these topics in continuing education short courses.

Koray Dogan is currently a Senior Solutions Consultant at Advanced Solutions Center of i2 Technologies Logistics Division. He received B.S. in Industrial Engineering from Bilkent University of Turkey, in 1991, M.S. and Ph.D. degrees in Industrial Engineering from Georgia Institute of Technology in 1992 and 1996, respectively. His development and consulting experience lies in various aspects of supply chain management; network design, demand management, inventory management, production and distribution planning and its integration to strategic network design problems.

(1.) i2 Technologies, Logistics Division, 50 Main Street, Suite 1000, White Plains, NY 10606, USA E-mail: KorayDogan@i2.com

(2.) School of Industrial and Systems Engineering, 765 Ferst Drive, Atlanta, GA 30332-0205, USA E-mail: marc.goetschalckx@isye.gatech.edu

References

(1.) Geoffrion, A.M. and Powers, R.F. (1995) Twenty years of strategic distribution system design: an evolutionary perspective, (implementation in OR/MS: an evolutionary view). Interfaces, 25, 105-128.

(2.) Thomas, D. and Griffin, P.M. (1996) Coordinated supply chain management. European Journal of Operational Research, 94, 1-15.

(3.) Vidal, C. and Goetschalckx, M. (1997) Strategic production-distribution models: a critical review with emphasis on global supply chain models. European Journal of Operational Research, 98, 1-18.

(4.) Aikens, C.H. (1985) Facility location models for distribution planning. European Journal of Operational Research, 22, 263-279.

(5.) van Roy, T. (1983) Cross decomposition for mixed integer programming. Mathematical Programming, 25, 46-63.

(6.) van Roy, T. (1986) A cross decomposition algorithm for capacitated facility location. Operations Research, 34, 145-163.

(7.) Geoffrion, A.M. and Graves, G.W. (1974) Multicommodity distribution system design by Benders decomposition. Management Science, 20, 822-844.

(8.) Cohen, M.A. and Lee, H.L. (1988) Strategic analysis of integrated production-distribution systems: models and methods. Operations Research, 36, 216-228.

(9.) Cohen, M.A. and Lee, H.L. (1989) Resource deployment analysis of global manufacturing and distribution networks, Journal of Manufacturing Operations Management, 2, 81-104.

(10.) Arntzen, B.C., Brown, G.G., Harrison, T. P. and Trafton, L.L. (1995) Global supply chain management at Digital Equipment Corporation. Interfaces, 25, 69-93.

(11.) Lee, C. (1991) An optimal algorithm for the multiproduct capacitated facility location problem with a choice of facility type. Computers and Operational Research, 18, 167-182.

(12.) Lee, C. (1993) A cross decomposition algorithm for a multiproduct-multitype facility location problem. Computers and Operational Research, 20, 527-540.

(13.) Goetschalckx, M., Nemhauser, G., Cole, M.H., Wei, R., Dogan, K. and Zang, X. (1994) Computer aided design of industrial logistic systems, in Proceedings of the Third Triennial Symposium on Transportation Analysis (TRISTAN III), pp. 151-178.

(14.) Dogan, K. and Goetschalckx, M. (1997) A primal decomposition method for the integrated design of multi-period production distribution systems. Research Report, School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA.

(15.) Dogan, K. (1996) A primal decomposition scheme for the design of strategic production distribution systems. Unpublished Ph.D. dissertation, Georgia Institute of Technology, Atlanta, GA.

                 Test case and formulation size parameters
Item                              1 season 3 seasons
Number of products                    12        12
Number of locations                    6         6
Number of manufacturing lines         81        87
Number of customers                  238       714
Number of transportation channels   1358      3738
Number of nodes                      355      1065
Number of arcs                      1470      4074
Number of capacitated arcs            88       102
Number of integer variables           88       102
Number of continuous variables     17640     48888