It looks like this:
using JuMP
using GLPK
using MathOptInterface
using GLPKMathProgInterface
using Gurobi
n = 12
env = Gurobi.Env()
m = Model(Gurobi.Optimizer)
set_optimizer_attribute(m, "NonConvex", 2)
@variable(m, x[i = 1:n, j = 1:n] >= 0)
@variable(m, y[i = 1:n, j = 1:n], Bin)
@variable(m, 1 <= nr[i = 1:n, j = 1:n] <= n, Int)
start = [18, 15, 17, 17, 20, 13, 15, 17, 17, 21, 18, 15]
pop = [100, 80, 130, 90, 1500, 1200, 200, 250, 800, 1000, 1100, 750]
growth = [0.2, 0.15, 0.22, 0.24, 0.12, 0.42, 0.1, 0.11, 0.14, 0.2, 0.05, 0.05]
mgrowth = zeros(n,n)
for i = 1:n
for j = 1:n
mgrowth[i,j] = (1 + growth[i])^(n-j)
end
end
@objective(m, Max, sum(x[i,j]*y[i,j] * pop[j] * mgrowth[i,j] for i =1:n, j = 1:n) )
for j = 1:n
@constraint(m, sum(y[i,j] for i = 1:n) == 1 )
end
for i = 1:n
@constraint(m, sum(y[i,j] for j = 1:n) == 1 )
end
for i = 1:n
@constraint(m, x[i,1]*y[i,1] <= start[i] )
end
for j = 2:n
@constraint(m, nr[1,j] >= sum(y[2,i] + y[3,i] for i = 1:j-1) )
@constraint(m, nr[1,j] * x[1,j] == sum(y[2,i]*x[2,i] + y[3,i]*x[3,i] for i = 1:j-1) )
@constraint(m, nr[2,j] >= sum(y[3,i] + y[1,i] + y[12,i] for i = 1:j-1) )
@constraint(m, nr[2,j] * x[2,j] == sum(y[3,i]*x[3,i] + y[1,i]*x[1,i] + y[12,i]*x[12,i] for i = 1:j-1) )
@constraint(m, nr[3,j] >= sum(y[2,i] + y[1,i] + y[4,i] + y[5,i] for i = 1:j-1) )
@constraint(m, nr[3,j] * x[3,j] == sum(y[2,i]*x[2,i] + y[1,i]*x[1,i] + y[4,i]*x[4,i] + y[5,i]*x[5,i] for i = 1:j-1) )
@constraint(m, nr[4,j] >= sum(y[1,i] + y[2,i] + y[3,i] for i = 1:j-1) )
@constraint(m, nr[4,j] * x[4,j] == sum(y[1,i]*x[1,i] + y[2,i]*x[2,i] + y[3,i]*x[3,i] for i = 1:j-1) )
@constraint(m, nr[5,j] >= sum(y[1,i] + y[2,i] + y[3,i] + y[4,i] + y[10,i] + y[11,i] + y[12,i] for i = 1:j-1) )
@constraint(m, nr[5,j] * x[5,j] == sum(y[1,i]*x[1,i] + y[2,i]*x[2,i] + y[3,i]*x[3,i] + y[4,i]*x[4,i]
+ y[10,i]*x[10,i] + y[11,i]*x[11,i] + y[12,i]*x[12,i] for i = 1:j-1) )
@constraint(m, nr[6,j] >= sum(y[2,i] + y[12,i] for i = 1:j-1) )
@constraint(m, nr[6,j] * x[6,j] == sum(y[2,i]*x[2,i] + y[12,i]*x[12,i] for i = 1:j-1) )
@constraint(m, nr[7,j] >= sum(y[8,i] for i = 1:j-1) )
@constraint(m, nr[7,j] * x[7,j] == sum(y[8,i]*x[8,i] for i = 1:j-1) )
@constraint(m, nr[8,j] >= sum(y[7,i] for i = 1:j-1) )
@constraint(m, nr[8,j] * x[8,j] == sum(y[7,i]*x[7,i] for i = 1:j-1) )
@constraint(m, nr[9,j] >= sum(y[10,i] + y[11,i] + y[12,i] for i = 1:j-1) )
@constraint(m, nr[9,j] * x[9,j] == sum(y[10,i]*x[10,i] + y[11,i]*x[11,i] + y[12,i]*x[12,i] for i = 1:j-1) )
@constraint(m, nr[10,j] >= sum(y[9,i] + y[11,i] for i = 1:j-1) )
@constraint(m, nr[10,j] * x[10,j] == sum(y[9,i]*x[9,i] + y[11,i]*x[11,i] for i = 1:j-1) )
@constraint(m, nr[11,j] >= sum(y[7,i] + y[8,i] + y[9,i] + y[10,i] + y[12,i] for i = 1:j-1) )
@constraint(m, nr[11,j] * x[11,j] == sum(y[7,i]*x[7,i] + y[8,i]*x[8,i] + y[9,i]*x[9,i] + y[10,i]*x[10,i] + y[12,i]*x[12,i] for i = 1:j-1) )
@constraint(m, nr[12,j] >= sum(y[1,i] + y[2,i] + y[3,i] + y[4,i] + y[5,i] + y[6,i] + y[7,i] + y[8,i] + y[9,i] for i = 1:j-1) )
@constraint(m, nr[12,j] * x[12,j] == sum(y[1,i]*x[1,i] + y[2,i]*x[2,i] + y[3,i]*x[3,i] + y[4,i]*x[4,i]
+ y[5,i]*x[5,i] + y[6,i]*x[6,i] + y[7,i]*x[7,i] + y[8,i]*x[8,i] + y[9,i]*x[9,i] for i = 1:j-1) )
end
optimize!(m)
println(JuMP.value.(x))
println(JuMP.value.(y))
println(JuMP.value.(nr))
println(JuMP.objective_value(m))
println(JuMP.termination_status(m))
As you can see, I added a mgrowth to the objective which is a shop-specific monthly growth in the customer size.