Accelerated Benders Decomposition and Local Branching for Dynamic Maximum Covering Location Problems

We detail the development of the branch-and-Benders-cut method proposed in Cordeau et al. (2019). While the process is nearly identical in the dynamic case as the static case, the development process is a prerequisite for the acceleration techniques we discuss below. As in Cordeau et al. (2019), we define $J_{s}=\left\{j\in J:\sum_{t=1}^{T}\sum_{i\in I}a_{ij}^{t}=1\right\}$ as the set of users which are covered by only one facility.

Like in the classical Benders decomposition (Benders 1962), we project out the (continuous) $z$ variables of model (1). The value of the $z$ variables in the objective function is replaced with an auxiliary variable, $\theta$. Then, in iteration $v$, the main problem (MP) can be written as

where $\pi^{r}$ and $\sigma^{r}$ are the optimal dual vectors associated with the Benders subproblem in iteration $r$. Since the Benders subproblem is feasible $\forall x\in\Omega$, it is not necessary to include feasibility cuts.

We denote —now and throughout this paper— $\tilde{x}$ as the candidate solution in iteration $v$. Unless otherwise specified, it is assumed that this solution is integer feasible, i.e. that $\tilde{x}\in\Omega$ and $\tilde{x}_{i}^{t}\in\{0,1\},1\leq t\leq T,i\in I$. For each $1\leq t\leq T$ and $j\in J$, let $I_{j}^{t}(x)=\sum_{i\in I}a_{ij}^{t}x_{i}^{t}$ denote the coverage of user $j$ in period $t$ by solution $x$. When there is no risk of confusion, we use $\tilde{I}_{j}^{t}=I_{j}^{t}(\tilde{x})$.

The Benders Primal Subproblem (BPS) is then given by

Let $\pi_{j}^{t}$ and $\sigma_{j}^{t}$ denote the dual variables associated with Constraints (4b) and (4c). The Benders Dual Subproblem (BDS) is then

Similar to Cordeau et al. (2019), this subproblem can be easily solved by inspection: If $\tilde{I}_{j}^{t}<1$, then the optimal solutions are $\tilde{\pi}_{j}^{t}=d_{j}^{t},\tilde{\sigma}_{j}^{t}=0$. If $\tilde{I}_{j}^{t}>1$, then $\tilde{\pi}_{j}^{t}=0,\tilde{\sigma}_{j}^{t}=d_{j}^{t}$. If $\tilde{I}_{j}^{t}=1$, then any solution $\tilde{\pi},\tilde{\sigma}\geq 0$ such that $\tilde{\pi}_{j}^{t}+\tilde{\sigma}_{j}^{t}=d_{j}^{t}$ will be optimal. The optimality cut associated with these values is then

The resulting cut, reformulated as a function of the variables $x$, is presented in Proposition 4.1. We set $\pi^{v+1},\sigma^{v+1}=\tilde{\pi},\tilde{\sigma}$, and the left-hand side of the cut forms the term $\operatorname{Opt}_{\pi^{v+1},\sigma^{v+1}}$ in (3b).

We remark that the names of the cuts (B0), (B1), and (B2) correspond to those of the equivalent cuts in Cordeau et al. (2019).

We note the following observations about these optimality cuts:

• •

  As in Cordeau et al. (2019), the cuts (B0) are dominated by the cuts (B1).

• •

  By definition of the set $J_{s}$, the condition $\tilde{I}_{j}^{t}\leq 1$ will always be satisfied for $j\in J_{s}$.

• •

  The optimality cuts  (6) are valid for both integer and fractional candidate solutions $\tilde{x}$. As described in Cordeau et al. (2019), optimality cuts can be generated for fractional solutions to improve the convergence rate.

