Hello,

I am working with an integer optimization problem and I am looking for ways to speed the code up. The math problem is the following

I am sampling K different outcomes of the random variable \omega^{k} which in the example is msimul and solving

\max_{I(\omega^{k})\forall k} \frac{1}{K} \sum_{k=1}^{K} F(I(\omega^{k}),\omega^{k})

subject to

I(\omega^{1})=I(\omega^{2})=I(\omega^{3})=\cdots=I(\omega^{K})

Find below a MWE

```
using JuMP, BenchmarkTools, Gurobi, Random, Distributions
S = 100
J = 77
Ι = 150 #this is an Iota \Iota #
Random.seed!(12)
σ = 3
ϕ = 1
κ = 0.95
mtemp = rand(Normal(0,1),Ι,J,S)
msimul = exp.(mtemp)
β = ones(J)
F=rand(Ι)*100000
function problem(msimul,σ,F,J,Ι,κ,ϕ,β)
m = (msimul).^(1-σ)
model = Model(Gurobi.Optimizer)
K = copy(S)
@variable(model, Ipost[1:Ι,1:J,1:K], Bin)
@variable(model, Itrue[1:Ι], Bin)
set_silent(model)
# Objective Function
println("Objective")
@objective(model,
Max,
ϕ/K*sum(β[j]*(m[:,j,k]'*Ipost[:,j,k]) for j=1:J for k=1:K) - F'*Itrue )
println("Constraints")
# Constraint 1 and 2 - non-anticipativity constraint and ex-post choice
for k=1:K
for j=1:J
@constraint(model, Ipost[:,j,k]'*Ipost[:,j,k] == 1)
@constraint(model, sum(Ipost[:,j,k]'*Itrue ) == 1)
end
end
println("Optimize")
optimize!(model)
Ipost = value.(Ipost)
Ipost = trunc.(Int,Ipost)
Iopt = value.(Itrue)
return trunc.(Int,Iopt), objective_value(model)
end
solution = @btime problem(msimul,σ,F,J,Ι,κ,ϕ,β)
```