Persisting Subproblem Infeasibility Benders Decomposition

Question

I'm currently working on an optimization problem where I'm using Benders decomposition to solve a complex problem involving the installation of charging stations. The master problem determines the placement of charging stations, while the subproblem involves optimal power flow. I've encountered an issue with subproblem infeasibility despite using feasibility cuts.

Here's an overview of my approach and the formulation of my feasibility cut:

• I've fixed decision variables $x_{i,t}$ in the master problem, which determine the number of charging stations to be installed at location bus $i$ of type $t$.
• The primal subproblem is presented in its canonical form, and I've provided the associated dual variables next to each constraint. I transformed each equality constrained to two inequality constraints.
• I've formulated a feasibility cut. This cut is derived from the extreme ray values of the subproblem's dual variables.

This is the linear subproblem where $P_k^g,Q_k^g,P_{ij},Q_{ij},V_i$ and $\delta_i$ are continuous decision variables: $$min \sum_{k\in G}\left( C_{1,k} P_{k}^{g} +C_{0,k}\right)$$ S.t. $$P_{k}^{g} \geq P_{k}^{gl} \quad \forall k\in G \qquad 𝛾_{k} $$ $$-P_{k}^{g} \geq -P_{k}^{gu} \quad \forall k\in G \qquad 𝜓_{k}$$ $$Q_{k}^{g} \geq Q_{k}^{gl} \quad \forall k\in G \qquad \epsilon_{k}$$ $$-Q_{k}^{g} \geq -Q_{k}^{gu} \quad \forall k\in G \qquad \zeta_{k}$$ $$V_{i} \geq v_{i}^{l} \quad \forall i\in N \qquad \eta_{i}$$ $$-V_{i} \geq - v_{i}^{u} \quad \forall i\in N \qquad \tau_{i}$$ $$-\delta_{i} +\delta_{j} \geq -\theta_{ij}^{\Delta u} \quad \forall ( i,j) \in E \qquad \xi_{ij}$$ $$\delta_{i} -\delta_{j} \geq \theta_{ij}^{\Delta l} \quad \forall ( i,j) \in E \qquad \omicron_{ij}$$ $$\sum_{k\in G_{i}} P_{k}^{g} -\sum_{( i,j) \in E_i} P_{ij} \geq \sum_{k\in L_{i}} P_{k}^{d} +\sum_{t\in T} P_{t}^{cs}\overline{x_{i,t}} \quad \forall i\in N \qquad \iota_{i}$$ $$-\sum_{k\in G_{i}} P_{k}^{g} +\sum_{( i,j) \in E_i} P_{ij} \geq -\sum_{k\in L_{i}} P_{k}^{d} -\sum_{t\in T} P_{t}^{cs}\overline{x_{i,t}} \quad \forall i\in N \qquad \kappa_{i}$$ $$\sum_{k\in G_{i}} Q_{k}^{g} -\sum_{( i,j) \in E_{i}} Q_{ij} \geq \sum_{k\in L_{i}} Q_{k}^{d} \quad \forall i\in N \qquad \lambda_{i}$$ $$-\sum_{k\in G_{i}} Q_{k}^{g} +\sum_{( i,j) \in E_{i}} Q_{ij} \geq -\sum_{k\in L_{i}} Q_{k}^{d} \quad \forall i\in N \qquad \mu_{i}$$ $$P_{ij} -( V_{i} -V_{j}) G_{ij} +B_{ij} ( \delta_{i} -\delta_{j}) \geq 0 \quad \forall ( i,j) \in E \qquad \rho_{ij}$$ $$-P_{ij} +( V_{i} -V_{j}) G_{ij} -B_{ij} ( \delta_{i} -\delta_{j}) \geq 0 \quad \forall ( i,j) \in E \qquad \sigma_{ij}$$ $$Q_{ij} -( V_{j} -V_{i}) B_{ij} +G_{ij} ( \delta_{i} -\delta_{j}) \geq 0 \qquad \forall ( i,j) \in E \qquad \upsilon_{ij}$$ $$-Q_{ij} +( V_{j} -V_{i}) B_{ij} -G_{ij} ( \delta_{i} -\delta_{j}) \geq 0 \quad \forall ( i,j) \in E \qquad \omega_{ij}$$

