I am currently working on solving a optimization problem using julia and gurobi. I need to solve this optimizaton problem around 500 times (each time uses a different set of data). However, instead of running a julia script file manually for each set of different data, I currently have a for loop which changes the data and calls a function (which I created myself) to perform the optimization. Every time a new iteration starts starts in the loop, the same julia session continues to be used. How can I start a new julia session for ever iteration in the for loop. I would like each iteration to start fresh. Also I am running my scripts using VS code.
You cannot start a new Julia section from within the Julia section. You could use DaemonMode.jl or loop in a bash script, but that would be unreasonable considering how Julia works.
But probably you don’t need to restart the Julia section at all, if you explain better what you want to do you might get some more proper advice.
Ok, thanks for reply. I am essentially trying to do a sensitivity analysis study using Julia (to build the model) and Gurobi (as the solver). So I have a range of values which need to be changed for the coefficients of the variables in the objective function. However instead of changing the value of the coefficients manually and then finding the solution, I have opted to have a for loop which changes the values of the coefficients for me and then find the solution. However, I am getting some weird results doing it this way. However, if I were to run over the simulations manually for those cases which produce weird results and do each run in a new julia session, the results now resemble what it should look like. So I wanted to see how I can start a new julia session in each iteration of the loop to see if the weird results go away. I am thinking along the lines of something in julia workspace is causing this. However, I noticed that you can’t clear the workspace again and also I have set a time limit for the gurobi solver (its a large MIQP problem) and don’t want anything affecting the solver when a new session should begin.
Looks like a use case for multiprocessing (or even multithreading).
If you put everything inside a function and run it, without any global variable, each run should be independent of the others. Can you post your code?
4 Likes
#options:
# 1 - Minimum Up and Down Time
# 2 - Startup Cost
# 3 - Shutdown Cost
# 4 - Ramp Rates
# 5 - DLE Gap
# 6 - Startup Time
# 7 - Gamma
# 8 - Beta
# 9 - Alpha
#SET OPTION HERE
#------------------
option = 3
opt = option
#-----------------
using JuMP
#using CPLEX
using Gurobi
using CSV
using DataFrames
using XLSX
#MIQP Function
#------------------------
function MIQP(nUnits, nHours, a, b, c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
#MODEL CONSTRUCTION
#--------------------
uc_miqp = Model(Gurobi.Optimizer)
set_optimizer_attribute(uc_miqp, "TimeLimit", 600)
#VARIABLES
#--------------------
@variable(uc_miqp, power[1:nHours, 1:nUnits])
@variable(uc_miqp, isOn[1:nHours, 1:(nZones*nUnits)] >= 0, integer=true)
@variable(uc_miqp, 0 <= startup[1:nHours, 1:nUnits] <= 1, integer=true)
@variable(uc_miqp, 0 <= rstartup1[1:nHours, 1:nUnits] <= 1, integer=true)
@variable(uc_miqp, 0 <= rstartup2[1:nHours, 1:nUnits] <= 1, integer=true)
@variable(uc_miqp, 0 <= shutdown[1:nHours, 1:nUnits] <= 1, integer=true)
#SET INITIAL value
#---------------------
for kk = 1:nHours
for jj = 1:nUnits
setvalue(power[kk, jj], init_power[kk, jj])
setvalue(startup[kk, jj], init_startup[kk, jj])
setvalue(rstartup1[kk, jj], init_rstartup1[kk, jj])
setvalue(rstartup2[kk, jj], init_rstartup2[kk, jj])
setvalue(shutdown[kk, jj], init_shutdown[kk, jj])
end
end
for kk = 1:nHours
for jj = 1:nUnits
for zz = 1:nZones
setvalue(isOn[kk, (jj + (nUnits*(zz-1)))], init_isOn[kk, (jj + (nUnits*(zz-1)))])
end
end
end
#COST EXPRESSIONS
#---------------------
#FUEL COST
#---------------------
@expression(uc_miqp, fuelcost1, sum((power.^2)*a')) #a*(P^2)
@expression(uc_miqp, fuelcost2, sum(power*b')) #b*P
@expression(uc_miqp, isOnCost, sum(isOn*(repeat(c, 1, nZones))')) #c*isOn_ll
#STARTUP COST#
#---------------------
@expression(uc_miqp, startupCost, sum(startup*ZeroStartUpCost'))
@expression(uc_miqp, rstartupCost1, sum(rstartup1*StartUpCost'))
@expression(uc_miqp, rstartupCost2, sum(rstartup2*StartUpCost'))
#SHUTDOWN COST
#----------------------
@expression(uc_miqp, shutdownCost, sum(shutdown*ShutDownCost'))
#OBJECTIVE FUNCTION
#---------------------
@objective(uc_miqp, Min, fuelcost1 + fuelcost2 + isOnCost + startupCost + rstartupCost1 + rstartupCost2 + shutdownCost)
#CONSTRAINTS
#---------------------
#Set the Upper Bounds for isOn depending on the number of operational nZones
for kk = 1:nHours
for jj = 1:nUnits
for zz = 1:nZones
@constraint(uc_miqp, isOn[kk, (jj + (nUnits*(zz-1)))] <= DLE[jj, zz])
end
end
end
#Ensure isOn_ll and isOn_hl are not equal to 1 at all times
#@constraint(uc_miqp, [kk = 1:nHours, jj = 1:nUnits], isOn_ll[kk, jj] + isOn_hl[kk, jj] <= 1)
for jj = 1:nUnits
for kk = 1:nHours
@constraint(uc_miqp, sum(isOn[kk, (jj:nUnits:(nZones*nUnits))]) <= 1)
end
end
#Load Demand constraint
#@constraint(uc_miqp, [kk = 1:nHours], sum(power[kk, :]) >= Load[kk])
for kk = 1:nHours
@constraint(uc_miqp, sum(power[kk, :]) == Load[kk])
end
#Spinning Reserve constraint
#@constraint(uc_miqp, [kk = 1:nHours], sum((maxPow.*isOn_ll[kk:kk, :])) + sum((maxPow.*isOn_hl[kk:kk, :])) >= Load[kk] + SR[kk])
for kk = 1:nHours
@constraint(uc_miqp, sum((repeat(maxPow, 1, nZones).*isOn[kk:kk, :])) >= Load[kk] + SR[kk])
end
#Generation Limit constraint
#@constraint(uc_miqp, [kk = 1:nHours, jj = 1:nUnits], power[kk, jj] <= minDLE[jj]*isOn_ll[kk, jj] + maxPow[jj]*isOn_hl[kk, jj])
#@constraint(uc_miqp, [kk = 1:nHours, jj = 1:nUnits], power[kk, jj] >= minPow[jj]*isOn_ll[kk, jj] + maxDLE[jj]*isOn_hl[kk, jj])
for jj = 1:nUnits
for kk = 1:nHours
@constraint(uc_miqp, power[kk, jj] <= sum(DLEZone[jj:jj, 2:2:(2*nZones)].*isOn[kk:kk, (jj:nUnits:(nZones*nUnits))]))
@constraint(uc_miqp, power[kk, jj] >= sum(DLEZone[jj:jj, 1:2:(2*nZones)].*isOn[kk:kk, (jj:nUnits:(nZones*nUnits))]))
end
end
#=
#Ensures that if a Unit is generating 0MW it should be off
#Needed for Units with a Minimum Generation Level of 0MW
for jj = 1:nUnits
for kk = 1:nHours
@constraint(uc_miqp, sum(isOn[kk, (jj:nUnits:(nZones*nUnits))]) <= power[kk, jj])
end
end
=#
#Enforce startup=1 when moving from on to off for the first hour
#@constraint(uc_miqp, [jj = 1:nUnits], isOn_ll[1, jj] + isOn_hl[1, jj] - startup[1, jj] <= 0)
for jj = 1:nUnits
@constraint(uc_miqp, sum(isOn[1, (jj:nUnits:(nZones*nUnits))]) - startup[1, jj] == 0)
end
#Enforce startup=1 when moving from on to off
#@constraint(uc_miqp, [kk = 2:nHours, jj = 1:nUnits], (-isOn_hl[(kk-1), jj] + isOn_hl[kk, jj]) + (-isOn_ll[(kk-1), jj] + isOn_ll[kk, jj]) - startup[kk, jj] <= 0)
for jj = 1:nUnits
for kk = 2:nHours
@constraint(uc_miqp, (sum(isOn[kk, (jj:nUnits:(nZones*nUnits))]) - sum(isOn[(kk-1), (jj:nUnits:(nZones*nUnits))])) - startup[kk, jj] <= 0)
@constraint(uc_miqp, (1 - sum(isOn[(kk-1), (jj:nUnits:(nZones*nUnits))])) - startup[kk, jj] >= 0)
end
end
#Enforce shutdown = 0 for the first hour
for jj = 1:nUnits
@constraint(uc_miqp, shutdown[1, jj] <= 0)
end
#Enforce shutdown = 1 when moving from on to off
#@constraint(uc_miqp, [kk = 2:nHours, jj = 1:nUnits], (isOn_hl[(kk-1), jj] - isOn_hl[kk, jj]) + (isOn_ll[(kk-1), jj] - isOn_ll[kk, jj]) - shutdown[kk, jj] <= 0)
for jj = 1:nUnits
for kk = 2:nHours
@constraint(uc_miqp, (sum(isOn[(kk-1), (jj:nUnits:(nZones*nUnits))]) - sum(isOn[kk, (jj:nUnits:(nZones*nUnits))])) - shutdown[kk, jj] <= 0)
@constraint(uc_miqp, (sum(isOn[(kk-1), (jj:nUnits:(nZones*nUnits))])) - shutdown[kk, jj] >= 0)
end
end
#Minimum Up Time constraint
#@constraint(uc_miqp, [kk = 1:nHours, jj = 1:nUnits; kk > nHours - MinimumUpTime[jj]], startup[kk, jj] - (sum(isOn_ll[kk:nHours, jj]) + sum(isOn_hl[kk:nHours, jj]))/length(kk:nHours) <= 0)
#@constraint(uc_miqp, [kk = 1:nHours, jj = 1:nUnits; kk <= nHours - MinimumUpTime[jj]], startup[kk, jj] - (sum(isOn_ll[kk:(kk+MinimumUpTime[jj]-1), jj]) + sum(isOn_hl[kk:(kk+MinimumUpTime[jj]-1), jj]))/length(kk:(kk+MinimumUpTime[jj]-1)) <= 0)
for jj = 1:nUnits
if StartUpTime[jj] > 0
for kk = 1:nHours
if kk <= 1
if kk > nHours - MinimumUpTime[jj]
sumidx = collect(kk:nHours)
else
sumidx = collect(kk:(kk+MinimumUpTime[jj]-1))
end
@constraint(uc_miqp, startup[kk, jj] - (sum(isOn[sumidx, (jj:nUnits:(nZones*nUnits))]))/length(sumidx) <= 0)
else
if kk <= StartUpTime[jj]
sumidx2 = collect(1:(kk-1))
if kk > nHours - MinimumUpTime[jj]
sumidx = collect(kk:nHours)
else
sumidx = collect(kk:kk+MinimumUpTime[jj]-1)
end
else
sumidx2 = collect((kk-StartUpTime[jj]):(kk-1))
if kk > nHours - MinimumUpTime[jj]
sumidx = collect(kk:nHours)
else
sumidx = collect(kk:(kk+MinimumUpTime[jj]-1))
end
end
@constraint(uc_miqp, startup[kk, jj] - (sum(isOn[sumidx, (jj:nUnits:(nZones*nUnits))]))/length(sumidx) <= 0)
@constraint(uc_miqp, startup[kk, jj] + (sum(isOn[sumidx2, (jj:nUnits:(nZones*nUnits))]))/length(sumidx2) <= 1)
end
end
else
for kk = 1:nHours
if kk > nHours - MinimumUpTime[jj]
sumidx = collect(kk:nHours)
else
sumidx = collect(kk:(kk+MinimumUpTime[jj]-1))
end
@constraint(uc_miqp, startup[kk, jj] - (sum(isOn[sumidx, (jj:nUnits:(nZones*nUnits))]))/length(sumidx) <= 0)
end
end
end
#Minimum Down Time constraint
#@constraint(uc_miqp, [kk = 2:nHours, jj = 1:nUnits; kk > nHours - MinimumDownTime[jj]], shutdown[kk, jj] + (sum(isOn_ll[kk:nHours, jj]) + sum(isOn_hl[kk:nHours, jj]))/length(kk:nHours) <= 1)
#@constraint(uc_miqp, [kk = 2:nHours, jj = 1:nUnits; kk <= nHours - MinimumDownTime[jj]], shutdown[kk, jj] + (sum(isOn_ll[kk:kk+MinimumDownTime[jj]-1, jj]) + sum(isOn_hl[kk:kk+MinimumDownTime[jj]-1, jj]))/length(kk:kk+MinimumDownTime[jj]-1) <= 1)
for jj = 1:nUnits
for kk = 2:nHours
if kk > nHours - MinimumDownTime[jj]
sumidx = collect(kk:nHours)
else
sumidx = collect(kk:(kk+MinimumDownTime[jj]-1))
end
@constraint(uc_miqp, shutdown[kk, jj] + (sum(isOn[sumidx, (jj:nUnits:(nZones*nUnits))]))/length(sumidx) <= 1)
end
end
#Ramp Up Rate constraint
#@constraint(uc_miqp, [kk = 2:nHours, jj = 1:nUnits], -power[(kk-1), jj] + power[kk, jj] <= RampUp[jj] + (maximum([(minPow[jj] - RampUp[jj]), 0]))*startup[kk, jj])
for jj = 1:nUnits
for kk = 2:nHours
@constraint(uc_miqp, -power[(kk-1), jj] + power[kk, jj] <= RampUp[jj] + (maximum([(minPow[jj] - RampUp[jj]), 0]))*startup[kk, jj])
end
end
#Ramp Down Rate constraint
#@constraint(uc_miqp, [kk = 2:nHours, jj = 1:nUnits], power[(kk-1), jj] - power[kk, jj] <= maximum([minPow[jj], RampDown[jj]]) - (maximum([(minPow[jj] - RampDown[jj]), 0]))*isOn_ll[kk, jj] - (maximum([(minPow[jj] - RampDown[jj]), 0]))*isOn_hl[kk, jj])
for jj = 1:nUnits
for kk = 2:nHours
@constraint(uc_miqp, power[(kk-1), jj] - power[kk, jj] <= maximum([minPow[jj], RampDown[jj]]) - (maximum([(minPow[jj] - RampDown[jj]), 0]))*sum(isOn[kk, (jj:nUnits:(nZones*nUnits))]))
end
end
#Ensure that Unit is off for at least MinimumDownTime before Real Start
for jj = 1:nUnits
for kk = 2:nHours
if kk <= MinimumDownTime[jj]
sumidx = collect(1:(kk-1))
else
sumidx = collect((kk-MinimumDownTime[jj]):(kk-1))
end
@constraint(uc_miqp, rstartup2[kk, jj] + (sum(isOn[sumidx, (jj:nUnits:(nZones*nUnits))]))/length(sumidx) <= 1)
end
end
#Real Startup constraint
#@constraint(uc_miqp, [kk = 1:nHours, jj = 1:nUnits; kk <= StartUpTime[jj]], rstartup[kk, jj] == startup[kk, jj])
#@constraint(uc_miqp, [kk = 1:nHours, jj = 1:nUnits; kk > StartUpTime[jj]], rstartup[(kk - StartUpTime[jj]), jj] == startup[kk, jj])
for jj = 1:nUnits
for kk = 1:nHours
if kk <= StartUpTime[jj]
#@constraint(uc_miqp, rstartup1[kk, jj] >= startup[kk, jj])
#@constraint(uc_miqp, rstartup1[kk, jj] <= startup[kk, jj])
@constraint(uc_miqp, rstartup1[kk, jj] == startup[kk, jj])
else
#@constraint(uc_miqp, rstartup1[kk, jj] >= 0)
#@constraint(uc_miqp, rstartup1[kk, jj] <= 0)
@constraint(uc_miqp, rstartup1[kk, jj] == 0)
#@constraint(uc_miqp, rstartup2[(kk - StartUpTime[jj]), jj] >= startup[kk, jj])
#@constraint(uc_miqp, rstartup2[(kk - StartUpTime[jj]), jj] <= startup[kk, jj])
@constraint(uc_miqp, rstartup2[(kk - StartUpTime[jj]), jj] == startup[kk, jj])
end
end
end
#SOLVE IT AND DISPLAY THE RESULTS
#--------------------------------
status = optimize!(uc_miqp) # solves the model
#Save Solution Files
isOn_Zone = value.(isOn)
isOn_Gen = zeros(nHours, nUnits)
for tt = 1:nHours
for jj = 1:nUnits
isOn_Gen[tt, jj] = sum(isOn_Zone[tt, (jj:nUnits:(nZones*nUnits))])
end
end
CSV.write("Option_$opt\\Group_$gp\\$var\\powerprofile.csv", DataFrame(value.(power)), writeheader=false)
CSV.write("Option_$opt\\Group_$gp\\$var\\commitmentprofile.csv", DataFrame(value.(isOn)), writeheader=false)
CSV.write("Option_$opt\\Group_$gp\\$var\\commitmentprofile_gen.csv", DataFrame(isOn_Gen), writeheader=false)
CSV.write("Option_$opt\\Group_$gp\\$var\\startupprofile.csv", DataFrame(value.(startup)), writeheader=false)
CSV.write("Option_$opt\\Group_$gp\\$var\\realstartupprofile1.csv", DataFrame(value.(rstartup1)), writeheader=false)
CSV.write("Option_$opt\\Group_$gp\\$var\\realstartupprofile2.csv", DataFrame(value.(rstartup2)), writeheader=false)
CSV.write("Option_$opt\\Group_$gp\\$var\\shutdownprofile.csv", DataFrame(value.(shutdown)), writeheader=false)
FuelCost_sol = zeros(nHours, nUnits)
StartUpCost_sol1 = zeros(nHours, nUnits)
StartUpCost_sol2 = zeros(nHours, nUnits)
ShutDownCost_sol = zeros(nHours, nUnits)
for jj = 1:nUnits
for kk = 1:nHours
FuelCost_sol[kk, jj] = a[jj]*(value.(power[kk, jj]^2)) + b[jj]*(value.(power[kk, jj])) + c[jj]*(sum(value.(isOn[kk, (jj:nUnits:(nZones*nUnits))])))
StartUpCost_sol1[kk, jj] = StartUpCost[jj]*(value.(rstartup1[kk, jj]))
StartUpCost_sol2[kk, jj] = StartUpCost[jj]*(value.(rstartup2[kk, jj]))
ShutDownCost_sol[kk, jj] = ShutDownCost[jj]*(value.(shutdown[kk, jj]))
end
end
CSV.write("Option_$opt\\Group_$gp\\$var\\FuelCost.csv", DataFrame(FuelCost_sol), writeheader=false)
CSV.write("Option_$opt\\Group_$gp\\$var\\StartUpCost1.csv", DataFrame(StartUpCost_sol1), writeheader=false)
CSV.write("Option_$opt\\Group_$gp\\$var\\StartUpCost2.csv", DataFrame(StartUpCost_sol2), writeheader=false)
CSV.write("Option_$opt\\Group_$gp\\$var\\ShutDownCost.csv", DataFrame(ShutDownCost_sol), writeheader=false)
#Return Objective Value
Obj_value = objective_value(uc_miqp)
ans = "Obj Value = $Obj_value"
open("Option_$opt\\Group_$gp\\$var\\Obj_Value.txt", "w") do file
write(file, ans)
end
return objective_value(uc_miqp)
end
#SET NUMBER OF ITERATIONS
#------------------------
num_iter = 10
#SET NUMBER OF GROUPS
#--------------------
num_gp = 9
#SET GENERATORS FOR EACH GROUP
#-----------------------------
gen_gp = [[1, 2], [3, 4], [5, 6, 7, 8], [9, 10, 11, 12, 13, 14], [15, 16, 17], [18, 19, 20], [21], [22, 23], [24, 25, 26, 27]]
#CREATE VARIABLES TO STORE RESULTS
#---------------------------------
test1 = zeros(num_iter, (num_gp*2))
test2 = zeros(num_iter, (num_gp*2))
test3 = zeros(num_iter, (num_gp*2))
test4 = zeros(num_iter, (num_gp*2))
test5 = zeros(num_iter, (num_gp*2))
test6 = zeros(num_iter, (num_gp*2))
test7 = zeros(num_iter, (num_gp*2))
test8 = zeros(num_iter, (num_gp*2))
test9 = zeros(num_iter, (num_gp*2))
#SYSTEM DATA
#--------------------
nUnits = 27
nHours = 168
data = DataFrame(CSV.read("GenData.csv"))
a = reshape((data[1:nUnits, :a]), (1, nUnits))
b = reshape((data[1:nUnits, :b]), (1, nUnits))
c = convert(Array{Float64}, reshape((data[1:nUnits, :c]), (1, nUnits)))
maxPow = reshape((data[1:nUnits, :maxPow]), (1, nUnits))
minPow = reshape((data[1:nUnits, :minPow]), (1, nUnits))
StartUpTime = reshape((data[1:nUnits, :StartUpTime]), (1, nUnits))
MinimumUpTime = reshape((data[1:nUnits, :MinimumUpTime]), (1, nUnits))
MinimumDownTime = reshape((data[1:nUnits, :MinimumDownTime]), (1, nUnits))
RampUp = reshape((data[1:nUnits, :RampUp]), (1, nUnits))
RampDown = reshape((data[1:nUnits, :RampDown]), (1, nUnits))
StartUpCost = convert(Array{Float64}, reshape((data[1:nUnits, :StartUpCost]), (1, nUnits)))
ZeroStartUpCost = reshape((data[1:nUnits, :ZeroStartUpCost]), (1, nUnits))
ShutDownCost = convert(Array{Float64}, reshape((data[1:nUnits, :ShutDownCost]), (1, nUnits)))
numZones = reshape((data[1:nUnits, :numZones]), (1, nUnits))
maxnumZones = maximum(numZones)
nZones = maxnumZones
DLE = zeros(Int8, nUnits, maxnumZones)
for jj = 1:nUnits
for zz = 1:maxnumZones
if zz <= numZones[jj]
DLE[jj, zz] = 1
else
DLE[jj, zz] = 0
end
end
end
DLE_data1 = DataFrame(CSV.read("DLEZone.csv", header=false))
DLEZone = convert(Matrix, DLE_data1)
load_srdata = DataFrame(CSV.read("Load_SRData.csv"))
Load = reshape((load_srdata[1:nHours, :Load]), (1, nHours))
SR = reshape((load_srdata[1:nHours, :SR]), (1, nHours))
#GET INITIAL SOLUTION
#--------------------
init_power = convert(Matrix, DataFrame(CSV.read("init_powerprofile.csv", header=false)))
init_isOn = convert(Matrix, DataFrame(CSV.read("init_commitmentprofile.csv", header=false)))
init_startup = convert(Matrix, DataFrame(CSV.read("init_startupprofile.csv", header=false)))
init_rstartup1 = convert(Matrix, DataFrame(CSV.read("init_realstartupprofile1.csv", header=false)))
init_rstartup2 = convert(Matrix, DataFrame(CSV.read("init_realstartupprofile2.csv", header=false)))
init_shutdown = convert(Matrix, DataFrame(CSV.read("init_shutdownprofile.csv", header=false)))
#SET LIMITS FOR EACH TEST
test1_limit = zeros(1, (2*num_gp))
test2_limit = zeros(1, (2*num_gp))
test3_limit = zeros(1, (2*num_gp))
test4_limit = zeros(1, (2*num_gp))
test5_limit = zeros(1, (2*num_gp))
test6_limit = zeros(1, (2*num_gp))
test7_limit = zeros(1, (2*num_gp))
test8_limit = zeros(1, (2*num_gp))
test9_limit = zeros(1, (2*num_gp))
for gp = 1:num_gp
test1_limit[((gp-1)*2 + 1)] = 1
test1_limit[((gp-1)*2 + 2)] = 10
test2_limit[((gp-1)*2 + 1)] = 0.8*StartUpCost[gen_gp[gp][1]]
test2_limit[((gp-1)*2 + 2)] = 1.2*StartUpCost[gen_gp[gp][1]]
test3_limit[((gp-1)*2 + 1)] = 0.8*ShutDownCost[gen_gp[gp][1]]
test3_limit[((gp-1)*2 + 2)] = 1.2*ShutDownCost[gen_gp[gp][1]]
test4_limit[((gp-1)*2 + 1)] = 0.8*RampUp[gen_gp[gp][1]]
test4_limit[((gp-1)*2 + 2)] = 1.2*RampUp[gen_gp[gp][1]]
test5_limit[((gp-1)*2 + 1)] = 0.1*RampUp[gen_gp[gp][1]]
test5_limit[((gp-1)*2 + 2)] = RampUp[gen_gp[gp][1]]
test6_limit[((gp-1)*2 + 1)] = 0
test6_limit[((gp-1)*2 + 2)] = 9
test7_limit[((gp-1)*2 + 1)] = 0.8*a[gen_gp[gp][1]]
test7_limit[((gp-1)*2 + 2)] = 1.2*a[gen_gp[gp][1]]
test8_limit[((gp-1)*2 + 1)] = 0.8*b[gen_gp[gp][1]]
test8_limit[((gp-1)*2 + 2)] = 1.2*b[gen_gp[gp][1]]
test9_limit[((gp-1)*2 + 1)] = 0.8*c[gen_gp[gp][1]]
test9_limit[((gp-1)*2 + 2)] = 1.2*c[gen_gp[gp][1]]
end
#RUN SIMULATION BASED ON OPTION
#------------------------------
#MIN UP AND DOWN TIME ANALYSIS
if option == 1
mkdir("Option_1")
for gp = 1:num_gp
mkdir("Option_1\\Group_$gp")
global temp_MinimumUpTime = zeros(1, nUnits)
global test_MinimumUpTime = zeros(1, num_iter)
global temp_MinimumDownTime = zeros(1, nUnits)
global test_MinimumDownTime = zeros(1, num_iter)
test_MinimumUpTime = collect(range(test1_limit[((gp-1)*2 + 1)], test1_limit[((gp-1)*2 + 2)], length = num_iter))
test_MinimumDownTime = collect(range(test1_limit[((gp-1)*2 + 1)], test1_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test1[iter, ((gp-1)*2 + 1)] = round(test_MinimumUpTime[iter])
var = round(test_MinimumUpTime[iter])
mkdir("Option_1\\Group_$gp\\$var")
temp_MinimumUpTime = deepcopy(MinimumUpTime)
temp_MinimumDownTime = deepcopy(MinimumDownTime)
for jj = gen_gp[gp]
temp_MinimumUpTime[jj] = round(test_MinimumUpTime[iter])
temp_MinimumDownTime[jj] = round(test_MinimumDownTime[iter])
end
#Run Simulation
global test1[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, b, c, maxPow, minPow, StartUpTime, temp_MinimumUpTime, temp_MinimumDownTime, RampUp, RampDown, StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end
CSV.write("Test1_Result.csv", DataFrame(test1), writeheader=false)
#STARTUP COST ANALYSIS
elseif option == 2
mkdir("Option_2")
for gp = 1:num_gp
mkdir("Option_2\\Group_$gp")
global temp_StartUpCost = zeros(1, nUnits)
global test_StartUpCost = zeros(1, num_iter)
test_StartUpCost = collect(range(test2_limit[((gp-1)*2 + 1)], test2_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test2[iter, ((gp-1)*2 + 1)] = round(test_StartUpCost[iter], digits=2)
var = round(test_StartUpCost[iter], digits=2)
mkdir("Option_2\\Group_$gp\\$var")
temp_StartUpCost = deepcopy(StartUpCost)
for jj = gen_gp[gp]
temp_StartUpCost[jj] = round(test_StartUpCost[iter], digits=2)
end
#Run Simulation
global test2[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, b, c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, temp_StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end
CSV.write("Test2_Result.csv", DataFrame(test2), writeheader=false)
#SHUTDOWN COST ANALYSIS
elseif option == 3
mkdir("Option_3")
for gp = 1:num_gp
mkdir("Option_3\\Group_$gp")
global temp_ShutDownCost = zeros(1, nUnits)
global test_ShutDownCost = zeros(1, num_iter)
test_ShutDownCost = collect(range(test3_limit[((gp-1)*2 + 1)], test3_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test3[iter, ((gp-1)*2 + 1)] = round(test_ShutDownCost[iter], digits=2)
var = round(test_ShutDownCost[iter], digits=2)
mkdir("Option_3\\Group_$gp\\$var")
temp_ShutDownCost = deepcopy(ShutDownCost)
for jj = gen_gp[gp]
temp_ShutDownCost[jj] = round(test_ShutDownCost[iter], digits=2)
end
#Run Simulation
global test3[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, b, c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, StartUpCost, ZeroStartUpCost, temp_ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end
CSV.write("Test3_Result.csv", DataFrame(test3), writeheader=false)
#RAMP RATE ANALYSIS
elseif option == 4
mkdir("Option_4")
for gp = 1:num_gp
mkdir("Option_4\\Group_$gp")
global temp_RampUp = zeros(1, nUnits)
global test_RampUp = zeros(1, num_iter)
global temp_RampDown = zeros(1, nUnits)
global test_RampDown = zeros(1, num_iter)
test_RampUp = collect(range(test4_limit[((gp-1)*2 + 1)], test4_limit[((gp-1)*2 + 2)], length = num_iter))
test_RampDown = collect(range(test4_limit[((gp-1)*2 + 1)], test4_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test4[iter, ((gp-1)*2 + 1)] = round(test_RampUp[iter], digits=3)
var = round(test_RampUp[iter], digits=3)
mkdir("Option_4\\Group_$gp\\$var")
temp_RampUp = deepcopy(RampUp)
temp_RampDown = deepcopy(RampDown)
for jj = gen_gp[gp]
temp_RampUp[jj] = round(test_RampUp[iter], digits=3)
temp_RampDown[jj] = round(test_RampDown[iter], digits=3)
end
#Run Simulation
global test4[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, b, c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, temp_RampUp, temp_RampDown, StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end
CSV.write("Test4_Result.csv", DataFrame(test4), writeheader=false)
#OPERATIONAL ZONE ANALYSIS
elseif option == 5
mkdir("Option_5")
for gp = 1:num_gp
mkdir("Option_5\\Group_$gp")
global temp_DLEZone = zeros(1, nUnits)
global test_DLEgap = zeros(1, num_iter)
test_DLEgap = collect(range(test5_limit[((gp-1)*2 + 1)], test5_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test5[iter, ((gp-1)*2 + 1)] = round(test_DLEgap[iter], digits=3)
var = round(test_DLEgap[iter], digits=3)
mkdir("Option_5\\Group_$gp\\$var")
temp_DLEZone = deepcopy(DLEZone)
for jj = gen_gp[gp]
temp_DLEZone[jj, 2] = DLEZone[jj, 2] + ((DLEZone[jj, 3] - DLEZone[jj, 2])/2) - (round(test_DLEgap[iter], digits=3)/2)
temp_DLEZone[jj, 3] = DLEZone[jj, 3] - ((DLEZone[jj, 3] - DLEZone[jj, 2])/2) + (round(test_DLEgap[iter], digits=3)/2)
end
#Run Simulation
global test5[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, b, c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, temp_DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end
CSV.write("Test5_Result.csv", DataFrame(test5), writeheader=false)
#STARTUP TIME ANALYSIS
elseif option == 6
mkdir("Option_6")
for gp = 1:num_gp
mkdir("Option_6\\Group_$gp")
global temp_StartUpTime = zeros(1, nUnits)
global test_StartUpTime = zeros(1, num_iter)
test_StartUpTime = collect(range(test6_limit[((gp-1)*2 + 1)], test6_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test6[iter, ((gp-1)*2 + 1)] = round(test_StartUpTime[iter])
var = round(test_StartUpTime[iter])
mkdir("Option_6\\Group_$gp\\$var")
temp_StartUpTime = deepcopy(StartUpTime)
for jj = gen_gp[gp]
temp_StartUpTime[jj] = round(test_StartUpTime[iter])
end
#Run Simulation
global test6[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, b, c, maxPow, minPow, temp_StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end
CSV.write("Test6_Result.csv", DataFrame(test6), writeheader=false)
#GAMMA ANALYSIS
elseif option == 7
mkdir("Option_7")
for gp = 1:num_gp
mkdir("Option_7\\Group_$gp")
global temp_a = zeros(1, nUnits)
global test_a = zeros(1, num_iter)
test_a = collect(range(test7_limit[((gp-1)*2 + 1)], test7_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test7[iter, ((gp-1)*2 + 1)] = test_a[iter]
var = test_a[iter]
mkdir("Option_7\\Group_$gp\\$var")
temp_a = deepcopy(a)
for jj = gen_gp[gp]
temp_a[jj] = test_a[iter]
end
#Run Simulation
global test7[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, temp_a, b, c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end
CSV.write("Test7_Result.csv", DataFrame(test7), writeheader=false)
#BETA ANALYSIS
elseif option == 8
mkdir("Option_8")
for gp = 1:num_gp
mkdir("Option_8\\Group_$gp")
global temp_b = zeros(1, nUnits)
global test_b = zeros(1, num_iter)
test_b = collect(range(test8_limit[((gp-1)*2 + 1)], test8_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test8[iter, ((gp-1)*2 + 1)] = test_b[iter]
var = test_b[iter]
mkdir("Option_8\\Group_$gp\\$var")
temp_b = deepcopy(b)
for jj = gen_gp[gp]
temp_b[jj] = test_b[iter]
end
#Run Simulation
global test8[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, temp_b, c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end
CSV.write("Test8_Result.csv", DataFrame(test8), writeheader=false)
#ALPHA ANALYSIS
elseif option == 9
mkdir("Option_9")
for gp = 1:num_gp
mkdir("Option_9\\Group_$gp")
global temp_c = zeros(1, nUnits)
global test_c = zeros(1, num_iter)
test_c = collect(range(test9_limit[((gp-1)*2 + 1)], test9_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test9[iter, ((gp-1)*2 + 1)] = round(test_c[iter], digits=2)
var = round(test_c[iter], digits=2)
mkdir("Option_9\\Group_$gp\\$var")
temp_c = deepcopy(c)
for jj = gen_gp[gp]
temp_c[jj] = round(test_c[iter], digits=2)
end
#Run Simulation
global test9[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, b, temp_c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end
CSV.write("Test9_Result.csv", DataFrame(test9), writeheader=false)
else
println("Choose an valid option!")
end
Yes. Just put all this inside a function, which takes in the params that are global inputs now. Then you can run it as many times as needed.
2 Likes
How it that different from what I’m doing above. Is there something wrong with having the function in the same script as the for loop?
Without copying and pasting your entire script, could you clarify at what “level” you want to run things with different inputs every time?
Whenever I call the MIQP function in the for loop. That line usually looks like the following:
#Run Simulation
global test2[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, b, c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, temp_StartUpCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
Then perhaps I dont understand your problem, unfortunately.
I think your problem has to do with over-using globals? Are they not updating in the loop appropriately?
an MWE, something I can copy and paste which isolates your problem, would be useful.
1 Like
The values are updating in the loop appropriately (I have checked this already), the problem seems to lies with the solution obtained when the MIQP function is called from one iteration of the loop to the next. I am getting a correct solution (I said “a” because there are multiple solutions which can work for this optimization problem), however, there are instances where the solution is completely different from the solution for the next iteration in the loop. While both are correct, for the sensitivity analysis study, I know for this case that the change in coefficient of the objective value (this has to do with sensitivity analysis) should not affect the solution significantly. However, there are instances where I am seeing major changes when I should not see any. So my initial thoughts were that I should clear the workspace so that the solver is not influenced by the solution produced from the last iteration (in the way an initial solution could may affect outcome achieved by the solver. This is why I manually ran the simulation for the irregular case by updating the values themselves and running them in separate instances of julia, hence my initial question of how to run a new julia session. When running in separate instances the solution resembles what it should have been like.
If you start Julia many times, you will notice that everything becomes slow because things need to be compiled every time. Also, you cannot easily pass information back and forth. Your script contains 40 globals, my bet is that they can all be avoided when passing the program’s state around a bit more careful. That is, move the logic into functions and give the functions name. For example, instead of
#STARTUP COST ANALYSIS
elseif option == 2
mkdir("Option_2")
for gp = 1:num_gp
mkdir("Option_2\\Group_$gp")
global temp_StartUpCost = zeros(1, nUnits)
global test_StartUpCost = zeros(1, num_iter)
test_StartUpCost = collect(range(test2_limit[((gp-1)*2 + 1)], test2_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test2[iter, ((gp-1)*2 + 1)] = round(test_StartUpCost[iter], digits=2)
var = round(test_StartUpCost[iter], digits=2)
mkdir("Option_2\\Group_$gp\\$var")
temp_StartUpCost = deepcopy(StartUpCost)
for jj = gen_gp[gp]
temp_StartUpCost[jj] = round(test_StartUpCost[iter], digits=2)
end
#Run Simulation
global test2[iter, ((gp-1)*2 + 2)] = MIQP(nUnits, nHours, a, b, c, maxPow, minPow, StartUpTime, MinimumUpTime, MinimumDownTime, RampUp, RampDown, temp_StartU pCost, ZeroStartUpCost, ShutDownCost, nZones, DLE, DLEZone, Load, SR, init_power, init_isOn, init_startup, init_rstartup1, init_rstartup2, init_shutdown, opt, gp, var)
end
end---
CSV.write("Test2_Result.csv", DataFrame(test2), writeheader=false)--
#SHUTDOWN COST ANALYSIS
elseif option == 3
write a function
function analyze_startup_costs(...)
[...]
end
and call that from the code via
elseif option == 2
analyze_startup_costs(...)
elseif option == 3
The next thing that you will walk into is that the functions receive too many variables, like happens in your MIQP function. To solve that, store you variables in structs or other appropriate data types. For instance,
struct ModelSettings
RampUp
RampDown
StartupCost
end
Implementing functions will take some time, but I guarantee that it will make it much easier to manage the complexity of the program.
4 Likes
I stripped down the code to include only the necessary sections dealing with the logic and loop statements at the end (i.e. the section where I called my MIQP function). Can you let me know if I am on the right path wrt to the simplification of the code:
using JuMP
#using CPLEX
using Gurobi
using CSV
using DataFrames
using XLSX
function MIQP(StartUpCost, var)
return var
end
function analyze_startup_costs(num_gp, num_iter, gen_gp, StartUpCost)
#SET LIMITS FOR EACH TEST
test2_limit = zeros(1, (2*num_gp))
for gp = 1:num_gp
test2_limit[((gp-1)*2 + 1)] = 0.8*StartUpCost[gen_gp[gp][1]]
test2_limit[((gp-1)*2 + 2)] = 1.2*StartUpCost[gen_gp[gp][1]]
end
#CREATE VARIABLES TO STORE RESULTS
#---------------------------------
test2 = zeros(num_iter, (num_gp*2))
for gp = 1:num_gp
test_StartUpCost = collect(range(test2_limit[((gp-1)*2 + 1)], test2_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
global test2[iter, ((gp-1)*2 + 1)] = round(test_StartUpCost[iter], digits=2)
var = round(test_StartUpCost[iter], digits=2)
temp_StartUpCost = StartUpCost
for jj = gen_gp[gp]
temp_StartUpCost[jj] = round(test_StartUpCost[iter], digits=2)
end
#Run Simulation
global test2[iter, ((gp-1)*2 + 2)] = MIQP(temp_StartUpCost, var)
end
end
println(test2)
end
#SET NUMBER OF ITERATIONS
#------------------------
const num_iter = 10
#SET NUMBER OF GROUPS
#--------------------
const num_gp = 9
#SET GENERATORS FOR EACH GROUP
#-----------------------------
const gen_gp = [[1, 2], [3, 4], [5, 6, 7, 8], [9, 10, 11, 12, 13, 14], [15, 16, 17], [18, 19, 20], [21], [22, 23], [24, 25, 26, 27]]
#SYSTEM DATA
#--------------------
const nUnits = 27
const nHours = 168
const data = DataFrame(CSV.read("GenData.csv"))
const StartUpCost = convert(Array{Float64}, reshape((data[1:nUnits, :StartUpCost]), (1, nUnits)))
analyze_startup_costs(num_gp, num_iter, gen_gp, StartUpCost)
I still don’t understand your problem fully. In particular, the line
global test2[iter, ((gp-1)*2 + 1)] = round(test_StartUpCost[iter], digits=2)
Why not just pass the test2 into the function and update it inside the function? Why does it need to be global?
1 Like
Test 2 is being used to store the solution values for each iteration of each group.
The following is used to store the coefficient value that will change in the for loop:
global test2[iter, ((gp-1)*2 + 1)]
While
global test2[iter, ((gp-1)*2 + 2)]
will store the corresponding solution obtained after calling the MIQP function.
I didn’t want to update in the function because I want one variable (test2) to store the solution for all iterations for all groups i.e. in the for loops.
Just to reiterate my problem, the code while it may not be the best, works and gives me solutions. However there are instances where the solution I get does not match what I know it should be. So my initial thoughts are something is affecting the gurobi solver from iteration to iteration. My initial thoughts were that it had to do something with the creation of a new Julia session since if I were to run the simulations manually each with a new Julia session and call the functions myself, the solutions is what I expect. However, by doing it the same in one Julia session is causing irregular solutions.
1 Like
Because arrays are mutable in julia, you can still pass test2 to a function and it will update it, and you will be able to see those changes.
I maintain that the solution is to use an even bigger function. Just do
function main()
# everything in your script now
end
and do not use globals. And see if that fixes the problem.
1 Like
OK will try it out
1 Like
Edit:
This is an x-y problem. You don’t really want to know how to restart julia (which is a fairly nuclear option as far as a Julia script goes), you want to know how to not have persistent state between runs.
(I get sarcastic when I’m tired. A poor excuse, sorry)
Ok so I split my code into 3 separate files: 1. Main.jl, 2. analyze_startup_costs.jl and 3. MIQP.jl.
I got rid of the global declaration and this is what I have now. Just want to know if I am using include correctly or would using it like this cause some confusion:
Main.jl:
using JuMP
#using CPLEX
using Gurobi
using CSV
using DataFrames
using XLSX
include("MIQP.jl")
include("analyze_startup_costs.jl")
#SET NUMBER OF ITERATIONS
#------------------------
const num_iter = 10
#SET NUMBER OF GROUPS
#--------------------
const num_gp = 9
#SET GENERATORS FOR EACH GROUP
#-----------------------------
const gen_gp = [[1, 2], [3, 4], [5, 6, 7, 8], [9, 10, 11, 12, 13, 14], [15, 16, 17], [18, 19, 20], [21], [22, 23], [24, 25, 26, 27]]
#SYSTEM DATA
#--------------------
const nUnits = 27
const nHours = 168
const data = DataFrame(CSV.read("GenData.csv"))
const StartUpCost = convert(Array{Float64}, reshape((data[1:nUnits, :StartUpCost]), (1, nUnits)))
analyze_startup_costs(num_gp, num_iter, gen_gp, StartUpCost)
analyze_startup_costs.jl:
function analyze_startup_costs(num_gp, num_iter, gen_gp, StartUpCost)
#SET LIMITS FOR EACH TEST
test2_limit = zeros(1, (2*num_gp))
for gp = 1:num_gp
test2_limit[((gp-1)*2 + 1)] = 0.8*StartUpCost[gen_gp[gp][1]]
test2_limit[((gp-1)*2 + 2)] = 1.2*StartUpCost[gen_gp[gp][1]]
end
#CREATE VARIABLES TO STORE RESULTS
#---------------------------------
test2 = zeros(num_iter, (num_gp*2))
for gp = 1:num_gp
test_StartUpCost = collect(range(test2_limit[((gp-1)*2 + 1)], test2_limit[((gp-1)*2 + 2)], length = num_iter))
for iter = 1:num_iter
var = round(test_StartUpCost[iter], digits=2)
temp_StartUpCost = StartUpCost
for jj = gen_gp[gp]
temp_StartUpCost[jj] = round(test_StartUpCost[iter], digits=2)
end
#Run Simulation
MIQP(temp_StartUpCost, var, iter, gp, test_StartUpCost, test2)
end
end
println(test2)
end
And MIQP.jl:
function MIQP(temp_StartUpCost, var, iter, gp, test_StartUpCost, test2)
test2[iter, ((gp-1)*2 + 1)] = round(test_StartUpCost[iter], digits=2)
test2[iter, ((gp-1)*2 + 2)] = var
end