# Design of a Distribution Network Using Primal-Dual Decomposition

Authors: J. A. Marmolejo; R. Rodríguez; O. Cruz-Mejia; J. Saucedo
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

A method to solve the design of a distribution network for bottled drinks company is introduced. The distribution network proposed includes three stages: manufacturing centers, consolidation centers using cross-docking, and distribution centers. The problem is formulated using a mixed-integer programming model in the deterministic and single period contexts. Because the problem considers several elements in each stage, a direct solution is very complicated. For medium-to-large instances the problem falls into large scale. Based on that, a primal-dual decomposition known as cross decomposition is proposed in this paper. This approach allows exploring simultaneously the primal and dual subproblems of the original problem. A comparison of the direct solution with a mixed-integer lineal programming solver versus the cross decomposition is shown for several randomly generated instances. Results show the good performance of the method proposed.


---

## Body

## 1. Introduction

The constant emphasis on customer satisfaction has highlighted the importance of designing distribution networks of firms [1]. Optimal network design plays an important role in the supply chain operation, as good logistics distribution network can save transportation costs as well as improve customer service levels [2]. Facility location can be a criterion for the design of distribution networks. Many organizations consider facility location as a strategic decision for having high material handling cost.Distribution costs in many industries constitute an important part of the total logistics expenditure. Consequently, the final price is strongly linked to the location of facilities where products are manufactured or stored; see Figure1. In particular, a cost to take into consideration is the fix cost of opening and setting up new facilities as temporal consolidation centers or mixing centers, that is, cross-docking or merge-in-transit centers [3].Figure 1
Three-echelon supply chain.Given that the complexity of the mathematical model falls in Np-hard, we propose an efficient method to solve the design of a soda bottling distribution system, depicted in Figure2. The proposed distribution network is constituted by plants, warehouses (cross-docking and merge-in-transit centers), and distribution centers. Commodities are produced at several capacitated plants and the demand of distribution centers is satisfied from warehouses.Figure 2
Distribution system of the company.The problem is to determine among the possible warehouses the ones to be opened to consolidate the demand and by which warehouse each distribution center is exclusively served. The objective is minimizing fix costs and total transportation costs and to establish a network of routes that enables the flow of products in order to satisfy some demand characteristics. The proposed network is shown in Figure3.Figure 3
Proposed distribution system.


## 2. Literature Review

