Hi everyone,

I’m working in my master’s thesis and trying to lower the running time that my code has, by parallelizing the optimization of different models created by scenarios, to then calculate some other values after saving the results of each model. My thesis model is to large to post it here, but I’ve created an example that does exactly the same in a much lower scale.

The aim is to optimize a stochastic model using the Progressive Hedging algorithm. It’s worth noting that no current library can help my because of the size of my model, so I had to implement the algorithm from scratch.

The example is this one:

```
using JuMP, Gurobi, DataFrames
const GUROBI_ENV = Gurobi.Env()
I = DataFrame(MatPrima = ["Madera", "Acabado", "Carpinetria"], Costo=[2,4,5.2], Escritorios=[8,4,2], Mesas=[6,2,1.5], Sillas=[1,1.5,0.5], Precio = [60,40,10])
D = DataFrame(Escenario = ["Bajo", "Medio", "Alto", "Promedio"], Escritorio = [50,150,250,150], Mesas = [20,105,250,125], Sillas = [180,220,500,300] )
# Precios
p = I[:,"Precio"]
# Costos
c = I[:,"Costo"]
testoBj = zeros(3)
testX1 = zeros(3,3)
x_barra = [0.0, 0.0, 0.0]
x_b_record = zeros(3)
X_record = zeros(3,3)
OBJ_Record = zeros(3)
dc1 = [1.0,0.0,-1.0] #zeros(3)
dc2 = [1.0,0.0,-1.0] #zeros(3)
dc3 = [1.0,0.0,-1.0] #zeros(3)
act_lambda = zeros(3,3)
lambda_record = zeros(3,3)
xb_mod = 0
l_mod = 0
xmod_record = 0
lmod_record = 0
limite = 0
lim_record = 0
global alpha = 0.001
Niter = 1000
for iter in 1:Niter
# Creating model for scenarios
for sc in 1:3 # 1: low, 2: medium, 3: high
dakota_PH = Model(() -> Gurobi.Optimizer(GUROBI_ENV))
set_optimizer_attribute(dakota_PH, "Method", 2)
set_optimizer_attribute(dakota_PH, "OutputFlag", 0)
# Variables
@variable(dakota_PH, x[1:3] >= 0) # wood, finish, carpentry
@variable(dakota_PH, y[1:3] >= 0) # desks, tables, chairs
# Constraints
@constraint(dakota_PH, c1[c in 1:3], x[c] >= sum(I[c,j+2]*y[j] for j in 1:3))
@constraint(dakota_PH, c2[i in 1:3], y[i] <= D[sc,i+1])
# Objective function
@objective(dakota_PH, Min, - ( sum(p[i]*y[i] - c[i]*x[i] for i in 1:3) )
+ (dc1[sc]*x[1] + dc2[sc]*x[2] + dc3[sc]*x[3])
+ (alpha/2)*sum( (x[i]-x_barra[i])^2 for i in 1:3)
)
# solving
optimize!(dakota_PH)
# Recording objective and results (X) of each scenary
testoBj[sc] = objective_value(dakota_PH)
for i in 1:3
testX1[sc,i] = value(x[i])
end
end
# Creating x_barra
x_barra = (testX1[1,:].+testX1[2,:].+testX1[3,:] )./3 # Promedio de X
for sc in 1:3 # Rows/Scenarios Lambda matrix
dc1[sc] = dc1[sc] + alpha*(testX1[sc,1]-x_barra[1])
dc2[sc] = dc2[sc] + alpha*(testX1[sc,2]-x_barra[2])
dc3[sc] = dc3[sc] + alpha*(testX1[sc,3]-x_barra[3])
end
act_lambda = reshape([dc1';dc2';dc3'],3,3) # ACTUAL lambda matrix
x_b_record = [x_b_record ; x_barra] # Saving x_barra
X_record = [X_record ; testX1] # Saving X
OBJ_Record = [OBJ_Record ; testoBj] # Saving Objective
lambda_record = [lambda_record ; act_lambda] # Saving Lambda
# Calculation of the Limit
xb_prev_and_act = last(x_b_record,6)
xb_prev = first(xb_prev_and_act,3) # Previous x_barra
xb_mod = sum((x_barra .- xb_prev).^2) # act - prev # Module X_Barra
prev_l = reshape([lambda_record[((size(lambda_record)[1])-5),:]; # Previous Lamda
lambda_record[((size(lambda_record)[1])-4),:];
lambda_record[((size(lambda_record)[1])-3),:]],(3,3))'
l_mod = sum((act_lambda .- prev_l).^2) # act - prev # Module Lambdas
limite = xb_mod + (1/alpha^2)*l_mod
xmod_record = [xmod_record ; xb_mod] # Saving xb_mod
lmod_record = [lmod_record ; l_mod] # Saving l_mod
lim_record = [lim_record ; limite] # Saving limit
end
```

I’m not fully aware how to achieve this parallelization using the Distributed.jl library. I know it will be helpful to send the models by scenary to different threads, but this is when I have no idea how to do it properly to get its results on a “SharedArray” (or Dicts, as in my thesis I’m using Dicts to save the variable values, as they have 3 to 4 dimensions with specific indeces).

If anyone has any clue or maybe can help me in any form, I’ll be very grateful.

Thanks!