Constrain 2D array variable to contain a value `x` number of times

Specific Domains Optimization (Mathematical)
1

I am trying to solve for a 2D array x of dimension 4x2 that satisfy the following:

  1. x contains integers between 0 and 5 (inclusive)
  2. For each column col of x:
    • 1 <= col[i] - col[i-2] <= 2 for i = 3, 4, 7, 8
    • col[i] % 3 == 0 for i = 1, 2, 5, 6
  3. Values that are a multiple of 3 appear twice, while all other values appear once in x.

The goal is to maximize the average (across columns) of the column-wise variance of x. Here is my attempt so far:

using JuMP
using Statistics
using Gurobi

rows = 4
cols = 2
possible_vals = 6
model = Model(optimizer_with_attributes(Gurobi.Optimizer))

@variable(model, 0 <= x[1:rows, 1:cols] <= 5, Int)
@variable(model, z[1:rows, 1:cols] >= 0, Int)

for j=1:cols
    for i=1:rows
        if (i - 1) % 4 >= 2
            @constraint(model, 1 <= x[i,j] - x[i-2,j] <= 2)
        else
            @constraint(model, x[i,j] - (3 * z[i,j]) == 0)
        end
    end
end

@variable(model, y[1:possible_vals], Int)
for i=1:possible_vals
    if (i - 1) % 3 == 0
        @constraint(model, y[i] == 2)
    else
        @constraint(model, y[i] == 1)
    end
    @constraint(model, y[i] == sum(x .== (i - 1)))
end

@expression(model, f, mean([var(x[:,j]) for j in 1:cols]))
@objective(model, Max, f)
optimize!(model)

The last constraint is causing the problem to be infeasible. Is there a different way I can formulate this?

2

I’d initially model it like this:

using JuMP
model = Model()
rows, cols, values = 4, 2, 0:5
@variables(model, begin
    x[1:rows, 1:cols, values], Bin
    z[1:rows, 1:cols] >= 0, Int
    σ²[1:cols]
    μ[1:cols]
end)
@expression(model, X[i=1:rows, j=1:cols],  sum(k * x[i, j, k] for k in values))

@constraints(model, begin
    # Must pick one value for each square
    [i=1:rows, j=1:cols], sum(x[i, j, k] for k in values) == 1
    # (2): I don't really know what these represent
    [j=1:cols, i=1:rows; (i - 1) % 4 >= 2], 1 <= X[i,j] - X[i-2,j] <= 2
    [j=1:cols, i=1:rows; (i - 1) % 4 < 2], X[i,j] == 3 * z[i,j]
    # (3): Values that are a multiple of 3 appear twice
    [k=values], sum(x[:, :, k]) == 1 + (mod(k, 3) == 0 ? 1 : 0)
    # Compute mean of each column
    [j=1:cols], μ[j] == sum(X[i, j] for i in 1:rows) / rows
    # Compute variance
    [j=1:cols], σ²[j] == sum((X[i, j] - μ[j])^2 for i in 1:rows) / (rows - 1)
end)
@objective(model, Max, sum(σ²) / cols)

but you have aa problem. Variance is a convex quadratic, and you want to maximize it, which makes this problem non-convex. That leaves you with really only one solver choice, for which you’ll need a license:
(Edit: I see you just edited your post to use Gurobi, so I guess you’re in luck.)

import Gurobi
set_optimizer(model, Gurobi.Optimizer)
set_optimizer_attribute(model, "NonConvex", 2)
optimize!(model)

julia> value.(X)
4×2 Matrix{Float64}:
 3.0  3.0
 0.0  0.0
 4.0  5.0
 2.0  1.0
1 Like
3

I just saw that Statistics.mean and Statistics.var work on JuMP variables, so this simplifies to

using JuMP, Statistics, Gurobi
rows, cols, values = 4, 2, 0:5
model = Model(Gurobi.Optimizer)
set_optimizer_attribute(model, "NonConvex", 2)
@variable(model, x[1:rows, 1:cols, values], Bin)
@variable(model, z[1:rows, 1:cols] >= 0, Int)
@expression(model, X[i=1:rows, j=1:cols],  sum(k * x[i, j, k] for k in values))
@constraints(model, begin
    [i=1:rows, j=1:cols], sum(x[i, j, k] for k in values) == 1
    [j=1:cols, i=1:rows; (i - 1) % 4 >= 2], 1 <= X[i,j] - X[i-2,j] <= 2
    [j=1:cols, i=1:rows; (i - 1) % 4 < 2], X[i,j] == 3 * z[i,j]
    [k=values; mod(k, 3) == 0], sum(x[:, :, k]) == 2
    [k=values; mod(k, 3) != 0], sum(x[:, :, k]) == 1
end)
@objective(model, Max, mean(var(X[:, j]) for j in 1:cols))
optimize!(model)
value.(X)
1 Like
4

Thank you for the reformulation (and for adjusting to my discovery of Gurobi after posting), this is so compact and clever!

1 Like