Network design problems with real scenarios are widely addressed in optimization problems. Most of them study four supply chain functions: location, production, inventory, and transportation with the aim of integrating them. In [4], an integrated optimization model of supply chain functions in a multiplant, multiproduct, and multicustomer supply chain with deterministic demand is developed and formulated as a capacitated location, production, and distribution problem.More recent papers include the study of real cases focusing on particular aspects of the location problem. In [5], the authors present a review in the context of supply chain management and its integration with other decisions in the context of network design. For example, in [6, 7], interesting aspects in the distribution network design are enumerated including the classical facility location problem solved by different techniques according to the specific objectives.Nowadays the competitive and ever-changing business environment makes the distribution network design more complicated. New features including specific conditions with suppliers, distribution centers, and customers are modeled. A seminal work was done by [8] in which an extended view of the distribution network design including suppliers, facilities (production and distribution), customers, and many kinds of transportation means is included.Research work by [9] provides a literature review about the main research papers published from 2005 to 2015 under the following keywords: distribution networks design, supply chain, logistics, and global. The review includes a description of the operations research techniques used to solve each problem.Numerous methods of solution have been used to solve facility location problems including [10], in which a new method for the solution of the problem addresses the optimal location of distribution centers between plants and customers. Because of the dimension of the problem, they develop an algorithm based on Benders decomposition (BD) [11] for solving a multicommodity distribution network. Reference [12] introduces a mixed-integer problem (MIP) to model a multiproduct distribution network solved with BD and two Tabu Search heuristics that made possible the convergence and solution quality. Another classical model is found in [13], in which a model solves a minimization function that includes fix costs in warehouses, distribution centers, and transportation costs for multicommodities from plants to warehouses and finally to customers. Similarly in [13], a triechelon, multicommodity system including production, distribution, and transportation planning is solved using Lagrangian Relaxation (LR) in a heuristic. In [14–16], a facility location problem is also solved but they consider the specific features to solve the problem.In [17] a distribution network design problem with 3 stages is solved optimizing the numbers of locations and capacities of plants and warehouses. At the same time the problem minimizes total costs and satisfies all demands. Given the complexity of the solution the authors use a Lagrangian based solution procedure for the problem. In [18] a variable neighborhood search (VNS) heuristic method integrating a Tabu Search procedure is used to solve a large scale production and distribution network design model.A different solution strategy is presented by [19] based on bilevel programming problem. In this paper a distribution network with distribution and production plants around the world is solved using an extended genetic algorithm. Another multicriteria proposal is given by [20] whose authors solve the problem with MIP (is not large scale problem) making tactical decisions for distributing the product to customers. In [21], a multiobjective genetic algorithm is developed to solve a stochastic production-distribution network in which the objective function is optimizing the costs and service level.In [22], an integrated distribution network design and site selection problem is analyzed. The setting is in the context of transportation planning faced by the freight-forwarding industry. The problem includes a strategic level multicommodity network design. In this problem each commodity is defined by a unique pair of origin and destination points and a required amount of flow and other considerations proper to the real world using BD.In [23, 24], also a BD approach is used, first in combination with an intelligent algorithm to improve the time solution for the main problem and later in a modified version that takes advantage of the mathematical formulation. In both cases the problem is in deterministic, multicommodity, and single period contexts.BD or primal decomposition methods exploit only the primal structure of the original problem. However, many mixed-integer programming problems have primal and dual easy-to-solve problems; for example, in [25–27], cross decomposition (CD) is used to solve a capacitated facility location problem by defining subproblems of transportation and location of plants.In conclusion, the reviewed research papers addressing the facility location and transportation planning problems have missed modeling cross-docking and merge-in-transit centers which are a novel feature in this paper.In this paper, we propose single-sourcing constraints ensuring that each distribution center is exclusively served by a single cross-docking or merge-in-transit center; see Figure2. This condition arises often in bottled drinks companies where the operating conditions restrict that only one cross-docking center serves one distribution center. Additional operating conditions from the real case are included in the constraints and explained in Section 3.The model considers binary variables for cross-docking locations and the allocation of cross-docks to distribution centers. Continuous variables for the flow of a single commodity from plants to operating cross-docks are defined. This is a common problem for bottling companies to define specific supply routes to each distribution center. Using the proposed mathematical model more efficient transportation routes will be generated.We assume that the design of the distribution network can be solved efficiently by using decomposition techniques, more specifically a primal-dual decomposition. This method was originally developed for linear mixed-integer programming problems but the approach is more general and not restricted to such problems.Many combinatorial optimization problems can be solved if the complexity of variables and constraints were removed. Some examples are the assignment problem, the facility location problem [28], the optimal power flow [29], and other mixed-integer programming problems.Primal-dual decomposition can offer also better computational time than traditional decomposition techniques, such as Benders decomposition and Lagrangian Relaxation. However, the key point for having good results using primal-dual decomposition techniques is the mathematical structure of the problem. By this we mean if the rows and columns of the coefficient matrix can be rearranged so that the matrix has block-angular form, then primal-dual decomposition method will generate better solutions.The paper is organized as follows. Section3 presents a mathematical programming model of the distribution problem. Section 4 describes the cross decomposition that is the solution methodology. In Sections 5 and 6 we present the computer implementation and experimental results. Conclusions are reported in Section 7.


## 3. Mathematical Model

