Strategies for dealing with many indices

This tutorial was generated using Literate.jl. Download the source as a .jl file.

Common modeling patterns may lead to variables that are defined over many indices. For example, a transshipment model may have x[p, f, m, t] to represent the flow of product p from factory f to market m in time t. If each set is represented by a range of integers, it can be easy to make little errors like x[t, p, f, m] that are hard to catch and debug. The purpose of this tutorial is to explain how to used typed indices to avoid such bugs.

Required packages

This tutorial uses the following packages:

julia> using JuMP

The problem

Consider a transshipment model of a single product from factories to markets. One way to model the flow variable is:

julia> F, M = 3, 3(3, 3)
julia> model = Model()A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none
julia> @variable(model, x[1:F, 1:M] >= 0)3×3 Matrix{VariableRef}:
 x[1,1]  x[1,2]  x[1,3]
 x[2,1]  x[2,2]  x[2,3]
 x[3,1]  x[3,2]  x[3,3]

Now the question is: what does x[1, 2] represent? Was it x[f, m] or x[m, f]? In this simple example, it's probably easy to remember the flow from a factory to a market, but it gets harder if there are more sets. Was time the first index or the last?

Mis-typing the order of two sets is a common error, and it can result in bugs that are surprisingly hard to find and debug. If you didn't already know, would you be able to spot the difference in these two constraints, especially if the constraint expression was large or you were looking through many constraints trying to find out why the solution wasn't what you expected?

julia> @constraint(model, [f in 1:F], sum(x[f, m] for m in 1:M) == 1)3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[1,1] + x[1,2] + x[1,3] = 1
 x[2,1] + x[2,2] + x[2,3] = 1
 x[3,1] + x[3,2] + x[3,3] = 1

julia> @constraint(model, [f in 1:F], sum(x[m, f] for m in 1:M) == 1)3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[1,1] + x[2,1] + x[3,1] = 1
 x[1,2] + x[2,2] + x[3,2] = 1
 x[1,3] + x[2,3] + x[3,3] = 1

There are two approaches you can use to avoid bugs like x[m, f]: keyword indexing and typed indices.

Keyword indexing

The first approach is to use keyword indexing. Because keyword indexing does not work for Array containers, we must explicitly choose the DenseAxisArray container.

julia> F, M = 3, 3(3, 3)
julia> model = Model()A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none
julia> @variable(model, x[f in 1:F, m in 1:M] >= 0, container = DenseAxisArray)2-dimensional DenseAxisArray{VariableRef,2,...} with index sets:
    Dimension 1, Base.OneTo(3)
    Dimension 2, Base.OneTo(3)
And data, a 3×3 Matrix{VariableRef}:
 x[1,1]  x[1,2]  x[1,3]
 x[2,1]  x[2,2]  x[2,3]
 x[3,1]  x[3,2]  x[3,3]

Accessing x in the correct order works:

julia> x[f = 1, m = 2]x[1,2]

But the wrong order errors:

julia>     x[m = 2, f = 1]Invalid index m in position 1. When using keyword indexing, the indices must match the exact name and order used when creating the container.

Keyword indexing means we can write constraints like this:

julia> @constraint(model, [f in 1:F], sum(x[f = f, m = m] for m in 1:M) == 1)3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[1,1] + x[1,2] + x[1,3] = 1
 x[2,1] + x[2,2] + x[2,3] = 1
 x[3,1] + x[3,2] + x[3,3] = 1

Typed indices

A second approach is to define new types for each index:

julia> struct Factory
           x::Int
       end
julia> struct Market
           x::Int
       end

Then, we define each axis as a vector of these types:

julia> factories = Factory.(1:F)3-element Vector{Main.Factory}:
 Main.Factory(1)
 Main.Factory(2)
 Main.Factory(3)
julia> markets = Market.(1:M)3-element Vector{Main.Market}:
 Main.Market(1)
 Main.Market(2)
 Main.Market(3)
julia> model = Model()A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none
julia> @variable(model, x[factories, markets] >= 0)2-dimensional DenseAxisArray{VariableRef,2,...} with index sets:
    Dimension 1, Main.Factory[Main.Factory(1), Main.Factory(2), Main.Factory(3)]
    Dimension 2, Main.Market[Main.Market(1), Main.Market(2), Main.Market(3)]
And data, a 3×3 Matrix{VariableRef}:
 x[Main.Factory(1),Main.Market(1)]  …  x[Main.Factory(1),Main.Market(3)]
 x[Main.Factory(2),Main.Market(1)]     x[Main.Factory(2),Main.Market(3)]
 x[Main.Factory(3),Main.Market(1)]     x[Main.Factory(3),Main.Market(3)]

Accessing x in the correct order works:

julia> x[Factory(1), Market(2)]x[Main.Factory(1),Main.Market(2)]

But the wrong order errors:

julia>     x[Market(2), Factory(1)]KeyError: key Main.Factory(1) not found

Typed indices means we can write constraints like this:

julia> @constraint(model, [f in 1:F], sum(x[Factory(f), Market(m)] for m in 1:M) == 1)3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[Main.Factory(1),Main.Market(1)] + x[Main.Factory(1),Main.Market(2)] + x[Main.Factory(1),Main.Market(3)] = 1
 x[Main.Factory(2),Main.Market(1)] + x[Main.Factory(2),Main.Market(2)] + x[Main.Factory(2),Main.Market(3)] = 1
 x[Main.Factory(3),Main.Market(1)] + x[Main.Factory(3),Main.Market(2)] + x[Main.Factory(3),Main.Market(3)] = 1

or like this:

julia> @constraint(model, [f in factories], sum(x[f, m] for m in markets) == 1)1-dimensional DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape},1,...} with index sets:
    Dimension 1, Main.Factory[Main.Factory(1), Main.Factory(2), Main.Factory(3)]
And data, a 3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[Main.Factory(1),Main.Market(1)] + x[Main.Factory(1),Main.Market(2)] + x[Main.Factory(1),Main.Market(3)] = 1
 x[Main.Factory(2),Main.Market(1)] + x[Main.Factory(2),Main.Market(2)] + x[Main.Factory(2),Main.Market(3)] = 1
 x[Main.Factory(3),Main.Market(1)] + x[Main.Factory(3),Main.Market(2)] + x[Main.Factory(3),Main.Market(3)] = 1