Solving Capacitated Facility Location Problem Using Lagrangian Decomposition and Volume Algorithm

In this research, we will focus on one variant of the problem: the capacitated facility location problem (CFLP). In many formulations of the CFLP, it is assumed that each demand point can be supplied by only one open facility, which is the simplest case of the problem.We consider the case where each demand point can be supplied bymore than one open facility.We first investigate a Lagrangian relaxation approach. *en, we illustrate in the problem decomposition how to introduce tighter constraints, which solve the CFLP faster while achieving a better quality solution as well. At the same time, we apply the volume algorithm to improve both the lower and the upper bound on the optimum solution of the original problem for the large problem size.


Introduction
ere are a limited number of test case problems for the CFLP that we know of in the literature, which can be used to evaluate the results of our main goal of this research, the Lagrangian decomposition algorithm. Most of the algorithms developed in the literature have been tested on simulated data.
Lagrangian relaxation or Lagrangian decomposition was introduced in the early 70's through the work of Held and Karp [1,2] and Held et al. [3]. ey realized that the relationship between the systematic traveling-salesman problem and the minimum spanning tree problem gives a sharp lower bound on the optimum solution.
ereafter, they used a procedure called subgradient optimization [3] and showed that it is effective for approximating the maximum of certain piecewise linear concave functions. e concept of Lagrangian relaxation is that we relax some constraints and then penalize their violation, and the motivation for the relaxation of these constraints is that many combinatorial optimization problems consist of an easy problem with some additional complicating constraints. Applying the Lagrangian relaxation to these problems will relax these complicating constraints and absorb them into the objective function. e resulting relaxed problem becomes much easier to solve. Barahona and Anbil [4] presented an extension to the subgradient algorithm that will produce an approximation cost per iteration; they called it the volume algorithm. In general, it produces a primal vector as well as a dual vector that can be used as a starting point for a more exact method. ey were able to present a successful experiment with linear programming problems.
Alenezy and Khalaf in [12] discussed the Lagrangian decomposition and volume algorithm procedures in detail. In this article, we use variants of it in solving our CFLP for large size problems. e results obtained in this work were compared with the results of best greedy and greedy weak solutions.
In order to evaluate the performance of our Lagrangian decomposition algorithm, we develop a number of greedy heuristics (filters) that provide good solutions for known problems in the literature [9].
is is to illustrate the effectiveness of these heuristics for comparative purposes when evaluating our Lagrangian decomposition approach on large problem instances.
In the extreme case of very large problems, it is not possible to have any strong constraints. However, we take a dynamic approach to the model representation, whereby an effective tightening constraint is dynamically added to the model representation.
In particular, in this paper, we illustrate all the results and the analysis of our experiments for both the greedy heuristics and the Lagrangian decomposition algorithm. We also include tests on the performance of our algorithm compared to a number of benchmark models.

