A semi-Lagrangian relaxation heuristic algorithm for the simple plant location problem with order

Abstract

The Simple Plant Location Problem with Order (SPLPO) is a variant of the Simple Plant Location Problem (SPLP) where the customers have preferences over the facilities that will serve them. In particular, customers define their preferences by ranking each of the potential facilities from the most preferred to the least preferred. Even though the SPLP has been widely studied in the literature, the SPLPO, which is a harder model to deal with, has been studied much less and the size of instances that can be solved is very limited. In this article, we study a Lagrangian relaxation of the SPLPO model and we show that some properties previously studied and exploited to develop efficient SPLP solution procedures can be extended for their use in solving the SPLPO. We propose a heuristic method that takes a Lagrangian relaxation solution as the starting point of a semi-Lagrangian relaxation algorithm. Finally, we include some computational studies to illustrate the good performance of this method, which quite often ends with the optimal solution.

Keywords:

• Simple plant location problem
• preferences
• Lagrangian relaxation
• semi-Lagrangian relaxation

1. Introduction

The Simple Plant Location Problem with Order (SPLPO) is a variant of the Simple Plant Location Problem (SPLP) where the customers have preferences over the facilities that will serve them. In particular, customers define their preferences by ranking each of the potential facilities from the most preferred to the least preferred. This kind of model is highly relevant in situations where the locator, the agent that decides where to locate the facilities, has no control over the customers to force them to attend the facility that is the most convenient from the optimised point of view of the locator. Therefore, this kind of problem can also be seen as a Stackelberg game or as a bilevel model (Kochetov et al., 2022). Both of them are different approaches for the SPLPO to the one used in this article and will not be discussed here.

In the public sector, some types of service centers must be located in strategic locations to distribute goods and make them available to target populations. The objective of the decision maker is to maximise some measure of the social welfare that consumers will after having access to the resources whereas users are interested in minimising things such as total distance travelled. As a consequence, the goals of the locator and the goals of the users may point to different directions. Another application can be found in the context of business-to-consumer (B2C) e-commerce. As mentioned in H. Zhang et al. (2019), e-commerce has logistics advantages over traditional B2C. However, although costs can be reduced considerably due to the integration of operations with smart technologies, there is still the problem of distributing the demanded goods to the customers from the distribution centers. Sometimes customers collect the products from these centers and therefore may have preferences about where to go. Ease of access or travel times are examples of variables that can influence this decision.

1.1. Literature review, gap, and contributions

The SPLPO was first introduced by Hanjoul and Peeters (1987), where the authors provided a linear formulation for this problem and suggested a heuristic procedure for very small instances. The main results found in the literature are on valid inequalities to strengthen the original formulation. Cánovas et al. (2007) provided a new family of valid constraints that were used in combination with a preprocessing analysis. Vasil’yev, Klimentova, and Kochetov (2009) studied a bilevel formulation and considered several reductions of the bilevel problem to integer linear programs. A branch-and-cut method was proposed by Vasilyev and Klimentova (2010). Another more general family of valid inequalities can be found in Vasilyev et al. (2013) with a polyhedral study. More recent works are Calvete et al. (2020), Casas-Ramírez & Camacho-Vallejo (2017), and B. Zhang et al. (2022), where very closely related problems that included preferences are studied.

Since most of these papers use exact methods that struggle to solve large instances, there is a clear gap in the literature regarding methods to solve the SPLPO more efficiently. One of the main contributions of this article is to develop a procedure to solve the SPLPO efficiently in a heuristic way by using Lagrangian and semi-Lagrangian relaxation techniques. This idea was inspired by the success of these methods in solving the SPLP, but it had never been used before to solve the SPLPO, a much more complex and difficult problem. For example, already in Cornuéjols et al. (1977) an application of Lagrangian relaxation was presented to solve the SPLP in the context of the location of bank accounts. They proposed a heuristic algorithm and studied the Lagrangian dual to obtain lower and upper bounds. They also provided a bound for the relative error of these methods. Beltrán, Tandoki, and Vial (2006) suggested, defined, and applied the technique of semi-Lagrangian relaxation to the p-median problem. Some years later, the study was extended to the SPLP obtaining very good results (Beltrán et al., 2012). The methodology takes advantage of the linear formulation of the problem. First, it splits equality constraints Ax = b into $Ax\le b$ and $Ax\ge b,$ and then relaxes the latter. The new model has the same objective value as the original problem (it closes the duality gap), but at the cost of making the new problem more difficult to solve. However, it has some good properties that can be exploited. Details will be given later in the article. Another reduced form of this method was proposed by Monabbati (2014), who called it a surrogate semi-Lagrangian relaxation. A new algorithm for the dual problem using Lagrangian heuristics for both the original and the surrogated version of the semi-Lagrangian relaxation can be found in (Jörnsten & Klose, 2016). Lagrangian and semi-Lagrangian approaches have been successfully used in other problems such as the optimal diversity management problem (ODMP) (Masone et al., 2019), the quadratic assignment problem (H. Zhang et al., 2016) and more recently for determining the locations of uncapacitated distribution centers for B2C e-commerce (H. Zhang et al., 2019). A literature review summarising the theoretical and practical advances for the SPLP, SPLPO, and other related problems with our work is shown in Table 1.

Table 1. Literature review table.

Year	Authors	Problem	Summary points
1977	Cornuéjols et al.	SPLP	An application of Lagrangian relaxation and a heuristic algorithm. The article provides lower and upper bounds on the value of the optimal solution with a worst-case analysis of these bounds.
1987	Hanjoul and Peeters	SPLPO	A linear formulation and a greedy heuristic for small instances.
2006	Cánovas et al.	SPLPO	A new family of valid constraints and a preprocessing analysis. Extensive computational experiments are carried out.
2006	Beltrán et al.	p-median	The semi-Lagrangian relaxation theory is introduced.
2009	Vasil’ev et al.	SPLPO	A bilevel facility location problem with preferences order is studied. New valid inequalities and lower bounds are introduced.
2010	Vasil’ev and Klimentova	SPLPO	Cutting plane method based on valid inequalities from Vasil’ev et al. (2009).
2012	Beltrán et al.	SPLP	Computational experiments on large-scale instances of the SPLP with unknown optimal values by using semi-Lagrangian relaxation approach combined with CPLEX.
2013	Vasil’ev et al.	SPLPO	A polyhedral study and a new family of facet-defining inequalities.
2014	Monabbati	SPLP	Surrogate version of semi-Lagrangian relaxation is introduced with an application to the SPLP.
2017	Casas-Ramírez and Camacho-Vallejo	p-median with order	Study of two different formulations. A hybrid heuristic algorithm based on scatter search (Laguna & Martí, 2006) and GRASP (Resende & Ribeiro, 2016) is proposed.
2019	H. Zhang et al.	SPLP	Mixed integer programming model for determining uncapacitated distribution centers location for B2C e-commerce. Model solved with semi-Lagrangian relaxation.
2019	Masone et al.	ODMP	ODMP is formulated as a large-scale p-median problem. A three stages method that uses Lagrangian relaxation, variable fixing and a dynamic programming is presented.

The SPLP has been widely studied in the literature and even today its study is intensive due to the different applications that can be considered, see for example (Souto et al., 2021). However, the SPLPO has been studied much less and the size of the instances that can be solved is much more limited. One of the reasons for this is that models that incorporate preferences are more difficult to deal with: the usual approaches are either to include a quadratic number of constraints that incorporate the preferences or to use bilevel optimisation. As the success of the second methodology has been even more limited than for the first one, our research takes the standard SPLPO formulation as the starting point (see Section 2). In this work, we follow the line of research given in (Cornuéjols et al., 1977) and (Beltrán et al., 2006) for the SPLP and the p-median problem, respectively, by studying for the first time in the literature some properties and advantages of using a Lagrangian and semi-Lagrangian approach for the SPLPO. This is exploited later when we propose to use the solution of a Lagrangian relaxation algorithm as the starting point of a semi-Lagrangian relaxation method, and then to apply a variable fixing heuristic to find good feasible solutions (often the optimal solution).

