Julia: Run stochastic models in parallel inside a for loop

Question

It's my first time posting here, so I apologize in advance for any error in the writing of the post.

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!