JuMP Optimization Help

Specific Domains Optimization (Mathematical)
1

I’m brand new to JuMP.jl and optimization and I’m running into a bit of a roadblock figuring out how to solve what I thought would be a fairly straightforward problem. Imagine 6 jars of marbles, divided into two groups. Each jar has a different number of marbles in it. I have purchased 9 new marbles and I would like to distribute them amongst the 6 jars in a way that helps even out the number of marbles in each jar, subject to some constraints:

  1. I cannot put more than 3 marbles in any one jar.

  2. As the jars are divided into two groups, I cannot allocate more than 6 of the new marbles to a group.

  3. I must use all of the 9 new marbles.

  4. I cannot move marbles from one jar to another.

The goal is to, after some number of rounds, end up with the same number of marbles in each jar. Ultimately, I want to know how many rounds of purchasing new marbles/distributing them amongst the jars it will take to have the same number of marbles in each jar, given the constraints.

The trainwreck code below was my first attempt at optimizing a single round of this experiment :blush::

using GLPK
using JuMP
using Statistics

# number of new marbles
N = 9
# numer of marbles currently in each of the 3 jars in group 1
m₁ = [5,3,2]
# numer of marbles currently in each of the 3 jars in group 2
m₂ = [5,10,12]

model = Model(GLPK.Optimizer)

# set up x₁ integer-valued variables for newly-allocated group 1 marbles
@variable(model, x₁[1:3], Int)
# set up x₂ integer-valued variables for newly-allocated group 2 marbles
@variable(model, x₂[1:3], Int)

# the objective is to minimize the standard deviation of all jars
@objective(model, Min, std(vcat(m₁ .+ x₁, m₂ .+ x₂)))

# we must use all of the new marbles
@constraint(model, sum(vcat(x₁, x₂)) == N)
# no single jar can receive more than 3 new marbles
@constraint(model, 0 .<= x₁ .<= 3)
@constraint(model, 0 .<= x₂ .<= 3)
# neither group can receive more than 6 of the new marbles
@constraint(model, sum(x₁) <= 6)
@constraint(model, sum(x₂) <= 6)

Executing this code results in this error:

ERROR: abs2 is not defined for type GenericAffExpr. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective.

And switching to @NLobjective leads to a whole new set of errors/problems. I’m hoping someone more experienced in this realm can steer me in the right direction on this. Thanks so much!

2

GLPK is able to handle the integrality constraints that you defined, however it cannot handle the non-linear objective that you want.

As expressed by the error message, nonlinear objectives must be entered with the other macro.

Any reason to that objective?

You have a few options:

Easier to optimize:

  1. change you objective function to something linear to keep using GLPK, you probably can model a absolute value related objective using linear programming techniques.
    like minimizing the largest absolute distance between two jars. or even minimize the number of marbles in the jar with most marbles…

Harder to optimize:

  1. change the solver to something that handles non-linear and integer (this class of problems is extremely hard to solve to optimality), some options are: Alpine.jl, Juniper.jl, Knitro, SCIP, Couenne (via GitHub - jump-dev/AmplNLWriter.jl: Julia interface to AMPL-enabled solvers)

  2. you might be able to obtain a reformulation or an approximation of your objective using quadratic terms so you can use Gurobi, Xpress or CPLEX

  3. you might be able to obtain a reformulation or an approximation of your objective using cones so you can use Mosek

There are certainly options and solvers that I am missing…

3 Likes
3

Any reason to that objective?

It’s just the first thing that occurred to me. I see though that trying to find a linear objective function that will serve a similar purpose seems like the best way to go so I’m going to have a go at that. Thank you very much for your feedback, I really appreciate it :slightly_smiling_face:

4

like minimizing the largest absolute distance between two jars.

I’m attempting to do this and cannot figure it out :confounded::

# number of new marbles
N = 9
# numer of marbles currently in each of the 3 jars in group 1
m₁ = [5,3,2]
# numer of marbles currently in each of the 3 jars in group 2
m₂ = [5,10,12]

# concatenate m₁ and m₂
m = hcat(m₁, m₂)

model = Model(GLPK.Optimizer)

# a matrix instead of two different variables
@variable(model, x[1:3, 1:2], Int)

fmax(y) = maximum(y)
fmin(y) = minimum(y)

register(model, :fmax, 1, fmax, autodiff=true)
register(model, :fmin, 1, fmin, autodiff=true)

# the objective is to minimize the spread between the max and min
julia> @objective(model, Min, fmax(m + x) - fmin(m + x))
ERROR: MethodError: no method matching isless(::GenericAffExpr{Float64,VariableRef}, ::GenericAffExpr{Float64,VariableRef})
Closest candidates are:
  isless(::Missing, ::Any) at missing.jl:87
  isless(::Any, ::Missing) at missing.jl:88
5

GLPK does not accept user defined functions.
You have to write that as a set of linear inequality constraints. For a simpler case:

@variable(model, x[1:3], Int)
@variable(model, y >= 0)
for (i,j) in [(1,2), (2,3), (3,1)]
    @constraint(model, x[i] - x[j] <= y)
    @constraint(model, x[j] - x[i] <= y)
end
@objective(model, Min, y)

y will be larger or equal to abs(x[j] - x[i]) since you are minimizing, equality will hold

1 Like
6

Thank you for your assistance. I realized that dropping the requirement that x be an Int really makes my life easier. This example is a simplified approximation of a real-world problem I’m trying to solve and in my real-world problem, it’s not absolutely critical that the algorithm allocates all 9 of the new marbles - if there’s one leftover I can just save it for next round and if it needs 1 extra, I can always buy 1 more marble :wink:. That being said, I’ve settled on something like this (for now):

using Ipopt
using JuMP

N = 9
m₁ = [5,3,2]
m₂ = [5,10,12]
m = hcat(m₁, m₂)

model = Model(Ipopt.Optimizer)

@variable(model, x[1:3,1:2])

@objective(model, Min, sum(((m+x) .- mean(m+x)).^2))
@constraint(model, 0 .<= x .<= 3)
@constraint(model, sum(x) <= 9)
@constraint(model, sum(x[:,1]) <= 6)
@constraint(model, sum(x[:,2]) <= 6)

optimize!(model)
objective_value(model)

round.(value.(x))

julia> round.(value.(x))
3×2 Array{Float64,2}:
 1.0   3.0
 3.0  -0.0
 3.0  -0.0

julia> sum(round.(value.(x)))
10.0

# And to show that the variance has been reduced:
julia> var(m)
15.766666666666666

julia> var(m .+ value.(x))
8.139999918409327

In this case, I can round the values of x to see how many marbles to add to each jar and I end up using 10 marbles total. That’s okay though, I usually have a spare marble laying around or, as mentioned, I can always go buy an extra one to make it all work. And, with that, I’m signing off to reward myself with a nice :beer: because I cannot believe I’ve spent all day attempting to optimize a marbles-in-jars problem :rofl:

1 Like