The rest of the article is organised as follows: The formulation for the SPLPO is given in Section 2. Section 3 provides a review of the Lagrangian and semi-Lagrangian relaxation methods. Our contributions explaining how to apply them to the SPLPO are explained in Sections 4 and 5. The complete method that we propose is presented in Section 6, where all the algorithms proposed in previous sections will be combined to build a heuristic procedure. In Section 7, we run a computational study that shows the good performance of our method. Finally, some conclusions are given in Section 8.

2. SPLPO formulation

Let $I=\{1,\dots m\}$ be a set of customers, and let $J=\{1,\dots n\}$ be a set of sites where facilities can be opened. We have costs $c_{ij}\ge 0$ for supplying the demand of customer i from facility j and costs $f_j\ge 0$ for opening a facility at site j. Each customer ranks fully the list of sites where facilities can be located from the most preferred to the least preferred and no ties are allowed in this ranking. We say that facility k is i-worse than facility j if customer i prefers facility j to facility k and it is written as $k{<}_{i}j.$ We define $W_{ij}=\{k\in J|k{<}_{i}j\}$ as the set of facilities strictly i-worse than j, its complement as ${\overline{W}}_{ij}$ (i-better than or equal to j,) and $W_{ij}^{\prime }:=W_{ij}\cup \left\{j\right\}.$

Let x_ij be a decision variable that represents the fraction of the demand required by customer i and supplied by facility j. Since no capacities are considered, this demand can always be covered completely by one single facility. Therefore, we can guarantee that there is an optimal solution where the values for variables x_ij are in {0, 1}. Let y_j be a binary variable such that y_j = 1 if a facility is open at location j and y_j = 0 otherwise.

The SPLPO formulation is as follows: (1) $$\begin{array}{c}\mathrm{Min}.\sum _{i\in I}\sum _{j\in J}{c}_{ij}{x}_{ij}+\sum _{j\in J}{f}_j{y}_j\\ \mathrm{s}.\mathrm{t}.\sum _{j\in J}{x}_{ij}=1,\hspace{1em}i\in I,\end{array}$$(1) (2) $$x_{ij}\le y_j,\hspace{1em}i\in I,\forall j\in J,$$(2) (3) $$\sum _{k\in {\overline{W}}_{ij}}{x}_{ik}\ge y_j,\hspace{1em}i\in I,j\in J,$$(3) (4) $$x_{ij}\ge 0,\hspace{1em}i\in I,j\in J,$$(4) (5) $$y_j\in \{0,1\},\hspace{1em}j\in J.$$(5)

Constraint (1)(1) $$\begin{array}{c}\mathrm{Min}.\sum _{i\in I}\sum _{j\in J}{c}_{ij}{x}_{ij}+\sum _{j\in J}{f}_j{y}_j\\ \mathrm{s}.\mathrm{t}.\sum _{j\in J}{x}_{ij}=1,\hspace{1em}i\in I,\end{array}$$(1) is the assignment constraint and ensures that every customer i will be supplied by exactly one facility. Constraint (2)(2) $$x_{ij}\le y_j,\hspace{1em}i\in I,\forall j\in J,$$(2) imposes that customers can only be serviced by open facilities. This family of constraints is usually called variable upper bounds (VUBs). Inequality (3)(3) $$\sum _{k\in {\overline{W}}_{ij}}{x}_{ik}\ge y_j,\hspace{1em}i\in I,j\in J,$$(3) models the customers’ preferences, that is, if a facility j is open, every customer i will be served by either j or by a facility which i prefers to j (an i-better facility in the terms that we defined earlier).

3. Lagrangian and semi-Lagrangian relaxation background

Let P be a problem of the form: $$\min \{f(x)=c^{t}x/Ax\le b,Cx\le d,x\in X\},$$ where:

• $c\in {\mathbb{R}}^n,x\in {\mathbb{R}}^n,A\in {\mathbb{R}}^{m_1\times n}$ and $C\in {\mathbb{R}}^{m_2\times n}.$

• Set X may contain integrality constraints.

• The family of constraints $Ax\le b$ is assumed to be complicated, that is, problem P without them is much easier to solve.

Constraint $Ax\le b$ can be moved to the objective function with coefficients vector (Lagrangian multipliers) λ that work as penalties when they are violated. It yields the Lagrangian relaxation problem $LR(\lambda )$ related to P with multipliers $\lambda \ge 0$ and with optimal value $v(LR(\lambda ))$: $$v(LR(\lambda )):=\min \{f(x)+{\lambda }^t(Ax-b)/Cx\le d,x\in X\}.$$

We have also the associated Lagrangian dual problem LD with optimal value v(LD): $$v(LD):=\underset{\lambda \ge 0}{\max }\min \{f(x)+{\lambda }^t(Ax-b)/Cx\le d,x\in X\}=\underset{\lambda \ge 0}{\max }v(LR(\lambda )).$$

Let v(P) be the optimal objective value for problem P. It is well-known that $v(P_L)\le v(LD)\le v(P),$ where P_L is the linear relaxation of problem P. Consider now the problem $P\prime$: $$\min \{f(x)=c^{t}x/Ax=b,Cx\le d,x\in X\},$$ where

• $c\in {\mathbb{R}}^n,x\in {\mathbb{R}}^n,A\in {\mathbb{R}}^{m_1\times n},$ and $C\in {\mathbb{R}}^{m_2\times n}.$

• Set X includes the integrality constraints.

• A, b, and c are nonnegative.

The semi-Lagrangian relaxation problem $\mathit{SLR}(\lambda )$ related to $P\prime$ and its optimal value $v(\mathit{SLR}(\lambda ))$ are defined as follows: $$v(\mathit{SLR}(\lambda )):=\min \{f(x)+{\lambda }^t(b-Ax)/Ax\le b,Cx\le d,x\in X\}.$$

After having split Ax = b into the two inequalities $Ax\le b$ and $Ax\ge b,$ inequality $Ax\ge b$ has been relaxed with a vector of multipliers $\lambda \ge 0$ and inequality $Ax\le b$ has been kept as a constraint. Its semi-Lagrangian dual problem and optimal value are then $$v(\mathit{SLD}):=\underset{\lambda \ge 0}{\max }v(\mathit{SLR}(\lambda )).$$

It is clear that $v(LR(\lambda ))\le v(\mathit{SLR}(\lambda ))$ since $LR(\lambda )$ is a relaxation of SLR(λ). Thus, $v(P{\prime }_L)\le v(LD)\le v(\mathit{SLD})\le v(P\prime ).$ Beltrán et al. (2006) proved that the semi-Lagrangian dual problem has the same optimal value as $P\prime ,$ that is, $v(P\prime )=v(\mathit{SLD}).$ They also proved that the objective function of $\mathit{SLR}(\lambda )$ is concave, non-decreasing on its domain, and that b—Ax is a subgradient at point λ. Moreover, there is a λ with each component i in an interval $[{\lambda }_i^{*},+\infty ),$ where we obtain the optimal solution of $\mathit{SLR}(\lambda ).$ For details, see Geoffrion (1974), Fisher (2004), Beltrán et al. (2006), and Beltrán et al. (2012).

4. A Lagrangian relaxation method for the SPLPO

In this section, we consider a Lagrangian relaxation for the SPLPO. Our first result shows that the proposed model is easy to solve.

If we relax constraints (1)(1) $$\begin{array}{c}\mathrm{Min}.\sum _{i\in I}\sum _{j\in J}{c}_{ij}{x}_{ij}+\sum _{j\in J}{f}_j{y}_j\\ \mathrm{s}.\mathrm{t}.\sum _{j\in J}{x}_{ij}=1,\hspace{1em}i\in I,\end{array}$$(1) and (3)(3) $$\sum _{k\in {\overline{W}}_{ij}}{x}_{ik}\ge y_j,\hspace{1em}i\in I,j\in J,$$(3) , then they are moved to the objective function with penalty coefficients (multipliers) when they are violated, thus obtaining the following Lagrangian relaxation problem $LR(\mu ,\lambda ):$ $$\begin{array}{c}{\min }_{(x,y)}\sum _{i,j}{c}_{ij}{x}_{ij}+\sum _j{f}_j{y}_j+\sum _i{\mu }_i\left(1-\sum _j{x}_{ij}\right)+\sum _{i,j}{\lambda }_{ij}\left[y_j-\sum _{k\in {\overline{W}}_{ij}}{x}_{ik}\right]=\hfill \\ {\min }_{(x,y)}\sum _{i,j}\left(c_{ij}-{\mu }_i\right)x_{ij}-\sum _{i,j}{\lambda }_{ij}\sum _{k\in {\overline{W}}_{ij}}{x}_{ik}+\sum _j\left(f_j+\sum _i{\lambda }_{ij}\right)y_j+\sum _i{\mu }_i\hfill \\ \mathrm{s}.\mathrm{t}.(2),(4),(5).\hfill \end{array}$$

