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

x
contains integers between 0 and 5 (inclusive)  For each column
col
ofx
:
1 <= col[i]  col[i2] <= 2
fori = 3, 4, 7, 8

col[i] % 3 == 0
fori = 1, 2, 5, 6

 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 columnwise 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[i2,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?