Network flow problem in JuMP

Specific Domains Optimization (Mathematical)
1

Dear all,
I have the shortest path problem coded in JuMP as follows:

using JuMP
import HiGHS
import LinearAlgebra
G = [
    0 100 30 0 0
    0 0 20 0 0
    0 0 0 10 60
    0 15 0 0 50
    0 0 0 0 0
]

n = size(G)[1]

shortest_path = Model(HiGHS.Optimizer)

@variable(shortest_path, x[1:n, 1:n], Bin)
@constraint(shortest_path, [i = 1:n, j = 1:n; G[i, j] == 0], x[i, j] == 0)
@constraint(
    shortest_path,
    [i = 1:n; i != 1 && i != 2],
    sum(x[i, :]) == sum(x[:, i])
)
@constraint(shortest_path, sum(x[1, :]) - sum(x[:, 1]) == 1)
@constraint(shortest_path, sum(x[2, :]) - sum(x[:, 2]) == -1)
@objective(shortest_path, Min, LinearAlgebra.dot(G, x))

optimize!(shortest_path)

that works perfectly. However, I would like to re-write this code creating only the variables where a arc existis, e.g., not create n x n binay variables and set to zero all where G[i,j] =0.
My idea is determine the indexes i and j whose G[i,j] >0, as follows:

ind = findall(G .> 0)
7-element Array{CartesianIndex{2},1}:
 CartesianIndex(1, 2)
 CartesianIndex(4, 2)
 CartesianIndex(1, 3)
 CartesianIndex(2, 3)
 CartesianIndex(3, 4)
 CartesianIndex(3, 5)
 CartesianIndex(4, 5)

That’s ok. But, now, how can I create the variable x[i,j] where i and j belongs to aray ind So I create only the necessary variables and constraints…
Thank ypou so much for the help.

2

Use something like:

ind = [(i, j) for i in 1:size(G, 1), j in 1:size(G, 2) if G[i, j] > 0]
model = Model()
@variable(model, x[ind])
x[(1, 2)]

or

model = Model()
@variable(model, x[i=1:size(G,1), j=1:size(G,2); G[i,j] > 0])
x[1, 2]
1 Like
3

Thank you very much.
How can I write the constraint: the sum of x over Index i equal to 1 (for each j), for example If I use the option 1 or 2?

4

You could use

ind = [(i, j) for i in 1:size(G, 1), j in 1:size(G, 2) if G[i, j] > 0]
model = Model()
@variable(model, x[ind])
@constraint(
    model, 
    [j=1:size(G, 2)], 
    sum(x[(i, j)] for (i, jj) in ind if jj == j) == 1,
)


model = Model()
@variable(model, x[i=1:size(G,1), j=1:size(G,2); G[i,j] > 0])
@constraint(
    model, 
    [j=1:size(G,2)], 
    sum(x[i, j] for (i, jj) in eachindex(x) if jj == j) == 1,
)
1 Like