The multiplier vectors μ and λ are unrestricted and nonnegative, respectively.

To obtain the best possible lower bound with this relaxation, we need to solve the following Lagrangian dual problem: $$LD=\underset{\mu \mathrm{\hspace{0.33em}}\text{free},\lambda \ge 0}{\max }LR(\mu ,\lambda ).$$

Suppose that each customer i ranks each facility j with a number $p_{ij}\in \{1,\dots ,n\}$ with 1 and n the most and the least preferred facilities, respectively. Since each multiplier λ_ij in $\sum _i\sum _j{\lambda }_{ij}\sum _{k\in {\overline{W}}_{ij}}{x}_{ik}$ in $LR(\mu ,\lambda )$ will be multiplied by a sum of p_ij variables x_ik with $k{\ge }_{i}j,$ then each x_ij will be multiplied by a sum of $n-p_{ij}+1$ values λ_ik with $k{\le }_{i}j.$ Therefore, $\left|W{\prime }_{ij}\right|=n-p_{ij}+1$ and $$\sum _{i,j}{\lambda }_{ij}\sum _{k\in {\overline{W}}_{ij}}{x}_{ik}=\sum _{i,j}\left(\sum _{k\in W{\prime }_{ij}}{\lambda }_{ik}\right)x_{ij}.$$

Then $LR(\mu ,\lambda )$ can be written as follows: $$\begin{array}{c}{\min }_{(x,y)}\sum _{i,j}(c_{ij}-{\mu }_i-\sum _{k\in W{\prime }_{ij}}{\lambda }_{ik})x_{ij}+\sum _j\left(f_j+\sum _i{\lambda }_{ij}\right)y_j+\sum _i{\mu }_i\hfill \\ \mathrm{s}.\mathrm{t}.(2),(4),(5).\hfill \end{array}$$

This problem is easy to solve analytically for fixed vectors μ and λ. If we define Λ_ij as $${\Lambda }_{ij}:=\sum _{k\in W{\prime }_{ij}}{\lambda }_{ik},$$ then we have the following result:

Theorem 4.1.