• •

  When evaluated at solution $\tilde{x}$, the value of the left-hand side of the optimality cut (6) (and, thus, the corresponding objective value) matches the value of $\tilde{x}$ in Equation (2). This derives directly from the definition of $\Gamma^{t}(\tilde{x}^{t})$. However, this evaluation is not correct (i.e. equal to $f(\tilde{x}$ in Equation (2)) in general if evaluated at a different solution.

  We present a simple example of this in Figure 1. We drop the superscript $t$ as there is only one time period. If we generate the optimality cuts for the candidate solution $\tilde{x}=(1,0,0)$, i.e. facility $i_{1}$ is opened whereas facilities $i_{2}$ and $i_{3}$ are closed, we have the following optimality cuts

  	(B0):	$\displaystyle 0x_{i_{1}}+8x_{i_{2}}+10x_{i_{3}}+20\geq\theta,$
  	(B1):	$\displaystyle 10x_{i_{1}}+8x_{i_{2}}+10x_{i_{3}}+10\geq\theta,$
  	(B2):	$\displaystyle 20x_{i_{1}}+16x_{i_{2}}+12x_{i_{3}}+0\geq\theta.$

  All three types of cuts correctly identify the objective value of solution $\tilde{x}$ as $20$, matching the value in Equation (2). However, if we then use these optimality cuts to calculate the objective value of the solution $\bar{x}=(0,0.75,0.75)$, we obtain the values $33.5,23.5,$ and $21$ (for B0, B1, and B2 cuts respectively), while Equation (2) gives $19.5$.

So far, we have presented a simple generalisation of the Benders procedure of Cordeau et al. (2019) to the dynamic MCLP. However, as we will see in Section 6.2, this method can exhibit slow convergence. We next explore improvements with respects to its convergence rate.

In this section, we demonstrate and discuss the application of known techniques to our particular dynamic MCLP. These approaches aim to strengthen the Benders optimality cuts and improve the dual bound for the decomposition. For the sake of simplicity, the acceleration techniques are presented separately in relation to the decomposition presented in Section 4.1. However, it is possible to incorporate multiple techniques simultaneously. For a discussion of acceleration techniques in the general context of Benders decomposition, we refer to Rahmaniani et al. (2017).

There exist classes of improvement methods which heavily depend on the structure of $\Omega$, but which do not affect the structure of the Benders subproblems. Since they involve $\Omega$, their use is application-specific. However, since they do not impact the Benders subproblems, they can be applied to the methods in this work without affecting their validity. An example of such a technique would be a primal heuristic providing warmstart solutions. This can be used to obtain a good lower bound on the optimal objective value, which otherwise can be difficult during the early iterations of the Benders decomposition method (Rahmaniani et al. 2017). Another example would be valid cuts, which can be added to the main problem to tighten its linear programming (LP) relaxation. Examples for both techniques can be found in Santoso et al. (2005), Codato and Fischetti (2006), Costa et al. (2012).

The remainder of the techniques discussed in this section modify the Benders subproblems in some way, whether by changing the set of users $J$ (which, in turn, changes the amount of computational work required for the generation of cuts), modifying the frequency or type of cuts, or by directly providing reformulations for the BPS (4) or the BDS (5). When an aspect of the technique depends on the structure of $\Omega$ (and thus may be application-specific), we include a discussion on this dependence.

We implement the local branching scheme within the framework of our branch-and-Benders-cut methods, as proposed in Rei et al. (2009). Since it is based on our branch-and-Benders-cut methods, any of the acceleration techniques presented in Section 4.2 can be applied to the local branching method as well. To that end, contrary to Section 4.2, we present this section using the multi-cut version of the Benders optimality cuts. Due to the central role of the time period in the proceeding developments, the multi-cut formulation allows for a more natural and interpretable explanation. Thus, at iteration $v$ of the Benders decomposition, we consider the following local branching main problem:

where $x^{s}$ for $1\leq s\leq u$ are feasible facility location decisions previously determined. The function $\operatorname{Dist}_{x^{s}}$ represents a distance metric, indicating that only solutions found at a distance ${\kappa}^{s}$ from solution $x^{s}$ may be considered. We note that Constraints (19b) are the Benders optimality cuts (7c) provided by Proposition 4.2.

Our goal, given an integer candidate solution $\tilde{x}$ and a threshold distance $\tilde{\kappa}$, is to find the optimal solution to our original problem (1) restricted to distance $\tilde{\kappa}$ around $\tilde{x}$. As described in Rei et al. (2009), we can then generate Benders optimality cuts (8) for this high-quality solution, thus improving the upper bound in the local branching main problem (19). Simultaneously, the solution may also improve upon the incumbent, thus also increasing the lower bound. We can then exclude the subdomain of distance $\tilde{\kappa}$ around $\tilde{x}$ from the local branching main problem (19) by adding a constraint of type (19c).

Formally, let $\hat{x}$ denote the optimal solution to the restricted subproblem centered around $\tilde{x}$, given by the solution to the following subproblem:

If $\hat{x}$ has a strictly higher objective than $\tilde{x}$, we generate Benders optimality cuts (8) for $\tilde{x}$, and we set $x^{s+1}=\tilde{x},\kappa^{s+1}=\tilde{\kappa}+1,s=s+1$ and repeat the restricted problem (20) centered around $\hat{x}$. If $\hat{x}=\tilde{x}$, we create a diversified subproblem which replaces the threshold distance $\tilde{\kappa}$ in Constraint (20b) with $\tilde{\kappa}^{\prime}>\tilde{\kappa}$, and adds the constraint $\operatorname{Dist}_{\tilde{x}}\left(x\right)\geq 1$. This guarantees that the resulting optimal solution $\hat{x}$ is different than $\tilde{x}$, and we create Benders optimality cuts (8) from $\hat{x}$.

By setting $x^{s+1}=\tilde{x},\kappa^{s+1}=\tilde{\kappa}+1$, we ensure that the feasible domain in the restricted subproblem (20) in iteration $v$ and the main problem in iteration $v+1$ are complementary. More specifically, the restricted subproblem only considers solutions with distance $\operatorname{Dist}_{\tilde{x}}\left(x\right)\leq\tilde{\kappa}$, whereas the main problem considers solutions with distance $\operatorname{Dist}_{\tilde{x}}\left(x\right)\geq\tilde{\kappa}+1$. Since the candidate solution $\tilde{x}$ is integer feasible, it is not necessary to consider solutions with distance $\operatorname{Dist}_{\tilde{x}}\left(x\right)\in(\tilde{\kappa},\tilde{\kappa}+1)$, as these correspond to fractional solutions.

We remark that the local branching main problem (19) and the restricted subproblem (20) can be defined via any distance metric and with any threshold distance $\tilde{\kappa}$. By increasing the size of the subdomain, we can remove a larger area from the feasible space of the local branching main problem (19), but at the cost of a restricted subproblem (20) that is harder to solve than that of a smaller size. We also note that any method can be used to find the optimal solution $\hat{x}$ to the restricted subproblem (20). In Fischetti and Lodi (2003), a general-purpose mixed-integer linear solver is used, while the branch-and-Benders-cut approach with the acceleration techniques proposed in Section 4.2 is a natural choice in our case. However, by carefully selecting our distance metric and threshold distance $\tilde{\kappa}$, we can derive an exact solution method for the restriced subproblem (20) which is more effective than the accelerated branch-and-Benders-cut method.

The distance metric typically used is the Hamming distance (see, e.g., Fischetti and Lodi 2003, Rei et al. 2009). This results in Constraint (20b) taking the form

However, we propose to use a modified distance metric in the dynamic case, which enforces that the Hamming distance in each time period must be within our threshold. This results in Constraint (20b) taking the form

The motivation behind this new distance metric derives from infrastructure contexts, where the facilities under consideration correspond to significant investments (e.g. warehouses, stores, etc.). In those contexts, facilities which are added in early time periods are likely to persist throughout the time horizon. In the case of the general Hamming distance, this incurs a repeated penalty in each time period, increasing the total distance. As an illustration, consider the example in Figure 1 and the two solutions $x^{t}=(1,0,0)$ and $\hat{x}^{t}=(1,1,0),1\leq t\leq T$, and a time horizon $T=4$. These solutions are at distance $4$ when considering the general Hamming distance, whereas they are only at distance $1$ using the new metric.

An important benefit of this new distance metric allows for an efficient solving method if we consider $\tilde{\kappa}=2$. To describe this method we first note that, for a candidate solution $\tilde{x}$ and time period $t$, there are only two sets of modifications which result in an integer solution at distance exactly $1$:

1. 1a)

   Add one facility which is currently not selected.

2. 1b)

   Remove one facility which is currently selected.

Likewise, at distance exactly $2$, there are only three sets of modifications:

1. 2a)

   Add two facilities which are currently not selected.

2. 2b)

   Remove one facility which is currently selected, and add one facility which is currently not selected.

3. 2c)

   Remove two facilities which are currently selected.

Since, by assumption $a_{ij}^{t}\in\left\{0,1\right\}$ and $d_{j}^{t}>0$, options 1b) and 2c) cannot lead to an increase in the objective function, and hence, can be disregarded.

As a consequence, there is a limited number of feasible solutions contained within the restricted subproblem (20) when considering $\tilde{\kappa}=2$. We can then use Proposition 5.2 to evaluate the quality of feasible solutions (as given in Equation (2)) based on Benders optimality cuts.