Capacitated Facility Location Problem
e facility location problem (FLP) seeks to locate a number of facilities to serve a number of customers; thus, there is a set of potential facility locations F; opening a facility at location i ∈ F has an associated nonnegative fixed cost f i and has either a limited or unlimited capacity S i of available supply. ere is a set of customers (clients) or demand points D that require service; customer j ∈ D has a demand d j that must be satisfied by the open facilities. If a facility at location i ∈ F is used to satisfy part of the demand of client j ∈ D, then there is a service or transportation cost incurred, which is often proportional to the distance from i to j, denoted by c ij .
Let F � a set of potential facility locations, D � a set of customers or demand points, m � the number of potential locations of the facilities; m � |F|, n � the number of customer; n � |D|, j d � the demand of customer j ∈ D(where d j ≥ 0), c ij � the unit cost of supplying the demand of customer j ∈ D from facility i ∈ F (where c ij ≥ 0), S i � the capacity of facility i (i.e., the upper limit on the total demand that can be supplied from facility i, where S i ≥ 0), p � the desired number of open facilities, f i � the fixed cost associated with opening facility i (where f i ≥ 0), x ij � the fraction of the demand of customer j supplied from facility i (where 0 ≤ x ij ≤ 1), and y i � 1 if facility i is open 0 otherwise . en, the formulation of the CFLP as a mixed-integer programming problem, which is referred to as IP found in [12], is as follows. e objective equation: (1) is equation ensures that the demand of each customer is satisfied: To ensure that the closed facility does not supply any customer and that the demand supplied from the facility does not exceed the capacity of the facility, is is the integrality constraint: is constraint tightens the lower bound representation of the CFLP. e value of p is determined as follows; we sort all of the facilities in descending order of capacity. en, p is such that, in the sorted order of facilities, Although the right-hand side of the constraint (5) need not be bounded, as a consequence of our experiments using p + 2 speeds up considerably the computational time with no detriment to the solution quality, e last inequality provides bounds on the allocation variables x ij .
In order to develop a Lagrangian heuristic for the CFLP, first we need to consider a linear programming relaxation for the IP problem, which is the same formulation (IP), except we replace inequalities (4) by We will denote the LP relaxation by P.

Lagrangian Decomposition and Volume Algorithm (LD and VA)
In this part, we first consider a Lagrangian relaxation of the above problem (P). en, we describe how to use the volume algorithm [7], which is an extension to the subgradient optimization.
In order to investigate the solution of large-sized problems, we followed the methods of Alenezy and Khalaf [12], at is, by decomposing CFLP into independent problems, which are easier to solve. Motivated by this method, let u i be the dual multiplier for the equation j in (1) and let c ij � d j c ij − u j . en, a lower bound L(u) is given by solving the following problem: subject to x ij ≤ y i , ∀i ∈ F, j ∈ D, (11)

Advances in Operations Research
3.1. e Lagrangian Decomposition (LD). It was widely reported that solving the L(u) above provides a good lower bound on the integer optimum solution. is is improved by using the volume algorithm. Motivated by this and by the fact that it is difficult to solve the L(u) for large size problems, we decompose the L(u) problem into m independent subproblems for each i ∈ F to compute a lower bound, we also at this point relax constraint (13), and we will return to these later in this section. e generic forms of each subproblem are Solving this to get a lower bound (LB) is easy, and it is easy too for the upper bound (UB). First, we set to 0 any variable x j with c j > 0. en, we assume that the rest of the variables are ordered such that , where n ′ ⊆ n and n � |D|. (16) Now, let k be the largest index such that k j�1 d j ≤ S, and let (17) Having solved these m independent subproblems, we next consider the minimum number of facilities that are needed to supply all the demand. In this way, we enhance the LB here by comparing the number of opened facilities denoted by h, after having solved the m subproblems, to the minimum needed p. If h < p, then we sort the unopened facilities by their fixed costs and open the cheapest facilities until we have p opened. e LB is suitably increased to account for these extra fixed costs.

e Volume Algorithm (VA).
To improve the lower bound we obtained from solving the decomposition of the problem above, we use the VA developed in [7] and the same algorithm found in [12], which stopped by a number of passes � 200, and next, we extend it to 2000. is algorithm can be formulated by the following steps: Step 1. Start with a vector u and solve (8)- (13) to obtain (x, y) and L(u) and set t � 1.
where α is a number between 0 and 1. In order to set the value of α, we solve the following one-dimensional problem: e value b was originally set to 0.1, and then every 100 iterations we check if L(u t ) had not increased by at least 1%, in which case, we divide b by 2; otherwise, we keep it as it is. When b becomes less than 10 − 5 , we keep it constant at this value. Step Step 4. Stopping criteria.
(1) v t j < 0.02 (2) (UB − L(u t ))/UB ≤ 0.02 (3) e number of passes equals 2000 e algorithm terminates when one of the previous criteria is met. If stopping rules are not satisfied, then set t � t + 1, and go to step 2. e formula of the step size s is as the one used in the subgradient method [13]: where λ is a number between 0 and 2. In order to set its value, we define three types of iterations: (i) Iteration E, which is the iteration with no improvement on the lower bound. A sequence of E iterations requires the need for a smaller step size. erefore, after a sequence of 20 E iterations, we multiply λ by 0.66.
ij , for all j and d � v t · w. If d < 0, this means that a larger step size would have given a smaller value for L(u t ), and we call this iteration Y. (iii) If d ≥ 0, we call this iteration T; a T iteration suggests the need for a larger step size, so we multiply λ by 1.1.