An optimal solution for $LR(\mu ,\lambda )$ can be obtained as follows:$$y_j=\{\begin{array}{cc}1,\hfill & \mathrm{ }\mathit{\text{if}}\mathrm{ }\sum _i\min (0,c_{ij}-{\mu }_i-{\Lambda }_{ij})+f_j+\sum _i{\lambda }_{ij}<0,\hfill \\ 0,\hfill & \mathit{\text{otherwise}},\hfill \end{array}$$and$$x_{ij}=\{\begin{array}{cc}1,\hfill & \mathit{\text{if}}\mathrm{ }{y}_j=1\mathrm{ }\mathit{\text{and}}\mathrm{ }{c}_{ij}-{\mu }_i-{\Lambda }_{ij}<0,\hfill \\ 0,\hfill & \mathit{\text{otherwise}}.\hfill \end{array}$$

Proof.

The solution to $LR(\mu ,\lambda )$ is straightforward. □

An issue with this model is that $LR(\mu ,\lambda )$ has the integrality property, that is, its optimal value is equal to the value of the linear relaxation of the original problem. Furthermore, the solutions obtained for $LR(\mu ,\lambda )$ during the search for a solution to LD may be infeasible for the SPLPO formulation. So, as an alternative, we propose to use the solution of this problem as a starting point for a semi-Lagrangian procedure that allows us to find feasible solutions to the SPLPO. Before introducing this method, we provide first the details of the subgradient method used to solve the Lagrangian problem that has been formulated.

4.1. Subgradient method for the Lagrangian dual

The subgradient method was originally proposed by Held and Karp (1971) and validated by Held et al. (1974). Given vectors of multipliers λ and μ, this method tries to optimise the Lagrangian dual problem LD by taking steps along a subgradient of $LR(\mu ,\lambda ).$ A sketch of the whole procedure is given in Algorithm 1.

Algorithm 1:

Subgradient Method (SGM)

1 begin

   input: An SPLPO problem.

   output: A lower bound for SPLPO, multipliers $(\mu ,\lambda ).$

2  $\mathtt{iter}:=0;$

3  $\mathtt{MAXiter}\in \mathbb{N};$

4  Compute an upper bound LR_ub for $LR(\mu ,\lambda )$ by applying a heuristic procedure;

5  $\beta :=2;$

6  $\mathtt{k}\in \mathbb{N};$

7  $\mathtt{q}\in (0,1)\subset \mathbb{R};$

8  $[{\mu }^{\mathtt{iter}},{\lambda }^{\mathtt{iter}}]:=[{\min }_j\{c_{ij}+f_j\},0];$

9  $L{R}_{\mathit{best}}^{\mathtt{iter}}:=LR({\mu }^{\mathtt{iter}},{\lambda }^{\mathtt{iter}});$

10  repeat

11   Find a subgradient $s^{\mathtt{iter}}=(s_{\mu }^{\mathtt{iter}},s_{\lambda }^{\mathtt{iter}})$ for $LR({\mu }^{\mathtt{iter}},{\lambda }^{\mathtt{iter}});$

12   if $\Vert s^{\mathtt{iter}}\Vert \ge {10}^{-6}$ then

13    ${\alpha }^{\mathtt{iter}}:=\beta \frac{L{R}_{ub}-LR({\mu }^{\mathtt{iter}},{\lambda }^{\mathit{iter}})}{\Vert s^{\mathtt{iter}}|{|}_2^2};$

14   end

15   $[{\mu }^{\mathit{iter}+1},{\lambda }^{\mathit{iter}+1}]:=[{\mu }^{\mathit{iter}}-{\alpha }^{\mathit{iter}}{s}_{\mu }^{\mathit{iter}},\max (0,{\lambda }^{\mathit{iter}}-{\alpha }^{\mathit{iter}}{s}_{\lambda }^{\mathit{iter}})];$

16   $L{R}_{\mathit{best}}^{\mathtt{iter}+1}:=\min \{L{R}_{\mathit{best}}^{\mathtt{iter}},LR({\mu }^{\mathtt{iter}+1},{\lambda }^{\mathtt{iter}+1})\};$

17   $\mathtt{iter}:=\mathtt{iter}+1;$

18   if $LR({\mu }^{\mathtt{iter}},{\lambda }^{\mathtt{iter}})$ does not improve after $\mathtt{k}$ iterations then

19    $\beta =\beta \times \mathtt{q};$

20   end

21  until Either $\Vert s^{\mathtt{iter}}\Vert <{10}^{-6}$ or $\mathtt{iter}=\mathtt{MAXiter};$

22 end

In our computational experiments, which will be showed later, an upper bound LR_ub for $LR(\mu ,\lambda )$ on line 1 of Algorithm 1 is found with another heuristic (see Algorithm 2). As seen on line 1 of Algorithm 1, at each iteration we need to solve an instance of $LR(\mu ,\lambda ),$ which we know how to do (see Theorem 4.1).

It is possible, due to an overestimation of LR_ub, that the value $LR(\mu ,\lambda )$ does not improve for several iterations. This can be overcome by setting the parameter β to an initial value (e.g., 2) and reducing it slowly by a factor $\mathtt{q}\in (0,1).$ Details on the step size formula on line 1 can be found in Poljak (1967), Conforti et al. (2014), or Guignard (2003).

Regarding Algorithm 2, the heuristic to find a feasible solution to the original problem, first we open a facility $j^{*}$ that supplies all customers with the lowest operating cost, and it is removed from set $J\prime :=J.$ Then, we do iteratively the following greedy procedure. Among the closed facilities, we open the one that would lead to the best objective function value (smallest cost) when we open that facility and reallocate customers to that facility if it is more preferred than their current allocation. That facility is removed from $J\prime$ and we repeat the procedure, now with one less facility in $J\prime ,$ until $J\prime$ is empty. As can be seen from computational experiments in Section 7, the heuristic in Algorithm 2 works quite well in terms of efficiency and effectiveness.

Algorithm 2:

Upper Bound Heuristic for SPLPO (UBH)

1 begin

  input: A bipartite graph $G(I\cup J,E:=\{\})$ and costs c_ij.

   output: An upper bound for the SPLPO.

2  Find $j^{*}\in J$ such that $\sum _i{c}_{i{j}^{*}}=\min \{\sum _i{c}_{i1},\dots ,\sum _i{c}_{in}\};$

3  $J\prime :=J\setminus \left\{j^{*}\right\};$

4  $E_{j^{*}}:=\{(i,j^{*})/i\in I\};$

5  $E_r:=E_{j^{*}};$

6  $T{C}_{j^{*}}:=\sum _i{c}_{i{j}^{*}};$

7  repeat

8   for $j\in J\prime$ do

9    $E_j:=\left\{\right\};$

10    for $i\in I$ do

11     Find k_pref such that $p_{i{k}_{\mathit{pref}}}=\min \{\{p_{ij}\}\cup {\{p}_{ik}/(i,k)\in E_r\}\};$

12     $E_j:=E_j\cup \{(i,k_{\mathit{pref}})\};$

13    end

14    Compute $T{C}_j:=\sum _{(i,j)\in E_j}{c}_{ij};$

15   end

16   Find j_r such that $T{C}_{j_r}=\min \{T{C}_j/j\in J\prime \};$

17   $J\prime :=J\prime \setminus \left\{j_r\right\},E_r:=E_{j_r};$

18   If $T{C}_{j_r}<T{C}_{j^{*}}$ then $T{C}_{j^{*}}:=T{C}_{j_r}$ and $j^{*}:=j_r$

19   $E:=E_{j^{*}};$

20  until Either $J\prime =\left\{\right\}$ or E = E_r;

21 end

5. A semi-Lagrangian relaxation for the SPLPO

Now we use the technique of semi-Lagrangian relaxation, as this allows to close the duality gap. The equality constraints (1)(1) $$\begin{array}{c}\mathrm{Min}.\sum _{i\in I}\sum _{j\in J}{c}_{ij}{x}_{ij}+\sum _{j\in J}{f}_j{y}_j\\ \mathrm{s}.\mathrm{t}.\sum _{j\in J}{x}_{ij}=1,\hspace{1em}i\in I,\end{array}$$(1) have been split into two inequalities $\sum _{j\in J}{x}_{ij}\le 1$ and $\sum _{j\in J}{x}_{ij}\ge 1$ to obtain the following model: $$\begin{array}{c}v(\mathit{SLR}(\gamma )):==\underset{(x,y)}{\min }\sum _i\sum _j{c}_{ij}{x}_{ij}+\sum _j{f}_j{y}_j+\sum _i{\gamma }_i\left(1-\sum _j{x}_{ij}\right)\\ =\underset{(x,y)}{\min }\sum _i\sum _j\left(c_{ij}-{\gamma }_i\right)x_{ij}+\sum _j{f}_j{y}_j+\sum _i{\gamma }_i\\ =\underset{(x,y)}{\min }\sum _j\left(\sum _i\left(c_{ij}-{\gamma }_i\right)x_{ij}+f_j{y}_j\right)+\sum _i{\gamma }_i\end{array}$$ $\text{subject}\mathrm{ }\text{to:}(2),(3),(4),(5),$ and $\sum _{j\in J}{x}_{ij}\le 1.$ Every component of the multipliers vector γ is nonnegative.

As mentioned in Section 3, the objective function is concave and non-decreasing in its domain. Therefore, its semi-Lagrangian dual problem $$v(\mathit{SLD}):=\underset{{}_{\gamma \ge 0}}{\max }v(\mathit{SLR}(\gamma ))$$ can be solved with an ascent method. Also, as pointed out earlier, there are values ${\gamma }_i^{*}$ for all component i of γ and a set $Q=\{\gamma /{\gamma }_i\in [{\gamma }_i^{*},+\infty )\}$ such that for any $\gamma \in Q$ the optimum of $\mathit{SLR}(\gamma )$ is met.

For the SPLP (no preferences), Beltrán et al. (2012) proved that there is a closed interval where the search of the multipliers could be done. We provide the next two results for the SPLPO.

Theorem 5.1.

Let $c{p}_i:={\max }_j\{c_{ij}+f_j\}$ and $cp:={(c{p}_1,\dots ,c{p}_m)}^t$ be the maximum of the costs for each costumer i associated to each facility j and the vector of these costs, respectively. If $\gamma \ge cp$ (component wise), then $\gamma \in Q.$

Proof.

As we know, the semi-Lagrangian relaxation closes the duality gap if $\sum _{j\in J}{x}_{ij}=1$ for all i. By hypothesis, $c{p}_i-{\gamma }_i\le 0$ for all $i\in I.$ If we choose $j\prime$ such that $c{p}_i=c_{ij\prime }+f_{j\prime },$ then $(c_{ij\prime }-{\gamma }_i)+f_{j\prime }\le 0,$ and this inequality is true for any j since $j\prime$ gives the maximum among all $c_{ij}+f_j.$ Therefore, it cannot happen that $\sum _{j\in J}{x}_{ij}=0$ at an optimal solution. □

The result of Theorem 5.1 for the SPLPO is weaker than that obtained by Beltrán et al. (2012) for the SPLP in the following sense. They choose cp_i equal to ${\min }_j\{c_{ij}+f_j\},$ but with this, there is no guarantee that the preference constraints hold. We can set $x_{ij\prime }=1$ and $y_{j\prime }=1$ for the particular $j\prime$ such that ${\min }_j\{c_{ij}+f_j\}=c_{ij\prime }+f_{j\prime },$ but in this way not all the possible combinations of x_ij and y_j that take into account all customers’ preferences are considered. Therefore, we are taking the risk of obtaining a nonoptimal solution.

Let $c_i^1\le \dots \le c_i^n$ be the sorted costs c_ij for each $i\in I.$ The superindex represents the index $j\in J$ corresponding to the ascending order. That is, $c_i^1$ is the smallest value, and so on.

Theorem 5.2.

If $\gamma <c^1$, then $\gamma \notin Q.$

Proof.

By hypothesis $c_i^1-{\gamma }_i>0.$ Then $c_i^j-{\gamma }_i>0$ for all $j\in J.$ In that case, $x_{ij}^{*}=0$ for all $j\in J$ for any optimal solution. Therefore, it cannot happen that $\sum _{j\in J}{x}_{ij}=1$ for an optimal solution and we conclude that $\gamma \notin Q.$ □

The interval of search for the components of γ is defined as $$B=\gamma /c^1<\gamma \le cp+ϵ,ϵ>0\in {\mathbb{R}}^m,$$ since, as mentioned in Beltrán et al. (2012) for the SPLP case, $B\cap Q\ne \varnothing .$

Using the sorted costs $c_i^1\le \dots ,\le c_i^n,$ each component γ_i of γ can be either in an interval of the form $I_j=(c_i^j,c_i^{j+1}]$ where $j\in \{1,\dots ,m-1\},$ or out of it. For the first case, there is an infinity number of values of γ_i that can belong to a single interval $(c_i^j,c_i^{j+1}].$ But each of them has the same effect on the optimal value of $\mathit{SLR}(\gamma )$ since only a change from $x_{ij}=1$ to $x_{i,j+1}=1$ implies a change on the solution. Hence, we just need a single γ_i representative for the interval. As γ goes to infinity all combined costs $(c_{ij}-{\gamma }_i)+f_j$ become negative, and hence the semi-Lagrangian relaxation problem can be as difficult to solve as the original SPLPO problem. Then it is always convenient to choose a value ${\gamma }_i\in I_j$ as small as possible, that is, at distance ϵ from the lower bound of the interval. Also, it is easy to check that $c_i^n$ is always less than cp_i if $f_j>0.$ These facts will be applied in the ascent method used later.

5.1. Dual ascent method for the semi-Lagrangian dual

A dual ascent method (that we have named DAM) modifies the current values of the multipliers in order to achieve a steady increase on $\mathit{SLR}(\gamma )$ (see Bilde & Krarup, 1977). In our case this is possible due to the aforementioned properties (see Algorithm 3).

Algorithm 3:

Dual Ascent Method (DAM)

1 begin

  input: Problem $\mathit{SLR}({\gamma }^0).$

  output: A lower bound for $\mathit{SLR}(\gamma ),$ multipliers γ.

2  For all $i\in I$ sort c_ij as $c_i^1\le \dots \le c_i^n;$

3  $ϵ>0;$

4  $\mathtt{iter}:=0;$

5  $\mathtt{MAXiter}\in \mathbb{N};$

6  for $i\in I$ do

7   $c{p}_i:={\max }_j\{c_{ij}+f_j\};$

8   if ${\gamma }_i^0<c_i^1$ then

9    ${\gamma }_i^0:=c_i^1+ϵ;$

10   else if ${\gamma }_i^0>c_i^n$ then

11    ${\gamma }_i^0:=c_i^n;$

12   else

13    Find a j_i such that ${\gamma }_i\in (c_i^{j_i},c_i^{j_i+1}];$

14    ${\gamma }_i^0:=c_i^{j_i}+ϵ;$

15   end

16  end

17  repeat

18   Solve $\mathit{SLR}({\gamma }^{\mathtt{iter}});$

19   Find a subgradient $s^{\mathtt{iter}}$ for $\mathit{SLR}({\gamma }^{\mathtt{iter}});$

20   if $s^{\mathtt{iter}}\ne 0$ then

21    for $i\in I$ do

22     if $s_i^{\mathit{iter}}=1$ then

23      $j_i:=j_i+1;$

24      ${\gamma }_i^{\mathtt{iter}+1}:=\min \{c_i^{j_i}+ϵ,c{p}_i\};$

25     end

26    end

27    $\mathtt{iter}:=\mathtt{iter}+1;$

28   end

29  until Either $s^{\mathtt{iter}}=0$ or $\mathtt{iter}=\mathtt{MAXiter};$

30 end

An important issue is the starting point for the DAM. We have tested three different starting vectors:

1. γ = 0.

2. $\gamma =\mu ,$ where μ is a vector of multipliers that are found using SGM for LD associated to the relaxed constraints (1)(1) $$\begin{array}{c}\mathrm{Min}.\sum _{i\in I}\sum _{j\in J}{c}_{ij}{x}_{ij}+\sum _{j\in J}{f}_j{y}_j\\ \mathrm{s}.\mathrm{t}.\sum _{j\in J}{x}_{ij}=1,\hspace{1em}i\in I,\end{array}$$(1) .

3. Since $LR(\lambda ,\mu )$ has the integrality property, another option is to use the dual values of the linear relaxation of these problems as a starting point for DAM.

Although solving the linear relaxation takes less time than solving SGM (less than 6 s in the largest instance tested and less than 1 s in the smallest one), we noticed that, under the same parameter settings, DAM was slower. We obtained better results with the second approach, that is, solving the Lagrangian problem on a first stage in order to obtain multipliers far enough from zero but not too large so that the problems to be solved on the next stage are more tractable.

On line 3 of Algorithm 3, we find the maximum of the costs cp_i for each costumer i associated to each facility j. On line 3, the procedure needs to solve a sequence of problems $\mathit{SLR}(\gamma ),$ but due to the preference constraints, they are not as easy to solve as in the Lagrangian relaxation case $LR(\mu ,\lambda )$ proposed before. However, if $c_{ij}-{\gamma }_i>0,$ some variables x_ij can be set to 0.

A subgradient for $\mathit{SLR}(\gamma )$ is computed on line 3 of Algorithm 3 as $$s=\left[1-\sum _{j\in J}{x}_{ij}\forall i\right]=\left[s_{\gamma }\right]\in {\mathbb{R}}^m,$$ where each component of s belongs to {0, 1} as $\sum _{j\in J}{x}_{ij}\le 1$ must be satisfied.

Finally, the multipliers are updated (increased) on line 3 by moving from the current interval I_j to the next one for each component i in γ.

6. Speeding up the search for the optimal solution

After a certain number of iterations of DAM that have been fixed beforehand, we apply a variable fixing heuristic VFH. This procedure takes the set of variables y_j that are equal to 1 for a given solution. It chooses a subset of them, and the variables in this subset are fixed to 1. Then the subproblem with these fixed variables is solved. This is a smaller problem and easier to solve. The method is described in Algorithm 4.

Algorithm 4:

Variable Fixing Heuristic (VFH)

1 begin

  input: A vector γ from a solution obtained with DAM.

  output: A feasible solution for SPLPO.

2  $Y_{\gamma }:=\{y_j/y_j=1,j\in J\};$

3  $\mathtt{ps}\in [0,1]$ (a percentage $100\times \text{ps}\%$ of elements from $Y_{\gamma }$);

4  sort $Y_{\gamma }$ such that $y_j\prec y_k$ if $\sum _i(c_{ij}+f_j)<\sum _i(c_{ik}+f_k);$

5  choose the first $\mathtt{ps}\times 100\%$ smallest elements from the sorted $Y_{\gamma }$ to create set $Y_{\gamma }^{\mathit{subset}};$

6  fix $y_j:=1$ for all $y_j\in Y_{\gamma }^{\mathit{subset}};$

7  solve the SPLPO problem with an exact method;

8 end

Line 4 of Algorithm 4 shows the criterion that we used to sort variables y_j. We also tried to sort it by writing $y_j\prec y_k$ if $\sum _i{p}_{ij}<\sum _i{p}_{ik},$ where values p_ij are the preferences as described earlier. We obtained similar results, but the option proposed in Algorithm 4 performed slightly better.

6.1. Accelerated dual ascent procedure (ADA)

The procedure we propose to solve the SPLPO is summarised in Algorithm 5, which is basically made up of three stages, each involving one algorithm. First, SGM (Algorithm 1, explained in Section 4). Then, DAM (Algorithm 3, explained in Section 5). Finally, VFH (Algorithm 4, explained in Section 6).

Algorithm 5:

Accelerated Dual Ascent Algorithm (ADA)

1 begin

  input: An SPLPO problem.

  output: A feasible solution for SPLPO.

2  ${\mu }_i^0:={\min }_j\{c_{ij}+f_j\}$ for all $i\in I,{\lambda }_{ij}:=0$ for all $i\in ,j\in J;$

3  run $\text{SGM}({\mu }^0,{\lambda }^0)$ during $\mathtt{sgm}\_\mathtt{iter}$ iterations and find the iteration $\mathtt{best}\_\mathtt{sgm}\_\mathtt{iter}$ where is the best value of $LR(\mu ,\lambda );$

4  ${\gamma }^0:={\mu }^{\mathtt{best}\_\mathtt{sgm}\_\mathtt{iter}};$

5  run DAM with multipliers vector ${\gamma }^0$ during $\mathtt{dam}\_\mathtt{iter}$ iterations;

6  run DAM during $\mathtt{vfh}\_\mathtt{iter}$ iterations. In each iteration run VFH;

7  keep the best solution in the search process as the solution of this heuristic.

8 end

7. Computational results

In this section, we present the computational results obtained after having applied the algorithms introduced in the previous sections. The experiments have been carried out on a desktop computer with Intel® Xeon® 3.40 GHz processor and 16 Gb of RAM under a Windows® 10 operating system. When an MIP problem needed to be solved we used FICO Xpress® version 8.0. Part of the instances were taken from Cánovas et al. (2007), which are based on Beasley’s OR-Library, see Beasley (1990). These are:

• 131, 132, 133, 134: m = 50 and n = 50.

• a75_50, b75_50, c75_50: m = 75 and n = 50.

• a100_75, b100_75, c100_75: m = 100 and n = 75.

Additional larger data sets were generated using the same algorithm proposed in (Cánovas et al., 2007, Section 4, p. 145, Appendix A p. 149):

• a125_100, b125_100, c125_100: m = 125 and n = 100.

• a150_100, b150_100, c150_100: m = 150 and n = 100.

The headers of the tables have the following meanings:

• Prob: Name of the problem.

• Opt: Optimal objective function value of the problem Prob.

• y: Number of opened facilities.

• bestUB: Best upper bound.

• $v(P_L):$ Linear relaxation value for a problem.

• v(SGM): Best value using the subgradient method for LD (Algorithm 1).

• ${\text{GAP}}_{\mathrm{o}}=\text{bestUB}-\text{Opt}:$ Absolute gap between bestUB and Opt of a problem.

• ${\text{GAP}}_{\mathrm{o}}\%=\frac{{\text{GAP}}_{\mathrm{o}}}{\text{Opt}}\times 100\%:$ Relative gap between bestUB and Opt of a problem.

• ${\text{GAP}}_{\text{LP}}:\frac{v(P_L)-v(\mathit{SGM})}{v(P_L)}\times 100\%.$ The relative gap between the linear relaxation and the best lower bound v(SGM).

• ${\text{GAP}}_{\mathrm{X}}=\frac{\text{Opt}-v(P_L)}{\text{Opt}}\times 100\%:$ Relative gap between Opt and the linear relaxation value $v(P_L)$ of a problem by using Xpress.

• iter (itbsol): Number of iterations (iteration where best solution was found).

• t (Tt): CPU time in seconds (total time).

• sol10%: First solution found that is within 10% or less from the optimal solution.

• LB: Lower bound.

• impt: Percentage of time improved measuring Xpress versus ADA.

• ${\text{GAP}}_{\mathrm{V}}=\frac{\text{bestUB}-\text{LB}}{\text{bestUB}}\times 100\%:$ Relative gap between bestUB and LB of a problem by using DAM with VHF.

First we show in Table 2 the results for UBH (Algorithm 2) for the first 40 instances. It finds the optimal solution for 6 of the first 16 problems and in the rest of the instances it is not too far from the optimal with an average gap of 4.65%. All runs performed took less than 1 s whereas Xpress needed between 3 and 12 s to solve these instances to optimality. The value obtained by UBH was used as an upper bound LR_ub for SGM (Algorithm 1).

Table 2. Computational results for heuristic algorithm UBH.

Prob	Opt	y	UBH
			bestUB	GAP${}_{\mathrm{o}}$	GAP${}_{\mathrm{o}}$%	y
131_1	1001440	6	1001440	0	0.00	6
131_2	982517	9	982517	0	0.00	9
131_3	1039853	10	1139012	99159	9.54	9
131_4	1028447	9	1049791	21344	2.08	7
132_1	1122750	8	1239961	117211	10.44	7
132_2	1157722	9	1199057	41335	3.57	9
132_3	1146301	6	1146301	0	0.00	6
132_4	1036779	5	1036779	0	0.00	5
133_1	1103272	7	1103272	0	0.00	7
133_2	1035443	5	1100713	65270	6.30	7
133_3	1171331	6	1208198	36867	3.15	4
133_4	1083636	9	1090582	6946	0.64	9
134_1	1179639	4	1179639	0	0.00	4
134_2	1121633	7	1205809	84176	7.50	5
134_3	1171409	6	1173693	2284	0.19	7
134_4	1210465	3	1248143	37678	3.11	1
a75_50_1	1661269	7	1787955	126686	7.63	1
a75_50_2	1632907	6	1779576	146669	8.98	4
a75_50_3	1632213	7	1738404	106191	6.51	3
a75_50_4	1585028	5	1709978	124950	7.88	5
b75_50_1	1252804	8	1343201	90397	7.22	10
b75_50_2	1337446	9	1403629	66183	4.95	5
b75_50_3	1249750	9	1368788	119038	9.52	8
b75_50_4	1217508	9	1348203	130695	10.73	10
c75_50_1	1310193	11	1375386	65193	4.98	11
c75_50_2	1244255	10	1316595	72340	5.81	11
c75_50_3	1201706	12	1386817	185111	15.40	13
c75_50_4	1334782	11	1440836	106054	7.95	5
a100_75_1	2286397	4	2459349	172952	7.56	4
a100_75_2	2463187	3	2476632	13445	0.55	1
a100_75_3	2415836	3	2467003	51167	2.12	3
a100_75_4	2380150	4	2476632	96482	4.05	1
b100_75_1	1950231	8	2061201	110970	5.69	7
b100_75_2	2023097	8	2180286	157189	7.77	8
b100_75_3	2062595	8	2133724	71129	3.45	8
b100_75_4	1865323	9	1994265	128942	6.91	7
c100_75_1	1843620	6	2090221	246601	13.38	7
c100_75_2	1808867	11	2005071	196204	10.85	13
c100_75_3	1820587	8	2019651	199064	10.93	7
c100_75_4	1839007	9	2046525	207518	11.28	8

Table 3 shows the performance of the subgradient method SGM (Algorithm 1) over the first 40 data sets. The algorithm was stopped when the parameter β was equal to 0 with a starting value of 2, and it was decreasing linearly at a rate of 0.005 if the solution would not change for 30 iterations. We also show the results for the first solution within a 10% from $v(P_L).$ The decrease in the average number of iterations and execution time for this case is considerable, 65 and 69%, respectively. This shows that SGM can be stopped early to obtain good objective function values.

Table 3. Computational results for subgradient method.

Prob	$v(P_L)$	Until β = 0					First solution < 10%
		v(SGM)	itbsol	iter	t	GAP${}_{\text{LP}}\%$	sol10%	iter	t	GAP${}_{\text{LP}}\%$
131_1	925492	881876	1474	1500	19.42	4.71	833518	613	8.08	9.90
131_2	925195	872136	1450	1500	19.35	5.73	833747	695	9.46	9.90
131_3	955447	929854	1488	1500	20.02	2.68	860178	636	8.09	10.00
131_4	933025	903486	1462	1500	19.33	3.17	841581	580	7.11	9.80
132_1	1007417	981624	1500	1500	19.18	2.56	907222	479	6.21	9.90
132_2	990513	957530	1389	1425	18.01	3.33	897874	445	6.01	9.40
132_3	1009054	974902	1500	1500	19.09	3.38	911244	658	8.17	9.70
132_4	966305	910344	1388	1434	18.38	5.79	869850	604	8.34	10.00
133_1	998199	958109	1395	1498	19.24	4.02	898706	593	8.28	10.00
133_2	971719	947626	1443	1500	19.23	2.48	882263	423	5.19	9.50
133_3	1023593	982907	1473	1500	19.20	3.97	927636	565	7.36	9.40
133_4	1001253	948311	1489	1500	19.00	5.29	902108	634	8.45	9.90
134_1	1226933	1013753	1192	1226	16.28	2.21	933610	375	5.31	9.90
134_2	1041770	998540	1165	1239	16.16	4.15	942697	530	7.02	9.50
134_3	1023070	994435	1155	1286	16.31	2.80	921339	415	5.00	9.90
134_4	1050134	1003356	980	1026	13.44	4.45	945217	439	6.12	10.00
a75_50_1	1201542	1169011	732	850	17.02	2.71	1096795	106	2.32	9.40
a75_50_2	1188359	1158500	829	896	18.08	2.51	1076464	107	2.15	9.40
a75_50_3	1189183	1167203	936	999	20.41	1.85	1073799	157	3.08	9.70
a75_50_4	1200068	1175405	854	937	19.16	2.06	1094899	106	2.11	8.80
b75_50_1	900617	879687	1118	1182	24.17	2.32	812083	334	7.28	9.80
b75_50_2	922048	896384	1111	1211	25.22	2.78	830962	381	8.24	9.90
b75_50_3	897073	882132	1208	1271	25.27	1.67	808366	326	7.33	9.90
b75_50_4	885392	856157	784	925	19.35	3.30	801482	281	6.16	9.50
c75_50_1	892837	879074	845	876	18.19	1.54	804437	304	6.40	9.90
c75_50_2	882933	870772	1120	1175	23.23	1.38	797745	317	6.27	9.60
c75_50_3	859591	852325	1013	1050	22.21	0.85	786773	347	7.13	8.50
c75_50_4	928064	891598	924	1084	22.37	3.93	837696	415	8.32	9.70
a100_75_1	1800032	1725539	635	721	41.08	4.14	1626622	99	6.44	9.60
a100_75_2	1793377	1723484	679	728	42.48	3.90	1618944	99	6.15	9.70
a100_75_3	1788395	1724444	664	745	42.39	3.58	1632328	103	6.13	8.70
a100_75_4	1795298	1702235	482	678	38.33	5.18	1620377	100	6.17	9.70
b100_75_1	1330138	1289189	1010	1116	63.01	3.08	1198182	187	11.31	9.90
b100_75_2	1353824	1284338	652	691	39.18	5.13	1225149	373	21.44	9.50
b100_75_3	1351603	1306286	987	1044	59.23	3.35	1217952	270	15.23	9.90
b100_75_4	1332525	1258497	432	679	38.24	5.56	1201093	157	6.22	9.90
c100_75_1	1250876	1168284	510	710	40.39	6.60	1127235	303	17.03	9.90
c100_75_2	1239182	1154640	533	738	42.33	6.82	1119427	424	24.18	9.70
c100_75_3	1232462	1173522	670	804	46.28	4.78	1115951	270	15.11	9.50
c100_75_4	1243861	1194966	895	942	54.30	3.93	1121515	356	20.47	9.80

After the time reported, SGM did not obtain the value $v(P_L)$ for any of the problems. However, as mentioned before, method SGM is useful to find initial multipliers for the dual ascent DAM.

Finally, we can see the results for algorithm ADA in Table 4. The routine uses the following parameters:

Table 4. Computational results for ADA.

Prob	Xpress					UBH			SGM		DAM with VFH					ADA
	Opt	$v(P_L)$	GAP${}_{\mathrm{X}}\%$	y	t	bestUB	y	t	LB	t	bestUB	LB	GAP${}_{\mathrm{V}}\%$	y	t	Tt	impt %	GAP${}_{\mathrm{o}}\%$
a75_50_1	1661269	1201542	27.67	7	16.01	1787955	1	0.00	832797	4.25	1662877	1290126	22.42	5	21.32	25.57	−59.73	0.10
a75_50_2	1632907	1188359	27.22	6	21.10	1779576	4	0.00	828462	4.04	1648421	1261845	23.45	6	28.40	32.44	−53.77	0.95
a75_50_3	1632213	1189183	27.14	7	18.23	1738404	3	0.01	783044	3.11	1632866	1286922	21.19	7	24.38	27.49	−50.78	0.04
a75_50_4	1585028	1200068	24.29	5	18.28	1709978	5	0.00	863178	3.22	1585028	1250104	21.13	5	28.14	31.36	−71.70	0.00
b75_50_1	1252804	900617	28.11	8	11.05	1343201	10	0.00	355739	2.49	1252804	951766	24.03	8	28.28	30.76	−178.48	0.00
b75_50_2	1337446	922048	31.06	9	20.07	1403629	5	0.01	376669	3.22	1337446	954023	28.67	9	35.01	38.23	−90.46	0.00
b75_50_3	1249750	897073	28.22	9	17.24	1368788	8	0.00	348636	3.35	1255947	975724	22.31	11	39.27	42.62	−147.16	0.50
b75_50_4	1217508	885392	27.28	9	12.20	1348203	10	0.00	414596	3.27	1222561	978717	19.95	10	35.27	38.54	−215.77	0.42
c75_50_1	1310193	892837	31.85	11	17.47	1375386	11	0.00	550106	4.13	1310193	972827	25.75	11	38.38	42.51	−143.33	0.00
c75_50_2	1244255	882933	29.04	10	11.25	1316595	11	0.00	533586	3.21	1244255	964720	22.47	10	36.05	39.26	−248.98	0.00
c75_50_3	1201706	859591	28.47	12	8.45	1386817	13	0.00	549787	4.42	1201706	934664	22.22	12	33.18	37.60	−345.05	0.00
c75_50_4	1334782	928064	30.47	11	21.10	1440836	5	0.00	581078	4.07	1356714	974828	28.15	12	41.49	45.46	−115.48	1.64
a100_75_1	2286397	1800032	21.27	4	61.22	2459349	4	0.56	1649771	16.10	2321812	1950826	15.98	4	74.02	90.12	−47.20	1.55
a100_75_2	2463187	1793377	27.19	3	131.28	2476632	4	0.58	1630225	16.33	2476632	1982594	19.95	3	106.12	122.45	6.72	0.55
a100_75_3	2415836	1788395	25.97	3	148.41	2467003	3	0.75	1634629	16.46	2415836	1970116	18.45	3	99.46	115.92	21.89	0.00
a100_75_4	2380150	1795298	24.57	4	124.19	2476632	3	0.66	1628551	17.39	2386981	1938879	18.77	4	99.22	116.61	6.10	0.29
b100_75_1	1950231	1330138	31.80	8	624.08	2061201	7	0.71	1092987	17.41	1984720	1511894	23.82	10	125.41	142.81	77.12	1.77
b100_75_2	2023097	1353824	33.08	8	727.14	2180286	8	0.62	1150334	18.04	2058495	1517533	26.28	11	143.32	161.36	77.81	1.75
b100_75_3	2062595	1351603	34.47	8	841.07	2133724	8	0.55	1165120	18.31	2062595	1536082	25.53	8	268.38	285.69	65.91	0.00
b100_75_4	1865323	1332525	28.56	9	546.26	1994265	7	0.65	1143446	15.35	1886598	1512185	19.85	7	141.30	155.65	71.32	1.14
c100_75_1	1843620	1250876	32.15	6	430.40	2090221	7	0.70	1077180	20.28	1843620	1456543	21.00	6	215.21	235.49	45.8	0.00
c100_75_2	1808867	1239182	31.49	11	396.21	2005071	13	0.62	1053868	19.41	1815373	1415172	22.05	9	171.31	190.72	51.86	0.36
c100_75_3	1820587	1232462	32.30	8	423.16	2019651	7	0.60	1054526	16.15	1820587	1405713	22.79	8	194.30	210.44	50.27	0.00
c100_75_4	1839007	1243861	32.36	9	583.15	2046525	8	0.61	1064170	21.01	1839007	1417402	22.93	9	227.11	248.12	57.45	0.00
a125_100_1	3041451	2392412	21.34	2	257.17	3070535	1	0.62	2235753	47.44	3070535	2539124	17.31	2	151.48	198.93	22.65	0.96
a125_100_2	3040248	2393448	21.27	2	302.08	3070535	1	0.86	2232859	48.11	3070535	2629291	14.37	2	201.08	249.19	17.51	1.00
a125_100_3	3055260	2362216	22.68	3	325.07	3070535	1	0.60	2227182	47.31	3070535	2593865	15.52	3	209.36	256.67	21.04	0.50
a125_100_4	3056428	2381167	22.09	2	334.16	3070535	1	0.75	2244017	54.21	3070535	2529404	17.62	2	208.39	262.59	21.42	0.46
b125_100_1	2640798	1794710	32.04	7	5297.03	2850664	5	0.86	1601985	58.48	2640798	2082726	21.13	7	989.06	1047.53	80.22	0.00
b125_100_2	2550592	1790869	29.79	7	2086.19	2808259	9	0.95	1600582	56.24	2568522	2046667	20.32	7	848.03	904.27	56.65	0.70
b125_100_3	2604906	1792133	31.20	4	4059.45	2782609	5	0.91	1607076	58.09	2604906	1977647	24.08	4	435.16	493.95	87.85	0.00
b125_100_4	2580595	1792581	30.54	7	2686.31	2778713	4	0.89	1587315	60.43	2637611	1975516	25.10	7	358.06	418.49	84.42	2.21
c125_100_1	2491714	1663455	33.24	10	7805.07	2669965	9	0.79	1491505	61.43	2491714	1953129	21.62	10	1685.04	1746.47	77.62	0.00
c125_100_2	2468480	1674102	32.18	9	6296.32	2674261	8	0.91	1487416	58.23	2518055	1957245	22.27	10	519.22	577.46	90.83	2.01
c125_100_3	2559381	1672005	34.67	8	11381.22	2685807	7	0.90	1496552	59.10	2559381	1957575	23.51	8	1177.38	1236.48	89.14	0.00
c125_100_4	2538550	1691148	33.38	8	6407.39	2683676	9	0.88	1503317	61.29	2561098	2016969	21.25	8	726.04	786.33	87.71	0.89
a150_100_1	3768087	2906284	22.87	1	371.42	3768087	1	1.68	2691219	57.44	3768087	3200250	15.07	3	338.24	395.67	−6.53	0.00
a150_100_2	3768087	2891681	23.26	1	457.32	3768087	1	0.85	2695660	56.25	3768087	3212545	14.74	2	329.48	384.74	15.65	0.00
a150_100_3	3741364	2882917	22.94	3	340.47	3768087	1	0.91	2691438	58.44	3741364	3185719	14.85	3	306.47	363.91	−7.18	0.00
a150_100_4	3768087	2911017	22.75	1	591.15	3768087	1	0.95	2697003	63.41	3768087	3241759	13.97	2	337.10	400.51	32.25	0.00
b150_100_1	3271859	2159537	34.00	4	15483.45	3487710	4	0.95	1916485	68.42	3271859	2651812	18.95	4	2357.34	2425.76	84.33	0.00
b150_100_2	3227987	2160753	33.06	5	9677.02	3457623	5	1.77	1935777	68.19	3227987	2546429	21.11	5	975.05	1043.24	89.22	0.00
b150_100_3	3150075	2148658	31.79	6	5986.11	3390435	8	0.99	1920876	69.24	3150075	2582861	18.01	6	1081.41	1150.65	80.78	0.00
b150_100_4	3342783	2190488	34.47	5	9596.04	3637438	8	0.96	1927555	69.01	3342783	2602414	22.15	5	2225.04	2294.05	76.09	0.00
c150_100_1	2979389	1988908	33.24	9	10008.29	3196861	5	0.87	1790894	71.12	2979389	2408302	19.17	9	1520.24	1591.37	84.10	0.00
c150_100_2	3109105	1985389	36.14	7	21734.35	3291990	6	1.98	1766696	72.39	3109105	2404613	22.66	7	949.00	1021.39	95.30	0.00
c150_100_3	2937767	1982388	32.52	8	10590.24	3079664	4	0.97	1770318	71.10	2954519	2379867	19.45	7	2369.41	2440.51	76.96	0.57
c150_100_4	3165327	1997014	36.91	10	65134.01	3189247	7	0.98	1775186	75.22	3176587	2424997	23.66	12	5842.23	5917.45	90.91	0.36

• sgm_iter: Maximum number of iterations for subgradient algorithm SGM,

• dam_iter: Maximum number of iterations for dual ascent algorithm DAM,

• vfh_iter: Maximum number of iterations for variable fixing heuristic VFH.

The values for the parameters were empirically set with a calibration process based on the results from computational experiments obtained with the given problems. We note that the number of iterations needed in each procedure in ADA is related to the size of the problem. The best results were obtained when sgm_iter took a value of at least the number of customers n and dam_iter no more than 10% of n. It must also be noted that only a very small number of VFH iterations are necessary. Different values for sgm_iter, dam_iter y vfh_iter were used for the four groups of data sets tested, see Table 5.

Table 5. Parameter settings of the ADA algorithm for the indicated method.

m	n	SGM				DAM		VFH
		sgm_iter	β	k	q	dam_iter	ε	vfh_iter	ps
75	50	50	2	30	0.005	3	0.01	2	0.25
100	75	100	2	30	0.005	7	0.01	2	0.25
125	100	170	2	30	0.005	10	0.01	2	0.25
150	100	170	2	30	0.005	12	0.01	2	0.25

The good performance of UBH is maintained for larger instances ($m=125,n=100$) and ($m=150,n=100$), where the gap with respect to the optimum is less than 10%. UBH reached the optimum in five of these problems a125_100_3, a125_100_4, a150_100_1, a150_100_2, and a150_100_4. Except for the first group (m = 75, n = 50) and a few instances in the rest, the times obtained with ADA when it finds the optimal solution are considerably lower than those obtained with Xpress. It can be observed particularly in the last group. For example, instance c_150_100_2 had a reduction of 95% on the total time. A summary of the improvement on the CPU times when ADA reached the optimal solution is shown in Table 6.

Table 6. Average improved times Xpress vs ADA when the optimal solution is met.

m	n	Average impt %
75	50	−179.75
100	75	48.16
125	100	83.71
150	100	54.40

From the columns associated with DAM and VFH, we observe the impact that the acceleration heuristic has on the process of searching for the optimum. Before applying VFH, DAM provides a lower bound (LB), which is approximately 21% away from the optimum obtained by Xpress. ADA is always close to the optimal values, in fact, the average GAP is only 0.43%. Only in the smallest instances ADA does not work so well, but those times are negligible anyway.

In Table 7, we show a rough comparison between the times obtained by ADA and those reported in Vasilyev et al. (2013) for the formulation and method proposed by the authors and by Cánovas et al. (2007). Even though we are aware that the comparison is not fair due to the difference in hardware used and the fact that ADA is a heuristic procedure, we have included this summary so that the reader has at least an idea of how much time can be reduced using our approach considering the fact that ADA either finds the optimum (rows in boldface) or is very close to it in the reported cases.

Table 7. Times comparison.

Prob	Opt	Time (s)
		XPR	C&B	ADA	$\mathrm{GA}{\mathrm{P}}_{\mathrm{o}}\%$
a75_50_1	1661269	166	125	25	0.10
a75_50_2	1632907	147	100	32	0.95
a75_50_3	1632213	163	111	27	0.04
a75_50_4	1585028	143	79	31	0.00
b75_50_1	1252804	8	55	31	0.00
b75_50_2	1337446	158	134	38	0.00
b75_50_3	1249750	144	77	42	0.50
b75_50_4	1217508	80	39	39	0.42
c75_50_1	1310193	11	74	43	0.00
c75_50_2	1244255	68	47	39	0.00
c75_50_3	1201706	42	24	38	0.00
c75_50_4	1334782	144	78	45	1.64
Average		106	78	36	0.30
a100_75_1	2286397	1040	319	90	1.55
a100_75_2	2463187	1678	534	122	0.55
a100_75_3	2415836	1508	1037	116	0.00
a100_75_4	2380150	1210	822	116	0.29
b100_75_1	1950231	5205	3103	143	1.77
b100_75_2	2023097	7798	4120	161	1.75
b100_75_3	2062595	7947	5634	286	0.00
b100_75_4	1865323	3400	1492	156	1.14
c100_75_1	1843620	4713	2575	235	0.00
c100_75_2	1808867	4258	2569	191	0.36
c100_75_3	1820587	2962	1676	210	0.00
c100_75_4	1839007	4452	2732	248	0.00
Average		3847	2218	173	0.62

The headers for the time columns are:

• XPR: Xpress with settings and formulation from Cánovas et al. (2007).

• C&B: Cut-and-Branch method with the setting from Vasilyev et al. (2013).

• ADA: The ADA algorithm with the setting given in Table 5.

The hardware used for the experiments in columns XPR and C&B was a PC with Intel Core 2 Duo 1.8 GHz processor and 1 Gb of RAM. We have highlighted those cases where ADA reached the optimum.

8. Conclusions

In this article, we have showed for the first time in the literature how semi-Lagrangian relaxation can be a promising technique to solve the SPLPO. We have proposed an algorithm (ADA) that exploits these ideas and that performed very well on large instances, finding the optimal solution in most of the tested cases.

In the heuristic method ADA the assignment and VUBs constraints in the SPLPO were relaxed to formulate a Lagrangian problem. We solved its dual with a subgradient method and used the resulting multiplier vector as the starting point of an ascent algorithm for the dual of a semi-Lagrangian formulation. A better starting point for the dual ascent algorithm would be the multipliers obtained in the subgradient algorithm by relaxing only the assignment constraints, the same family of constraints relaxed in our proposed SPLPO semi-Lagrangian problem. However, the sequence of linear subproblems that need to be solved in the subgradient method are not as easy to solve as those when our proposed relaxation is used. We used the variable fixing heuristic VFH to speed up the search for the optimum since the subproblems in ADA become more and more difficult as the number of iterations grows.

In summary, this article reports the following findings:

1. A theoretical analysis of properties of Lagrangian and semi-Lagrangian formulations for the SPLPO is done.

2. A heuristic UBH to obtain fast reasonably good solutions for the SPLPO.

3. We combine the results obtained in the theoretical analysis in a heuristic procedure that finds good solutions, often an optimal solution.

4. The computational experiments carried out show that ADA performs very well on large instances.

A limitation of our study, which could be improved in future works, is that, since the last stage of ADA is a heuristic, finding an optimal solution is not guaranteed, and, as in any heuristic procedure, the parameter setting of the algorithms is an issue that needs special attention. In our case, the relationship between the size of a particular instance and the proposed values is clear.

Disclosure Statement

No potential conflict of interest was reported by the author(s).