Algorithmic Methods for Concave Optimization Problems

Xu Zhang¹, Xue Tian², Chen Wang², Tao Li³

¹Beijing Datang Dingwang Trading Company, Beijing, China.
²Beijing Wuzi University, Beijing, China.
³Lehigh University, Bethlehem, USA.

DOI: 10.4236/ajibm.2017.77067                

Abstract

In this thesis, we reformulate the original non-linear model for the LMRP. Firstly, we introduced a set of parameters to represent the non-linear part of the cost increase for a facility space allocated potential additional costs and new set of decision variables, indicating how many customers each equipment distribution. The algorithms are tested on problems with 5 to 500 potential facilities and randomly generated locations. Then using actual data to validate this new method is better. Our work was motivated by the modeling approach used in the Maximum Expected Covering Location Problem (MEXCLP). We compare new method and Lagrangian relaxation method to solve LMRP with constant customer demand rate and equal standard deviation of daily demand.

Keywords

LMRP, Lagrangian Method, Linearization Method

1. Introduction

Our work was motivated by the modeling approach used in the Maximum Expected Covering Location Problem (MEXCLP). MEXCLP is introduced by Mark S. Daskin in 1983 [1] [2] . The location model with risk pooling (LMRP) is introduced by Mark S. Daskin, Collette R. Coullard, and Zuo-Jun Max Shen in 2002 [3] . This problem chooses facility locations in order to minimize the total cost of building the facilities, transporting goods from facilities to customers and holding inventory to take advantage of economies of scale and protect against uncertain demand. So there are four parts in the objective function; the construction cost, the transportation cost, the cycle stock cost and the safety stock cost. This problem can be solved quite efficiently using Lagrangian relaxation when the ratio of the demand variance to mean is the same for every customer. In our thesis, we consider a further special case in which we assume all facilities share the same customer demand rate and standard deviation of daily demand. Ozsen [4] considered the interdependence between capacity and inventory management in the LMRP. The Lagrangian sub-problem is also a non-linear integer program. They proposed an efficient algorithm for the continuous relaxation of this sub-problem.

This model is a kind of covering problem; it decides the number of vehicles in each location in order to maximize the expected number of demands that can be covered, given that vehicles may be unavailable (in use). The model assumes that there is an equal probability that a vehicle is busy at any location. As the objective function is the expected number of demands, the decision variables that choose the number of vehicles in each location appear in an exponential term. This makes the objective function non-linear, just like the LMRP problem. Das- kin introduces a set of parameters to represent the increase in the expected coverage for each additional vehicle, as well as a set of binary decision variables to indicate whether the customer is covered a specific numbers of times. By using the sum of all th0e benefits of adding a new vehicle to represent the expected coverage, he changes the problem into one that is linear and easy to solve. So we apply the same idea to convert the LMRP into a linear mixed-integer programming problem and compare it with the Lagrangian method of Mark S. Daskin, Collette R. Coullard, and Zuo-Jun Max Shen to see if it will give us a more efficient method [5] [6] . Daskin and Teo [7] presented a stochastic version of the LMRP problem, and they developed a Lagrangian method for this problem. They also discussed the influence of changing the key parameters.

2. Model Formulation

2.1. Maximum Expected Covering Location Problem

As we mentioned in the Introduction chapter, our approach for solving concave binary minimization problems is inspired by a reformulation strategy that is sometimes used to solve other binary optimization problems in which the objective function contains a non-linear function of the sum of the binary variables. The basic idea is to introduce auxiliary parameters and binary variables and use their product to represent the none-linear part, and use these to linearize the objective function.

One model that uses this approach is the maximum expected covering location problem (MEX-CLP) by Mark S. Daskin in 1983. The MEXCLP chooses locations of facilities that can sometimes be unavailable (e.g., because the ambulance located there is busy on another call). A demand node is covered by a facility if it is within a certain coverage radius of it. The goal of the MEXCLP is to locate at most P facilities to maximize the total expected coverage of the demand nodes.

The MEXCLP assumes that the probability that a facility is unavailable at any time is given by q. It also assumes that facility unavailability are independent, so if there are n facilities that cover a demand node, then the probability that all of them are unavailable is given by qⁿ. Since the number of covering facilities, n, is not known a priori, we have to express it in terms of the decision variables as ${\displaystyle {\sum }_{j\in J}{a}_{ij}{X}_j}$ , where $a_{ij}$ is a parameter that equals 1 if facility j covers demand node i and 0 otherwise. Then the model can be formulated as follows:

Parameters:

$J$ set of potential facilities, indexed by j,

$I$ set of customer nodes, indexed by i,

$q$ the probability that a facility is unavailable at any time,

$P$ the maximum number of the facilities can be chosen,

$h_i$ the demand generated at node i,

$a_{ij}=\{\begin{array}{l}1\text{ifafacilityat}j\text{cancoverdemandsatcustomernode}i\\ \text{0}\text{otherwise}\end{array}$

Decision Variables;

$X_j$ the number of facilities to be built at j

Then the model can be formulated as follows:

$\text{Maximize}\underset{i\in I}{{\displaystyle \sum }}{h}_i\left(1-q^{{\displaystyle {\sum }_{j\in J}{a}_{ij}{X}_j}}\right)$

$\text{subject}\text{to}\underset{j\in J}{{\displaystyle \sum }}{X}_j\le P$

$X_j\in 0,1,\forall j\in J$

2.2. Linearization of Maximum Expected Covering Location Problem

In the original formulation, the probability that the demand of a customer node i is covered is given by $1-q^{{\displaystyle {\sum }_{j\in J}{a}_{ij}{X}_j}}$ , which is a non-linear function of X_j. Instead of computing the probability directly, proposes adding up the benefits of each new facility. We now summarize his approach.

Figure 1 shows how we compute the probability by adding up benefits. The x-axis is the number of facilities that cover the demand node and the y-axis is the probability. Assume there are k facilities that can cover the demand node. We first compute the bene t of adding the $\left(n+1\right)$ st facility (assuming we already have n facilities ) ( $n=0,1,\cdots ,P-1$ ) and then we can add them up from n = 1 to k.

The availability probability for n facilities is $\left(1-q^{n+1}\right)-\left(1-q^n\right)=q^n\left(1-q\right)$ . We introduce a new variable $Z_{jk}$ to represent the number of times covered, which we define to be 1 if demand node i is covered k or more times, and 0 if not. The model then can be formulated as follows.

$\text{Maximize}\underset{k=1}{\overset{P}{{\displaystyle \sum }}}\underset{i\in I}{{\displaystyle \sum }}{h}_i{q}^{k-1}\left(1-q\right)Z_{jk}$

$\text{subject}\text{to}\underset{j\in J}{{\displaystyle \sum }}{a}_{ij}{X}_j-\underset{k=1}{\overset{P}{{\displaystyle \sum }}}{Z}_{ik}\ge 0,\forall i\in I$

$X_j\in 0,1,\forall j\in J$

Figure 1. Adding up benefits to compute probability.

$Z_{jk}\in 0,1,\forall j\in J;k=1,2,\cdots ,P.$

In this model, we add up the benefits to replace the none-linear part of the objective function and that gives us a linear formulation. In what follows we propose a similar method to reformulate the LMRP model as a linear one.

2.3. The Location Model with Risk Pooling

The LMRP model is an extension of the UFLP that considers uncertain demand. Besides the fixed cost of opening locations and the variable transportation cost, it also includes the cost of cycle stock and safety stock. As a result, the LMRP is structured much like the UFLP model, with two extra non-linear terms in the objective function. Despite its concave objective function, the LMRP problem can be solved by Lagrangian relaxation quite efficiently, just like the UFLP, assuming that the ratio of the customer demand rate and the standard deviation of daily demand are constant. We use the following notations:

Parameters:

$I$ set of retailers, indexed by i,

$J$ set of candidate DC sites, indexed by j,

$u_i$ mean daily demand of retailer i, for each $i\in I$

${\sigma }_i^2$ variance of daily demand of retailer i, for each $i\in I$ ,

$f_j$ fixed (daliy) demand of locating a DC at candidate site j, for each $j\in J$ ,

$K_j$ fixed cost for DC j to place an order from the supplier, including fixed components of both ordering and transportation costs, for each $j\in J$ ,

$d_{ij}$ cost per unit to ship between retailer i and candiddate DC site j, for each $i\in I$ and $j\in J$

$\theta$ a constant parameter that captures the safety stock costs at candidate sites.

Decision Variables:

$X_j=\{\begin{array}{l}1\text{ifwelocateatcandidatesite}j\\ 0\text{if}\text{not}\end{array}$

$Y_{ij}=\{\begin{array}{l}\text{1}\text{ifdemandsatretailer}i\text{areassignedtoaDCatcandidatesite}j\\ 0\text{ifnot}\end{array}$

Then the model is formulated as follows.

$\text{Minimize}\underset{j\in J}{{\displaystyle \sum }}\left\{f_j+\underset{i\in I}{{\displaystyle \sum }}{d}_{ij}{Y}_{ij}+K_j\sqrt{\underset{i\in I}{{\displaystyle \sum }}{u}_i{Y}_{ij}}+\theta \sqrt{\underset{i\in I}{{\displaystyle \sum }}{\sigma }_i^2{Y}_{ij}}\right\}$

$\text{subjectto}\underset{j\in J}{{\displaystyle \sum }}{Y}_{ij}=1,\forall i\in I$

$Y_{ij}\le X_j,\forall i\in I,\forall j\in J$

$X_j\in \left\{0,1\right\},\forall j\in J$

$Y_{ij}\in \left\{0,1\right\},\forall i\in I,\forall j\in J$

2.4. Linearization of LMRP

To make the objective function linear, we introduce a new parameter ${\gamma }_{jk}$ to represent the cost of safety and cycle stock cost that k retailers are assigned to DC j, that is

${\gamma }_{jk}=K_j\sqrt{ku}+\theta \sqrt{k{\sigma }^2}$

Also we introduce a new decision variable

$Z_{jk}=\{\begin{array}{l}1,\text{ifexactly}k\text{retailersareassignedtoDC}j,\\ 0,\text{if}\text{not}\end{array}$

To associate $Z_{jk}$ with its meaning using linear constraints, we add the constraints

$\underset{k=0}{\overset{\left|I\right|}{{\displaystyle \sum }}}k{Z}_{jk}=\underset{i\in I}{{\displaystyle \sum }}{Y}_{ij},\forall j\in J$

$\underset{k=0}{\overset{\left|I\right|}{{\displaystyle \sum }}}{Z}_{jk}=1,\forall j\in J$

The second constraint says that only one of the $Z_{jk}$ can be equal to 1 for each j and the first constraint makes sure that the 1 appears when $k={\displaystyle {\sum }_{i\in I}{Y}_{ij}}$ which is just how we define the meaning of $Z_{jk}$ .

So the linear model is:

$\text{Minimize}{\displaystyle {\sum }_{j\in J}\left\{f_j{X}_j+{\displaystyle {\sum }_{i\in I}{d}_{ij}{Y}_{ij}}+{\displaystyle {\sum }_{k\in J}{\gamma }_{jk}{Z}_{jk}}\right\}}$ (2.1)

$\text{subject}\text{to}{\displaystyle {\sum }_{j\in J}{Y}_{ij}}=1,\forall i\in I$ (2.2)

$Y_{ij}\le X_j,\forall i\in I,\forall j\in J$ (2.3)

${\displaystyle {\sum }_{k=0}^{\left|I\right|}k{Z}_{jk}}={\displaystyle {\sum }_{i\in I}{Y}_{ij}},\forall j\in J$ (2.4)

${\displaystyle {\sum }_{k=0}^{\left|I\right|}{Z}_{jk}}=1,\forall j\in J$ (2.5)

$X_j\in \left\{0,1\right\},\forall j\in J$ (2.6)

$Y_{ij}\in \left\{0,1\right\},\forall i\in I,\forall j\in J$ (2.7)

$Z_{jk}\in \left\{0,1\right\},\forall j\in J,\forall k=0,\cdots ,I$ (2.8)

From these two formulations, we can see although the second method is linear, it has many more constraints than the original formulation. On the other hand, it can be solved by an o―the-shelf MIP solver and does not require Lagrangian relaxation as in the original LMRP. So it’s hard to say which computation time would be shorter only by looking at the models. We will test randomly generated examples and compare the solution time of the two methods in Chapter 4.

2.5. The Lagrangian Relaxation Method for the LMRP

Similar to the UFLP, we solve the LMRP by relaxing the assignment constraints Equation (2.2) to obtain the following Lagrangian sub-problem:

$\begin{array}{l}\text{Minimize}\underset{j\in J}{{\displaystyle \sum }}\left\{f_j+\underset{i\in I}{{\displaystyle \sum }}{d}_{ij}{Y}_{ij}+K_j\sqrt{\underset{i\in I}{{\displaystyle \sum }}{u}_i{Y}_{ij}}+\theta \sqrt{\underset{i\in I}{{\displaystyle \sum }}{\sigma }_i^2{Y}_{ij}}\right\}\end{array}$ $+\underset{i\in I}{{\displaystyle \sum }}{\lambda }_I\left(1-\underset{j\in J}{{\displaystyle \sum }}{Y}_{ij}\right)=\text{Minimize}\underset{j\in J}{{\displaystyle \sum }}\{f_j+\underset{i\in I}{{\displaystyle \sum }}\left(d_{ij}-{\lambda }_i\right)Y_{ij}$ $+K_j\sqrt{\underset{i\in I}{{\displaystyle \sum }}{u}_i{Y}_{ij}}+\theta \sqrt{\underset{i\in I}{{\displaystyle \sum }}{\sigma }_i^2{Y}_{ij}}\}+\underset{i\in I}{{\displaystyle \sum }}{\lambda }_i$

$\text{subject}\text{to}{Y}_{ij}\in X_j,\forall i\in I,\forall j\in J$

$X_j\in \left\{0,1\right\},\forall j\in J$

$Y_{ij}\in \left\{0,1\right\},\forall i\in I,\forall j\in J$

Although the sub-problem is a concave integer minimization problem, it can be solved relatively efficiently, using a sorting method developed by Mark S. Daskin, Collette R. Coullard, and Zuo-Jun Max Shen in 2003. The algorithm relies on the assumption that the ratio of the demand variance to the demand mean is a constant for all retailers. That is, for all $i\in I$ , ${\sigma }_i^2/u_i=\gamma \le 0$ . Then we can collapse two square root terms into one and apply the sorting algorithm to solve the resulting sub-problem.

The optimal objective function value of the Lagrangian sub-problem gives us a lower bound of the original problem; then we need an upper bound. There are many ways to find a feasible solution to get the upper bound; in this paper, we use a simple algorithm to generate the solution from the sub-problem result. This is shown in the appendix.

Finally, we recursively update $\lambda$ to get a smaller gap between the lower and upper bound. Our stopping condition in the computational tests in this thesis is when the number of iterations is over 500 or the gap is less than or equal to 5 percent of the upper bound. There is no limit for CPU time since the first stopping condition includes it.

3. Computational Results

3.1. Testing on Random Instances

We implemented the Lagrangian method in C++ and the linearization method in AMPL with CPLEX version 12.4.0.0. Table 1 is the comparison of the solution time for the linearization and Lagrangian methods. The linearization method has a similar solution time as the Lagrangian method when the number of facilities is small. However, the solution time increases faster with the number of

Table 1. Average solution time comparison.

facilities for the linearized method than it does for the Lagrangian method. So for larger scale problems (which are more practical) the Lagrangian method will have better performance.

In our experiments, for one specific data set, CPLEX gets stuck when it tries to solve the linearization problem. It takes over 90 seconds while the other samples with the same data scale only need 2.91 seconds on average. When we use the Lagrangian method to solve the same problem, the method also stopped because the number of iterations is over 500. The reason that the Lagrangian method can’t solve this kind of data set in a small number of iterations is that the Lagrangian relaxation’s optimal value can’t reach the original problem’s optimal value, and the gap is over 5 percent, but the reason why the gaps are large for this data set is not clear. Similarly, we still can’t understand why CPLEX also gets stuck for this data set.

The gap between the upper bound and the lower bound is not large; it is on average 4.3 percent. However, there is a significant cant increase in the gap when the data scale grows larger.

Figures 2-4 and Table 1 are the comparison of the average solution time for the Lagrangian method and the linearization method. In the figure, the lower line (blue) represents the solution time for the Lagrangian method and the upper line (orange) is for the linearization method.

Figure 5 and Figure 6 and Table 2 and Table 3 are the comparison of the maximum solution time for the Lagrangian method and the linearization method.

Figure 2. Average solution time for data scale 5 to 100.

Figure 3. Average solution time for data scale 150 to 450.

Figure 4. Bounds from Lagrangian method.

Figure 5. Maximum solution time for data scale 5 to 100.

Figure 6. Maximum solution time for data scale 150 to 450.

Table 2. Maximum solution time comparison.

Table 3. Data distribution for LMRP.

Based on the results, we can see that the gap between the lower and upper bounds is acceptable, even in the cases that stopped due to the 500-iterations limit. Additionally, when the scale is not large, we can get results from the Lagrangian method with a tiny gap, say 0.5 percent. So it can be concluded that the Lagrangian method is reliable.

The linearization method becomes slower than the Lagrangian method on average when the data scale is large. From Table 4 and Table 5, we see that the maximum solution time of the ten samples has a similar trend as the average solution time. However, there do exist in-stances for which the linearization method runs faster. The solve-speed depends on the instances.

For further research, we will test on larger scale problems and real world instances.

3.2. Testing on Benchmark Instances

As in real life, the data such as distance, fixed cost of opening a new facility and the demand of different places are not always independent, so it is necessary to compare the computation time not only on random data sets, but also on examples that come from more realistic instances.

Our data comes from “An inventory-location model: Formulation, solution algorithm and computational results. Annals of Operations Research” and we use two data sets. For the 88-node dataset, representing the 50 largest cities in the 1990 US census along with the 48 capitals of the continental US minors duplicates, the mean demand was obtained by dividing the population data by 1000 and rounding the result to the nearest integer. Fixed facility location costs were obtained by dividing the facility location costs by 100. For the 150-node dataset, representing the 150 largest cities in the continental US for the 1990 census,

Table 4. Sample test data for LMRP.

Table 5. Sample test data for LMRP.

the mean demand was obtained in the same manner. The fixed facility costs were all set to 100, one thousandth of the value in the dataset given by Mark S. Daskin, Collette R. Coullard, and Zuo-Jun Max Shen in 2002. These changes were made to allow us to deal with smaller numbers.

For the 88-node dataset, the solution time for the Lagrangian method is 0.203 s and it takes CPLEX 2.435 s. For the 150-node dataset, the solution time for the Lagrangian method is 0.539 s and it takes CPLEX 19.673 s.

We see that the solution time for both methods is a little bit smaller than the average of the random samples and the Lagrangian method is still much faster han the linearization method. So the randomness of the initial instances may not have much influence on the comparison of these two methods.

4. Conclusions

Our linearization of the LMRP requires longer solution time on average than the Lagrangian method does. However, it performs better in some special instances.

For future research, we would like to determine under what conditions the linearization method will have a shorter solution time than the Lagrangian method does.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Daskin, M.S. (1982) Application of an Expected Covering Model to Emergency Medical Service System Design. Decision Sciences, 13, 416-439. https://doi.org/10.1111/j.1540-5915.1982.tb00159.x
[2]	Daskin, M.S. (1983) A Maximum Expected Covering Location Model: Formulation, Properties and Heuristic Solution. Transportation Science, 17, 48-70. https://doi.org/10.1287/trsc.17.1.48
[3]	Daskin, M.S., Coullard, C.R. and Shen, Z.-J.M. (2002) An Inventory-Location Model Formulation, Solution Algorithm and Computational Results. Annals of Operations Research, 110, 83-106. https://doi.org/10.1023/A:1020763400324
[4]	Ozsen, L., Coullard, C.R. and Daskin, M.S. (2008) Capacitated Warehouse Location Model with Risk Pooling. Naval Research Logistics, 55, 295-312. https://doi.org/10.1002/nav.20282
[5]	Shen, Z.-J.M. (2007) Integrated Supply Chain Design Models: A Survey and Future Research Directions. Journal of Industrial and Management Optimization, 3. https://doi.org/10.3934/jimo.2007.3.1
[6]	Shen, Z.-J.M., Coullard, C.R. and Daskin, M.S. (2003) A Joint Location-Inventory Model. Transportation Science, 37, 40-50. https://doi.org/10.1287/trsc.37.1.40.12823
[7]	Snyder, L.V., Daskin, M.S. and Teo, C.-P. (2007) The Stochastic Location Model with Risk Pooling. European Journal of Operational Research, 179, 1221-1238. https://doi.org/10.1016/j.ejor.2005.03.076