LetK be the set of manufacturing plants. An element k∈K identifies a specific plant of the company. Let I be the set of the potential cross-docking warehouses. An element i∈I is a specific cross-docking warehouse. Finally, let J be the set of distribution centers; a specific distribution center is any j∈J. Let Z denote the set of integers {0,1}.Parameters. Consider the following:Qk = capacity of plant k.βi = capacity of cross-docking warehouse i.Fi = fixed cost of opening cross-docking warehouse in location i.Gki = transportation cost per unit of the product from plant k to the cross-docking warehouse i.Cij = cost of shipping the product from cross-dock i to the distribution center (CeDis) j.dj = demand for the distribution center j.Decision Variables. We have the following sets of binary variables to make the decisions about the opening of the cross-docking warehouse and the distribution for the cross-docking warehouse to the distribution center:(1)Yi=1If  locationiis  used  as  a  cross-docking  warehouse,0otherwise,Xij=1if  cross-dockisupplies  the  demand  of  CeDisj,0otherwise.Wki is the amount of product sent from plant k to the cross-dock i which is represented by continuous variables.We can now state the mathematical model as a (P) problem. See [30]. Consider(2)minWki,Yi,XijZ=∑k∈K∑i∈IGkiWki+∑i∈IFiYi+∑i∈I∑j∈JCijdjXij.Subject to Constraints
Capacity of the plant is as follows:(3)∑i∈IWki≤Qk,∀k∈K.
Balance of product is as follows:(4)∑j∈JdjXij=∑k∈KWki,∀i∈I.
Single cross-docking warehouse to distribution center is as follows:(5)∑i∈IXij=1,∀j∈J.
Cross-docking warehouse capacity is as follows:(6)∑j∈JdjXij≤βiYi,∀i∈I.
Demand of items is as follows:(7)pYi≤∑k∈KWki,∀i∈I(8)p=min⁡dj(9)Wki≥0,∀i∈I,∀k∈K(10)Yi∈Z,∀i∈I(11)Xij∈Z,∀i∈I,∀j∈J.The objective function (2) considers in the first term the cost of shipping the product from plant k to the cross-docking warehouse i. The second term contains the fix cost to open and operate the cross-docking warehouse i. The last term incorporates the cost of fulfilling the demand of the distribution center j. Constraint (3) implies that the output of plant k does not violate the capacity of plant k. Balance constraint (4) ensures that the amount of products that arrive to a distribution center j is the same as the products sent from plant k. The demand of each distribution center j will be satisfied by a single cross-docking warehouse i, and this is achieved by constraint (5). Constraint (6) bounds the amount of products that can be sent to a distribution center j from an opened cross-docking warehouse i. Constraint (7) guarantees that any opened cross-docking warehouse i receives at least the minimum amount of demand requested by a given distribution center j. Constraint (8) ensures that the minimum demand of each distribution center j is considered. Finally, constraints (9), (10), and (11) are the nonnegative and integrality conditions.


## 4. Cross Decomposition Approach

