Objective function value mysteriously negative with StochasticPrograms.jl

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)

Hi @matt_w, welcome to the forum.

Are you trying to solve a min-max problem? I don’t think StochasticPrograms is setup for that. Why wouldn’t an optimal solution just be for y = 0?

You might want to try GitHub - joaquimg/BilevelJuMP.jl: Bilevel optimization in JuMP instead.

You’ll need to define one lower level, which contains every scenario and the objective is the expected cost: Multiple second level problems · Issue #1 · joaquimg/BilevelJuMP.jl · GitHub

Hi @odow, thanks for the quick response.

You are correct, I am attempting to code a min-max problem. I thought StochasticPrograms.jl would work because of the example of one used in the documentation here: https://martinbiel.github.io/StochasticPrograms.jl/dev/manual/decisions/

I am expecting y[1,3]=1 to be the optimal decision in my code above as it should reduce the max-flow capable in the second-stage most of the time.

I’ll take a look at your suggestion using BilevelJuMP.jl to solve the discrete scenario case, but I was hoping to use StochasticPrograms.jl so I could extend my code to sample from a distribution instead of a discrete number of predefined scenarios: https://martinbiel.github.io/StochasticPrograms.jl/dev/manual/data/#Sampling

It looks like StochasticPrograms just converts the Max into a Min by multiplying the objective by -1. that likely explains your negative objective.

But that’s not quite the same model as the adversarial problem you’re trying to solve, so I don’t think you can use StochasticPrograms here (regardless of the possible options or modeling tricks).

1 Like