I have a two-stage multi-scenario MILP in JuMP with an inventory/order-flow structure.
The first-stage variables are shared across all scenarios. Each scenario then has its own continuous flow variables and binary timing variables. The objective and epigraph constraints depend on the shared first-stage variables. For fixed first-stage variables, the scenario-specific part is essentially a feasibility recourse problem.
A runnable extensive-form model builder with example data is here:
By default, the script builds a small two-scenario instance using scenarios 0 and 1. The data directory also contains additional scenarios, but the larger extensive form is the scale at which solving directly becomes difficult.
It can be run with:
julia OptBMExtensive.jl
The script builds the extensive-form model and reports its size; it does not attempt to solve the large instance directly. I have not included my attempted Benders implementation because I am mainly trying to understand the correct formulation, decomposition, and cut structure, rather than debug a particular implementation.
Model structure
The sets are i \in I for initial items, t,k \in T for periods, and s \in S for scenarios.
The shared first-stage variables are x_i and \theta, with 0 \le x_i \le 1 for all i \in I, and \theta \ge 0.
The objective is \min \; \theta.
The scenario-dependent epigraph constraints are \theta \ge \sum_{i \in I} v^0_{s,i,1} x_i for all s \in S.
For each scenario, the remaining variables are scenario-specific. They include order quantities y_{s,t} \ge 0, usage variables 0 \le u^0_{s,i,t} \le 1 and u^1_{s,k,t} \ge 0, stock variables q^0_{s,i,t} \ge 0 and q^1_{s,k,t} \ge 0, and binary timing variables \delta_{s,t} \in \{0,1\} and \gamma_{s,t} \in \{0,1\}.
Here y is an order quantity, u is usage, q is remaining stock, \delta indicates whether a period is an order period, and \gamma indicates whether initial items may be used in that period.
The initial-item flow constraints are q^0_{s,i,1} + u^0_{s,i,1} = x_i for all s \in S and i \in I, and q^0_{s,i,t} + u^0_{s,i,t} = q^0_{s,i,t-1} for all s \in S, i \in I, and t > 1.
The ordered-item flow constraints are q^1_{s,k,k} = y_{s,k} for all s \in S and k \in T, and q^1_{s,k,t} + u^1_{s,k,t} = q^1_{s,k,t-1} for all s \in S, k \in T, and t > k.
For compactness, define the following period-level expressions:
A^0_{s,t} = \sum_{i \in I} r^0_{s,i,t} q^0_{s,i,t},
A^1_{s,t} = \sum_{k<t} r^1_{s,k,t} q^1_{s,k,t},
B^0_{s,t} = \sum_{i \in I} v^0_{s,i,t} u^0_{s,i,t},
B^1_{s,t} = \sum_{k<t} v^1_{s,k,t} u^1_{s,k,t}.
The per-period balance constraint is A^0_{s,t} + A^1_{s,t} + B^0_{s,t} + B^1_{s,t} - p_{s,t} y_{s,t} = d_{s,t} for all s \in S and t \in T.
The timing logic is represented with big-M constraints:
y_{s,t} \le M^y_{s,t}\delta_{s,t} for all s \in S and t \in T.
\sum_{i \in I} u^0_{s,i,t} + \sum_{k<t} u^1_{s,k,t} \le M^u_{s,t}(1-\delta_{s,t}) for all s \in S and t \in T.
\sum_{k \le t} q^1_{s,k,t} \le M^q_{s,t}(1-\gamma_{s,t}) for all s \in S and t \in T.
\sum_{i \in I} u^0_{s,i,t} \le M^0_{s,t}\gamma_{s,t} for all s \in S and t \in T.
The interpretation of these binary constraints is:
-
The binary variable \delta_{s,t} enforces order/use exclusivity in each scenario and period. If \delta_{s,t}=1, then y_{s,t} may be positive but all usage in that period is forced to zero. If \delta_{s,t}=0, then y_{s,t}=0 and usage may occur.
-
The binary variable \gamma_{s,t} enforces a priority rule between ordered stock and initial-item usage. If \gamma_{s,t}=1, then initial items may be used, but remaining ordered stock at the end of the period is forced to zero. If \gamma_{s,t}=0, then initial-item usage is forced to zero.
So the intended logic is: a period cannot both order and use inventory, and initial items can only be used once ordered stock has been exhausted by the end of that period.
Therefore, for fixed first-stage variables x, each scenario subproblem is a MILP feasibility problem: Q_s(x)=0 if scenario s is feasible for x, and Q_s(x)=+\infty otherwise.
Approximate JuMP structure
The JuMP structure is roughly:
model = Model()
@variable(model, 0 <= initial_qty[1:n_items] <= 1)
@variable(model, theta >= 0)
for s in scenarios
@variable(model, order_qty[1:n_periods] >= 0)
@variable(model, 0 <= use_initial[1:n_items, 1:n_periods] <= 1)
@variable(model, use_order[1:n_periods, 1:n_periods] >= 0)
@variable(model, initial_stock[1:n_items, 1:n_periods] >= 0)
@variable(model, order_stock[k = 1:n_periods, t = k:n_periods] >= 0)
@variable(model, order_period[1:n_periods], Bin)
@variable(model, initial_use_period[1:n_periods], Bin)
# Coupling to the first-stage variables
@constraint(model, [i = 1:n_items],
initial_stock[i, 1] + use_initial[i, 1] == initial_qty[i]
)
@constraint(model, [i = 1:n_items, t = 2:n_periods],
initial_stock[i, t] + use_initial[i, t] == initial_stock[i, t - 1]
)
# Ordered-item flow
@constraint(model, [k = 1:n_periods],
order_stock[k, k] == order_qty[k]
)
@constraint(model, [k = 1:n_periods, t = k+1:n_periods],
order_stock[k, t] + use_order[k, t] == order_stock[k, t - 1]
)
# Scenario balance constraints and big-M timing constraints follow here.
@constraint(model,
theta >= sum(initial_value[s, i, 1] * initial_qty[i] for i in 1:n_items)
)
end
@objective(model, Min, theta)
Formulation/decomposition question
I am unsure what the right decomposition should be for this structure.
One possible decomposition is:
- master problem: shared variables x and \theta;
- scenario subproblem: all scenario-specific variables for fixed x.
However, after fixing x, each scenario subproblem still contains the binary variables \delta and \gamma. Therefore, the scenario subproblem is a MILP feasibility problem, not a continuous LP feasibility problem. Classical LP-dual Benders feasibility cuts do not seem directly applicable.
Another possibility is to move the timing binaries \delta and \gamma into the master. Then the scenario subproblem becomes continuous, so Farkas/dual feasibility cuts should in principle be available. However, this makes the master much larger and more combinatorial, and on larger instances this has not worked well.
So my main question is:
What master/subproblem split would you recommend for this model, and what valid feasibility cuts would that split imply?
More specifically:
-
Should the master contain only the shared first-stage variables x and \theta, leaving the binary timing variables in the scenario subproblems? If so, what kind of logic-based Benders or feasibility cuts would be valid for this structure?
-
Or should the master include the timing binaries \delta and \gamma so that each scenario subproblem is continuous? If so, would the right cuts be LP Farkas feasibility cuts, and what would those cuts look like structurally?
-
Are there stronger cuts or reformulations that exploit the flow/timing structure, rather than just excluding one infeasible binary pattern with a no-good cut?
-
Once the formulation is chosen, would you recommend implementing it in JuMP as an outer-loop decomposition, or as branch-and-Benders using lazy constraints?
I would especially appreciate advice on what should be in the master, what should be in the subproblem, and what the valid feasibility cuts would look like.
, I’m now resting. I’ll see if I can share some insights after recovery.