Many of the large scale mixed-integer linear programming problems are too complex to be solved directly with commercial software. However, when the computational complexity grows exponentially according to the instance size, decomposition techniques usually offer better solutions. As it is seen before, in this cases Benders decomposition is used as well as the Lagrangian Relaxation. In general terms, Benders decomposition generates good results. However, the master problem of Benders can be difficult to solve and require very large computation time [31, 32].In this paper, with very large instances including several plants, cross-docks, and distribution centers, the optimal design of the distribution network involves a very large number of integer binary variables that generate large computational time for reaching the optimal solution.Because Benders method sometimes produces better bounds than Lagrange method but the solution of its master problem involves large amount of computational time, in this paper we use cross decomposition to obtain automatically the best primal and dual bounds and an approximate optimal solution. Cross decomposition can obtain exact and finite solutions for mixed-integer models if the continuous part of the problem is convex and linear. This decomposition method unifies Benders decomposition and Lagrangian Relaxation into a single framework that involves iterative solutions to a primal subproblem (Benders) and a dual subproblem (Lagrange).Cross decomposition [28] consists basically in a first stage of solving subproblems: primal and dual. The dual subproblem generates the binary variables (Y¯i and X¯ij) that will be fixed for the primal subproblem. Additionally, the dual subproblem provides a lower bound of the optimal value for the objective function (ΦDS). The primal subproblem generates the Lagrangian multipliers (λ¯i) that will be fixed for the dual subproblem. At the same time, the primal subproblem produces an upper bound of the optimal value for the objective function (ΩPS). In each iteration a primal and a dual convergence test will be performed (CTP and CTD). If any of these tests fail, there will be a need to solve the master problem; see Figure 4:CTP = IfΩPS<v¯ go to solve the primal subproblem, and if not, go to solve Dual Master Problem.CTD = IfΦDS<v_ go to solve the dual subproblem, and if not, go to solve Dual Master Problem.Figure 4
Cross decomposition algorithm.v
¯ is the least upper bound known and v_ is the largest lower bound known at the current step of the algorithm. Consequently, this method can be used to reduce cpu time of the original problem [29]. For large scale instances, that is, problems with large number of plants, cross-docks, and distribution centers, the direct solution is very complicated and cpu time increases exponentially. The exponential growth of the cpu time is related to selecting what cross-dock must supply the demand for the distribution center.In this work we use only the Dual Master Problem because the subgradient method is implemented; see [33, 34]. A good but not necessary optimal set of multipliers is obtained by subgradient optimization. The starting multipliers can be set equal to zero. The objective is to limit as much as possible the use of any master problem. The master problem is more difficult to solve than the primal or dual subproblems. To generate the dual subproblem we relax the cross-dock balance constraint (4). λ¯i is the Lagrangian multiplier of dualized constraint (4). By dualizing this constraint we obtain a dual subproblem that is less expensive to solve. This relaxation also speeds up the solution of this subproblem. Primal and dual convergence tests are used to check a bound improvement and to verify when an optimal solution is reached.The procedure generates a primal bound (ΩPS) and a dual bound (ΦDS) with corresponding optimal values. A disadvantage is that, for nonconvex problem, the convergence cannot be guaranteed because there is no monotonic improvement of the bounds. For this reason, the procedure includes a convergence test that makes sure of obtaining a better bound. Next we describe the primal and dual subproblems obtained from the original problem.Dual Subproblem (DS). Consider the following:(12)minWki,Yi,XijΦDS=∑k∈K∑i∈IGkiWki+∑i∈IFiYi+∑i∈I∑j∈JCijdjXij+∑i∈Iλ¯i∑k∈KWki-∑j∈JdjXijsubject to constraints (3), (5)–(11).Primal Subproblem (PS). Consider the following:(13)minWkiΩPS=∑k∈K∑i∈IGkiWki+∑i∈IFiY¯i+∑i∈I∑j∈JCijdjX¯ijsubject  to:∑i∈IWki≤Qk,∀k∈K∑j∈JdjX¯ij=∑k∈KWki,∀i∈I∑j∈JdjX¯ij≤βiY¯i,∀i∈IpY¯i≤∑k∈KWki,∀i∈Ip=min⁡djWki≥0,∀i∈I,∀k∈K.Dual Master Problem (DMP). Consider the following:(14)maxλiΨDMP=ΦDSsubject to:λi≥0.


### 4.1. Cross Decomposition Algorithm

Step 1 (initialize).
Select initial valuesλi for the Lagrangian multipliers and set up and apply the corresponding CTD. The starting multipliers can be set equal to zero. Either stop (the algorithm terminates when v¯-v_≤ϵ) or go to Step 4 or set up the dual subproblem.Step 2.
Solve the dual subproblem (DS) that is a lower bound. Apply CTP forY¯i and X¯ij. Either stop (the algorithm terminates when v¯-v_≤ϵ) or go to Step 4 or set up the primal subproblem.Step 3.
Solve the primal subproblem (PS) that is an upper bound. Apply CTD forλ¯i. Either stop (the algorithm terminates when v¯-v_≤ϵ) or go to Step 4 or set up the dual subproblem corresponding to the optimal solution of the current (PS) and go to Step 2.Step 4.
Solve the master problem. In this work we solve the Dual Master Problem (DMP). Find new values for the Lagrangian multipliers that are held fix in DS. Set up the corresponding subproblem and go, respectively, to Step2. In this case we use subgradient method to solve DMP.


## 4.1. Cross Decomposition Algorithm

Step 1 (initialize).
Select initial valuesλi for the Lagrangian multipliers and set up and apply the corresponding CTD. The starting multipliers can be set equal to zero. Either stop (the algorithm terminates when v¯-v_≤ϵ) or go to Step 4 or set up the dual subproblem.Step 2.
Solve the dual subproblem (DS) that is a lower bound. Apply CTP forY¯i and X¯ij. Either stop (the algorithm terminates when v¯-v_≤ϵ) or go to Step 4 or set up the primal subproblem.Step 3.
Solve the primal subproblem (PS) that is an upper bound. Apply CTD forλ¯i. Either stop (the algorithm terminates when v¯-v_≤ϵ) or go to Step 4 or set up the dual subproblem corresponding to the optimal solution of the current (PS) and go to Step 2.Step 4.
Solve the master problem. In this work we solve the Dual Master Problem (DMP). Find new values for the Lagrangian multipliers that are held fix in DS. Set up the corresponding subproblem and go, respectively, to Step2. In this case we use subgradient method to solve DMP.


