I am trying out Julia to model a stochastic network interdiction problem with the code below, where the objective of the first stage is to minimize the expected value of the second stage. However, I’ve encountered two (probably) related issues. 1. The interdiction decision variable, y, is never non-zero. 2. The objective function value is negative when it should only be positive as it is the sum of the flow into the last node.
using StochasticPrograms using GLPK # Network Interdiction Model M = 100 # Big M B = 1 # Interdiction Budget u = [ # Capacities of edges # 1 2 3 4 Nodes 0 10 10 0 #1 0 0 0 M #2 # Dummy flow to network sink 0 0 0 M #3 # Dummy flow to network sink M 0 0 0 #4 # Dummy flow from network sink to source ] m = size(u) # Number of nodes # Define edges susceptible to interdiction for future scenarios Itilda1 = zeros(4,4) Itilda1[1,2]=1 # Node 1 to Node 2 Itilda2 = zeros(4,4) Itilda2[1,3]=1 # Node 1 to Node 3 # Define a discrete set of scenarios ξ₁ = @scenario I[i in 1:m, j in 1:m] = Itilda1 probability = .2 ξ₂ = @scenario I[i in 1:m, j in 1:m] = Itilda2 probability = .8 @stochastic_model simple_model begin @stage 1 begin # Make interdiction decisions @decision(simple_model, y[i=1:m,j=1:m], Bin) # Interdiction decisions by edge @objective(simple_model, Min, 0) # Minimize the second stage result @constraint(simple_model, sum(y) <= B) # Number of Interdictions Below Interdiction Budget end @stage 2 begin @known(simple_model, y) # Interdiction decisions by edge from stage 1 @uncertain I[i in 1:m, j in 1:m] # Realization of the interdiction effectiveness @recourse(simple_model, x[i=1:m,j=1:m] >= 0) # Flow decisions by edge @objective(simple_model, Max, sum(x[1:m,4])) # Maxmize flow @constraint(simple_model, [i = 1:m, j = 1:m], x[i, j] <= u[i, j]*(1-(y[i,j]*I[i,j]))) # Capacity constraint @constraint(simple_model, [i = 1:m], sum(x[i, 1:m]) == sum(x[1:m, i])) # Flow balance constraint end end sp = instantiate(simple_model, [ξ₁, ξ₂], optimizer = LShaped.Optimizer) set_optimizer_attribute(sp, MasterOptimizer(), GLPK.Optimizer) set_optimizer_attribute(sp, SubProblemOptimizer(), GLPK.Optimizer) optimize!(sp) optimal_decision(sp) objective_value(sp)