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)[1] # 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)