## 5. Computer Implementation

In this section we compare the computational implementation for the direct solution versus the solution obtained using cross decomposition. Both were solved using a commercial software GAMS [35]. We randomly generated 30 instances according to the structure and complexity of the real case instance. The size of an instance is given by the number of manufacturing plants (|K|), the number of cross-docking warehouses (|I|), and the number of distribution centers (|J|). |J| are chosen randomly in set {254,…,1000} according to a uniform distribution. |K| and |I| follow the same proportion of the real case instance.In Table1, the instances generated with different numbers of plants, cross-docks, and distribution centers are observed. Additionally, an instance based on a real case with |K|=44, |I|=56, and |J|=254 was solved. We test the solution method under different circumstances to evaluate the performance under different complexity instances. The full-scale model and the decomposition strategy proposed were implemented in GAMS using the solver CPLEX [36] for the mixed-integer programming (DMP, DS) and the linear problems (PS). All mathematical models were carried out on AMD Phenom II N970 Quad-Core with a 2.2 GHz processor and 4 GB RAM. We set GAMS parameter OPTCR at 0.0015; that is, the relative termination tolerance is within 0.15% of the best possible solution. Additionally, the size of all MIP models was reduced through presolver phase of CPLEX. The cross decomposition algorithm stops when the values of the lower and upper bounds are equal, except for a small tolerance ϵ = 0.15%:(15)ϵ=UB-LBUB·100%.Table 1
Size and characteristics of problem instances.
Instances
|
K
|
|
I
|
|
J
|
Continuous variables
Binary variables
Constraints
INST-1
78
99
450
7,722
44,649
825
INST-2
58
73
332
4,234
24,309
609
INST-3
117
148
673
17,316
99,752
1,234
INST-4
110
140
635
15,400
89,040
1,165
INST-5
115
147
665
16,905
97,902
1,221
INST-6
173
220
998
38,060
219,780
1,831
INST-7
51
65
293
3,315
19,110
539
INST-8
106
135
614
14,310
83,025
1,125
INST-9
128
163
738
20,864
120,457
1,355
INST-10
124
157
713
19,468
112,098
1,308
INST-11
88
112
510
9,856
57,232
934
INST-12
66
84
382
5,544
32,172
700
INST-13
129
165
747
21,285
123,420
1,371
INST-14
81
103
466
8,343
48,101
856
INST-15
172
219
992
37,668
217,467
1,821
INST-16
80
102
463
8,160
47,328
849
INST-17
146
186
842
27,156
156,798
1,546
INST-18
113
144
654
16,272
94,320
1,199
INST-19
105
134
606
14,070
81,338
1,113
INST-20
145
185
837
26,825
155,030
1,537
INST-21
98
125
568
12,250
71,125
1,041
INST-22
61
78
354
4,758
27,690
649
INST-23
66
84
381
5,544
32,088
699
INST-24
122
156
706
19,032
110,292
1,296
INST-25
172
219
993
37,668
217,686
1,822
INST-26
68
86
391
5,848
33,712
717
INST-27
141
179
813
25,239
145,706
1,491
INST-28
114
145
659
16,530
95,700
1,208
INST-29
53
67
305
3,551
20,502
559
INST-30
140
178
807
24,920
143,824
1,481


## 6. Experimental Results