Computing LB.
In [12], we found that there are two approaches to solve the L(u) formulation of the problem above to obtain a LB. To keep this paper self-contained, we list the following approaches: e first approach is to decompose L(u) into m independent subproblems. Solve them as explained above. en, update this LB using the VA as mentioned previously, but in two ways for step 1. One way is to set the value of the vector u to 0 and continue using the VA to improve the LB. e other way is to set the vector u to values that is denoted by "warm start duals." e warm start duals are the values of the duals of the relaxed constraints (1), obtained from solving the greedy weak representation of the CFLP as a solution.
e second approach is first to remove constraints (11) from the L(u) formulation to reduce the size of the problem to make it possible to solve. en, solve it without the decomposition technique above, to obtain a LB. To improve this LB again, we apply the VA with the same two ways we discussed above.
Notice that the constraints (13) are included in the L(u) formulation, but only the part p ≤ i∈F y i is applicable in obtaining the LB as we explained before. e other part is redundant here since it only has an impact on obtaining the UB in the next approaches. e lower bound in this approach does not always have integer values y, while this is not the case in the decomposition approach; hence, the LB from this approach is often worse than the LB from the decomposition approach.

Computing UB.
To compute the UB, one way is to first remove constraints (3). Solve the P formulation for a LP solution. en, use a technique called randomized rounding (RR) with a new technique called the unit cost technique below, to treat the fractional solutions of y as a probability distribution and keep opening them randomly to get enough capacity. We keep updating this UB using the VA every 50 passes, until we meet one of the stopping criteria. We choose the passes to be 50 because the UB changes slowly, and from our experiments, 50 passes usually show some improvement in the UB. is is one way of obtaining the UB. e other way is similar to the above, except after we solve the P problem; we call the upper bound of the greedy weak representation. We call this strategy the "warm start UB." is UB results in a smaller step size of the VA. erefore, both the LB and the UB have converged faster in our experiments. en, after 50 passes, we call the RR technique below as before. We keep using the VA and the RR technique until we meet one of the stopping criteria.
A third technique called the "unit cost technique" (UC) selects the y ′ s to open according to their costs. It chooses the ones with the minimum unit cost. We call this technique once at the beginning before calling the RR. en after 50 passes, we switch to call the RR and keep calling the RR every 50 passes to update the UB until one of our stopping criteria is met. Next, we discuss these techniques in more detail.

Discussion.
Combining the approaches discussed above, which were used to obtain the LB and the UB, we got eight different approaches to solve the CFLP. Figure 1 below is a tree chart to simplify these eight approaches that are denoted by methods LR1-LR8. ese methods are used in solving different sizes of the CFLP. To decide which methods are used to solve which problem, we classify the problems into three classes: the small, the medium, and the large.

