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

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?

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

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

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

1 Like