Run stochastic models in parallel in for loop

Specific Domains
1

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!

2
3

Hi Pablo!

First of all, I’m not sure you want to use distributed computing (with separate memory), maybe multithreading (with shared memory) is closer to what you need here.

Now onto storage of results: the simplest way is to create a vector of dicts or named tuples that each thread can fill, like so:

results = Vector{Dict{Int,String}}(undef, 10)
Threads.@threads for i in 1:10
    d = Dict(1 => "hello", 2 => "goodbye")
    results[i] = d
end

As long as your threads do not write to the same array element, you don’t need a thread-safe data structure. Note however that dictionaries are not thread-safe by default (even if your threads write to different keys!).

Finally, before you use multithreading, be sure to optimize your sequential code as much as possible following the official performance tips or my own summary. If your sequential code is not efficient, multithreading will not bring a significant speedup.