The design of the distribution network studied in this paper was undertaken using an algorithm of cross decomposition described in earlier sections. Table2 illustrates the cpu times of 30 instances by proposed decomposition strategy. It is also shown that the approximate optimal solution is very close to the optimal/best found integer feasible solution. The maximum gap was 0.35%:(16)gap⁡=GAMS  Solution-CD  SolutionGAMS  Solution·100%.Table 2
Computational statistics.
Instances
|
K
|
|
I
|
|
J
|
GAP (%)
cpu time (seconds)
INST-1
78
99
450
0.20
954
INST-2
58
73
332
0.20
549
INST-3
117
148
673
0.29
1734
INST-4
110
140
635
0.33
1259
INST-5
115
147
665
0.36
1632
INST-6
173
220
998
0.25
3500
INST-7
51
65
293
0.25
495
INST-8
106
135
614
0.24
1476
INST-9
128
163
738
0.26
1847
INST-10
124
157
713
0.33
1810
INST-11
88
112
510
0.35
921
INST-12
66
84
382
0.28
769
INST-13
129
165
747
0.34
1930
INST-14
81
103
466
0.26
1093
INST-15
172
219
992
0.30
3180
INST-16
80
102
463
0.22
852
INST-17
146
186
842
0.28
2031
INST-18
113
144
654
0.29
1027
INST-19
105
134
606
0.32
921
INST-20
145
185
837
0.32
2160
INST-21
98
125
568
0.25
1328
INST-22
61
78
354
0.34
465
INST-23
66
84
381
0.28
643
INST-24
122
156
706
0.24
1759
INST-25
172
219
993
0.30
3420
INST-26
68
86
391
0.32
865
INST-27
141
179
813
0.30
2090
INST-28
114
145
659
0.32
1426
INST-29
53
67
305
0.27
539
INST-30
140
178
807
0.31
2647In Figures5, 6, 7, 8, 9, and 10 are shown the results of six instances randomly generated. These results were obtained by the cross decomposition of the original problem.Figure 5
Comparison between the objective functions obtained by GAMS and by the proposed cross decomposition.Figure 6
Comparison between the objective functions obtained by GAMS and by the proposed cross decomposition.Figure 7
Comparison between the objective functions obtained by GAMS and by the proposed cross decomposition.Figure 8
Comparison between the objective functions obtained by GAMS and by the proposed cross decomposition.Figure 9
Comparison between the objective functions obtained by GAMS and by the proposed cross decomposition.Figure 10
Comparison between the objective functions obtained by GAMS and by the proposed cross decomposition.Looking at the results in Figures5, 6, and 8 it can be observed that the number of iterations and the performance of the lower and upper bounds pick up well the complexity of small instances.At the same time, Figures7, 9, and 10 show an increment in iterations required and the gap of the solution. These instances were selected to test the performance of the cross decomposition in very large and complex scenarios.Tables1 and 2 present the structure of the generated instances and, aditionally, the computational statistics of each one.


## 7. Conclusions

In this paper a primal-dual method is used to design a distribution network for bottled drinks company. Cross decomposition is a good method for solving large scale mixed-integer programming problems, especially when the resulting primal and dual subproblems are easy-to-solve, as in this case. This work proposes a decomposition scheme that reduced the computational time while maintaining the convergence of primal and dual solutions.If the duality gap of the Lagrangian Relaxation is small, the algorithm converges quickly to optimal or near-optimal solutions. Otherwise, the method needs other algorithms in order to obtain an exact solution. However, this method can be used as a heuristic which produces a feasible solution.Computational tests are presented using 30 random instances and real case data. The results show that the proposed solution strategy obtains a maximum gap of 0.35%. For these kinds of problems, we can often obtain an acceptable gap between approximate optimal solution (cross decomposition solution) and the optimal solution (CPLEX solution). For problems with a large duality gap, it is recommended to use a branch and bound algorithm to reduce this gap. For the problem studied in this paper there was no need to use heuristic to eliminate the duality gap. Because of this the use of convergence tests is recommended.Future research can be directed towards developing new procedures to obtain specific Lagrangian multipliers and improve the quality of the lower and upper bound. The solution of large scale mathematical problems using traditional methods takes large computational times. For this reason, the use of cross decomposition techniques allows the solution in shorter computational time. In this paper the cpu time was <3500 seconds. Additionally, cross decomposition methods can be used in parallel computing which decreases cpu time.The model proposed in this paper assumes deterministic parameters and does not consider a decomposition by stages. Solving the model by stages potentially can generate even faster solutions but at the same time produce a Lagrangian multiplier per subproblem. Large number of Lagrangian multipliers imply a method more sensitive to numerical stability.Future research for distribution networks using decomposition techniques can focus on the implementation on a grid computing platform that takes advantage of supercomputer nodes. This approach can offer a better usage of computational resources. The design of the distribution network can be advanced using innovative concepts of collaboration in supply chains as it is the vendor managed inventory, which is a way to integrate production and supply decisions reducing the delay of information.



---

*Source: 7851625-2016-02-08.xml*