# Persisting Subproblem Infeasibility Benders Decomposition

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?

• When you say $P_g$ is a decision variable, do you mean $P^g_k?$ Also, are the $Q^g$ and $Q_{ij}$ symbols variables?
– prubin
Aug 8, 2023 at 19:11
• We need to know which elements of the primal subproblem are master problem variables.
– prubin
Aug 8, 2023 at 19:13
• Yes, sorry for the incomplete information. The decision variables for the subproblem are $P_k^g, Q_k^g, P_{ij}, Q_{ij}, V_i$ and $\delta_i$. $x_{i,t}$ is the only variable from the master problem that is present in the subproblem. The other symbols in the subproblem are data. Aug 8, 2023 at 19:21

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.

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.

• Thank you for the comments, I was indeed calling the wrong function (see answer). Aug 11, 2023 at 13:11