Computational Results
e selected model collection was drawn from Beasley [14], and we selected two instances where they fail to find feasible solutions when they use the factor � 1.5. e factor value is defined as i∈F S i / j∈D d j , which is used to rescale the capacity values. ey used heuristics that produce a feasible integer solution and used a Lagrangian relaxation and the volume algorithm to obtain a lower bound on the optimal value. For the CFLP, they were able to solve instances up to 1000 × 1000.
We use instances generated as in the work of Beasley [7] and the work of Ravi and Sinha [9], which are as follows: (i) We generate the demand points j and the facility points i uniformly at random in a  [10,160] distribution and are then rescaled so that i∈F S i / j∈D d j was set equal to the parameter factor; this was either set at 10 and in some cases 1. Next, we report the results of our experiments for two algorithms using a number of approaches denoted by methods LR1-LR8. e first are the approaches of the Lagrangian relaxation algorithm without decomposition (LR1-LR4), and we denote this algorithm by LR and VA. e second are the approaches of the Lagrangian relaxation with the decomposition algorithm (LR5-LR8), and we denote this algorithm by LD and VA.
Recall that, we classify the problems into three classes according to the size. e classes are small, medium, and large. In the next subsections, we report the results and the analysis of each class of the problems. Table 1 below, we report the solutions for solving the small problems of size 100-300. By comparing the solutions of the two algorithms explained above, we notice that both the LB and the UB of the Lagrangian decomposition with the volume algorithm LR5-LR8 are much better in both quality and running time.

e Small Size Problems. In
ere is a big difference in the LB between the two algorithms, and the reason is that one of the advantages of the LD algorithm is forcing y ′ s to be 1. In addition, we were able to 4 Advances

in Operations Research
Small and medium problems   Advances in Operations Research 5 obtain a closed LB as greedy's of the strong representation of the problem, but a better UB and running time. Also, the error values are less than greedy's in some of our methods LR5-LR8.
In the next tables, the column "passes" presents the number of times we go through the algorithm looking for a better LB and UB. e last row corresponds to the best greedy solution we were able to obtain; we include it in these tables for comparison purposes. e last column corresponds to the error percentage values, which can be obtained by (UB − LB)/UB. is value shows the percentage gap between the lower and the upper bounds of the problem. We use this value as a stopping criterion if it is less than 2%. Also, the stopping criteria here in most of the results are the number of passes, and in the rest, they are the error values.   Table 2 below lists all the results of solving the medium problems sized 400-900. We used the Lagrangian decomposition algorithm only because the size of the problem is too large to solve without decomposing them. e results of problem sized 500 × 500 shows that we were able to obtain a better LB and UB compared to those of greedy's solution. Also, the gap between them is less than the gap of greedy's solution. In problems sized 500-900, we were able to obtain a feasible integer solution (UB) within a gap of less than 1% of the lower bound LB in most of our solutions and in others less than 2.5%. e stopping criterion was the error value in most of the results.

e Large Size Problems.
In Table 3, we illustrate all the results of solving large-sized problems using the Lagrangian decomposition with the volume algorithm, which is the goal of this research. We were able to solve large instances of the CFLP. e largest problem solved was of size 3000 × 3000. In some of the problems, we were able to compare the solution with the greedy's solution. Table 3 shows that we obtained a better LB and UB than the greedy's. Also, the gap between them is much less than the gap between the UB and the LB of the greedy's. In all of the problem solution, the gap is less than 2.63%.

Conclusion
In this article, we have investigated solving very large instances of the CFLP. We first presented a general algebraic model for the CFLP. In addition, the created Lagrangian decomposition representation solves large instances of the CFLP. We have improved the lower bound of this technique by both exploiting the decomposition approach as well as introducing a tightening constraint, which is a new positive addition to the Lagrangian decomposition approach for solving the CFLP. In the case of very large problems, we have applied a randomized rounding technique along with our unit cost technique, which provided a good UB in a reasonable computational time. e thrust of this paper research was to be able to solve very large instances of the CFLP. We have illustrated that our algorithm scales up. Our results show that increasing the problem size leads to a small relative error between the LB and the UB without too much burden on the computational time. In order to process such large model instances, we have exploited a sparse technology technique in implementing our algorithm. us, we have developed a sparse data structure that we exploit in our algorithm implementation.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e author declares that there are no conflicts of interest.