Multiplying the RHS of each constraint with its dual variable, should give the objective function of the dual of the subproblem:

$$ max \sum_{k\in G} ( P_{k}^{gl} \gamma_{k} -P_{k}^{gu} \psi_{k} +Q_{k}^{gl} \epsilon_{k} -Q_{k}^{gu} \zeta_{k}) +\sum_{i\in N}( v_{i}^{l} \eta_{i} - v_{i}^{u} \tau_{i}) +\sum_{i\in N} ((\sum_{k\in L_{i}} P_{k}^{d} +\sum_{t\in T} P_{t}^{cs}\overline{x_{i,t}}) \iota_{i} -(\sum_{k\in L_{i}} P_{k}^{d} +\sum_{t\in T} P_{t}^{cs}\overline{x_{i,t}}) \kappa_{i}) +\sum_{i\in N}(\sum_{k\in L_{i}} Q_{k}^{d} \lambda_{i} -\sum_{k\in L_{i}} Q_{k}^{d} \mu_{i}) +\sum_{( i,j) \in E}( \theta_{ij}^{\Delta l} \omicron_{ij} -\theta_{ij}^{\Delta u} \xi_{ij}) $$

I've implemented this approach using JuMP with Julia and Gurobi as the solver. I retrieve the dual extreme ray values using JuMP.shadow_price(ConstraintRef). As for now, I follow the classical approach without lazy constraints to test the implementation. My understanding is that I should directly incorporate the extreme ray values obtained from the unbounded dual of the subproblem into the objective function of the dual of the subproblem. This results in the following feasibility cut:

$$\sum_{k\in G}\left( P_{k}^{gl} \overline{\gamma_{k}} -P_{k}^{gu} \overline{\psi_{k}} +Q_{k}^{gl} \overline{\epsilon_{k}} -Q_{k}^{gu} \overline{\zeta_{k}}\right) +\sum_{i\in N}\left( v_{i}^{l} \overline{ \eta_{i}} - v_{i}^{u} \overline{\tau_{i}}\right) +\sum_{i\in N}\left(\left(\sum_{k\in L_{i}} P_{k}^{d} +\sum_{t\in T} P_{t}^{cs} x_{i,t} \right)\overline{ \iota_{i}} -\left(\sum_{k\in L_{i}} P_{k}^{d} +\sum_{t\in T} P_{t}^{cs} x_{i,t}\right) \overline{\kappa_{i}}\right) +\sum_{i\in N}\left(\sum_{k\in L_{i}} Q_{k}^{d} \overline{\lambda_{i}} -\sum_{k\in L_{i}} Q_{k}^{d} \overline{\mu_{i}}\right) +\sum_{( i,j) \in E}\left( \theta_{ij}^{\Delta l} \overline{\omicron_{ij}} -\theta_{ij}^{\Delta u} \overline{\xi_{ij}}\right) \leq 0$$

However, the master problem still propagates the same infeasible solution to the subproblem in subsequent iterations.

I'm looking for insights into why this feasibility cut is not effectively cutting away the infeasible solution from the master problem. Are there any potential mistakes in the formulation or any suggestions on how to improve the feasibility cut?

Answer 1

After some trial and error, I fixed the following two problems:

Problem 1: the formulation provided here is indeed correct as mentioned above, but I forgot to add bounds on the variables $P_{ij}$ and $Q_{ij}$. This results in additional dual variables that should be added to the feasibility cut.

Problem 2: the function JuMP.shadow_price() returns the wrong values (see Benders Subproblem Infeasibility). Use JuMP.dual() to get the correct dual extreme ray values.

Answer 2

Disclaimer #1: I did not check your formulation of the dual problem and the Benders cut carefully. It looks plausible, but the density of the notation made my eyes water.

Disclaimer #2: I'm not a JuMP user. You'll need to check what I'm about to say with someone who is.

I suspect you may be calling the wrong function. When a solver detects an unbounded LP, it typically is sitting at an extreme point of the feasible region looking at an edge in a desirable direction that extends forever. The dual solution (which is what I suspect the shadow_price() function returns, based on its name) is distinct form the dual ray (the direction of that edge). If I'm right, you need to find a function that returns the ray.