updated jl algorithms, added argparses, modified definition of c

config_model.yml

@@ -5,13 +5,11 @@ data_path: 'D:/ULg_PhD_work/datasets/resite_ip/'
 # Spatial resolution (in degrees) of the potential sites.
 spatial_resolution: 0.25
 # Start time and end time of the analysis.
-time_slice: ['2018-01-01T00:00', '2018-01-31T23:00']
+time_slice: ['2011-01-01T00:00', '2020-12-31T23:00']
 # Technologies to deploy.
-#regions: ['GB', 'NL', 'FR', 'DE', 'DK', 'NO', 'PL', 'IE', 'IT', 'SE', 'FI', 'ES', 'GR', 'PT', 'BE', 'LT', 'LV', 'EE', 'HR']
-regions: ['BE', 'NL', 'LT', 'LV', 'EE', 'HR']
+regions: ['GB', 'NL', 'FR', 'DE', 'DK', 'NO', 'PL', 'IE', 'IT', 'SE', 'FI', 'ES', 'GR', 'PT', 'BE', 'LT', 'LV', 'EE', 'HR']
 technologies: ['wind_offshore']
-#deployments: [[80], [60], [57], [36], [35], [30], [28], [22], [20], [20], [15], [13], [10], [9], [6], [4], [3], [1], [1]]
-deployments: [[6], [60], [4], [3], [1], [1]]
+deployments: [[80], [60], [57], [36], [35], [30], [28], [22], [20], [20], [15], [13], [10], [9], [6], [4], [3], [1], [1]]
 
 siting_params:
   # Defines how \alpha is considered in space and time.
@@ -21,23 +19,27 @@ siting_params:
     smoothing: 'mean' # 'median'
     norm: 'min' # 'max'
   # Time-window length used to compute the criticality indicator. Integer value.
-  delta: 3                          # \in \mathbb{N}
+  delta: 1                          # \in \mathbb{N}
   # Threshold
-  c: 1                              # < n, \in \mathbb{N}
+  c: 0.1                            # \in [0, 1]
   # Solution method: BB or HEU or RAND or GRED.
   solution_method:
     SA:
       set: True
       neighborhood: 1
-      no_iterations: 100
-      no_epochs: 100
-      initial_temp: 200.
-      p: 0.05
-      no_runs: 30
-      no_runs_init: 20
-      algorithm: 'MIR' # 'MIR', 'SGH'
+      no_iterations: 2000
+      no_epochs: 500
+      initial_temp: 100.
+      p: 0.1
+      no_runs: 100
+      no_runs_init: 100
+      algorithm: 'SGH' # 'MIR', 'SGH'
     SGH:
       set: False
-      p: 0.05
-      no_runs: 30
-      algorithm: 'SGH'
\ No newline at end of file
+      p: 0.1
+      no_runs: 100
+      algorithm: 'SGH'
+    DGH:
+      set: False
+      no_runs: 6
+      algorithm: 'DGH'
\ No newline at end of file

src/jl/SitingHeuristics.jl

@@ -38,14 +38,14 @@ function main_SA(index_dict, deployment_dict, legacy_index_list,
   elseif init_sol_algorithm == "SGH"
     x_init_alg, LB_init_alg = Array{Float64, 2}(undef, R_init, L), Array{Float64, 1}(undef, R_init)
     for r = 1:R_init
-      if (div(r, 10) > 0) & (mod(r, 10) == 0)
+      if (div(r, 5) > 0) & (mod(r, 5) == 0)
         @info "$(Dates.format(now(), "HH:MM:SS")) Run $(r)/$(R) of $(init_sol_algorithm)"
       end
-      x_init_alg[r, :], LB_init_alg[r] = randomised_threshold_greedy_heuristic_partition(D, c, n_partitions, p,
-                                                                                         index_dict, legacy_index_list)
+      x_init_alg[r, :], LB_init_alg[r] = randomised_greedy_heuristic_partition(D, c, n_partitions, p,
+                                                                               index_dict, legacy_index_list)
     end
     LB_init_alg_best = argmax(LB_init_alg)
-    x_init = x_init_alg[LB_init_alg_best, :]
+    x_init = round.(value.(x_init_alg[LB_init_alg_best, :]))
 
   else
     println("No such algorithm available.")
@@ -54,7 +54,7 @@ function main_SA(index_dict, deployment_dict, legacy_index_list,
 
   println("Initial solution retrieved. Starting local search.")
   for r = 1:R
-    if (div(r, 10) > 0) & (mod(r, 10) == 0)
+    if (div(r, 5) > 0) & (mod(r, 5) == 0)
       @info "$(Dates.format(now(), "HH:MM:SS")) Run $(r)/$(R) of LS"
     end
     x_sol[r, :], LB_sol[r], obj_sol[r, :] = simulated_annealing_local_search_partition(D, c, n_partitions,
@@ -88,11 +88,48 @@ function main_SGH(index_dict, deployment_dict, legacy_index_list,
   if algorithm == "SGH"
     x_sol, LB_sol = Array{Float64, 2}(undef, R, L), Array{Float64, 1}(undef, R)
     for r = 1:R
-      if (div(r, 10) > 0) & (mod(r, 10) == 0)
+      if (div(r, 5) > 0) & (mod(r, 5) == 0)
         @info "$(Dates.format(now(), "HH:MM:SS")) Run $(r)/$(R)"
       end
-      x_sol[r, :], LB_sol[r] = randomised_threshold_greedy_heuristic_partition(D, c, n_partitions, p,
-                                                                               index_dict, legacy_index_list)
+      x_sol[r, :], LB_sol[r] = randomised_greedy_heuristic_partition(D, c, n_partitions, p,
+                                                                     index_dict, legacy_index_list)
+    end
+
+  else
+    println("No such algorithm available.")
+    throw(ArgumentError)
+  end
+
+  return x_sol, LB_sol
+
+end
+
+function main_DGH(index_dict, deployment_dict, legacy_index_list,
+                  criticality_matrix, c,
+                  no_runs=100, algorithm="DGH")
+
+  index_dict = Dict([(convert(Int64, k), convert(Int64, index_dict[k])) for k in keys(index_dict)])
+  deployment_dict = Dict([(convert(Int64, k), convert(Int64, deployment_dict[k])) for k in keys(deployment_dict)])
+  legacy_index_list = Vector{Int64}(vec(legacy_index_list))
+
+  D  = convert.(Float64, criticality_matrix)
+  c = convert(Float64, c)
+  R = convert(Int64, no_runs)
+  algorithm = string(algorithm)
+
+  W, L = size(D)
+
+  P = maximum(values(index_dict))
+  n_partitions = Vector{Int64}(vec([deployment_dict[i] for i in 1:P]))
+
+  if algorithm == "DGH"
+    x_sol, LB_sol = Array{Float64, 2}(undef, R, L), Array{Float64, 1}(undef, R)
+    for r = 1:R
+      if (div(r, 2) > 0) & (mod(r, 2) == 0)
+        @info "$(Dates.format(now(), "HH:MM:SS")) Run $(r)/$(R)"
+      end
+      x_sol[r, :], LB_sol[r] = greedy_heuristic_partition(D, c, n_partitions,
+                                                          index_dict, legacy_index_list)
     end
 
   else

src/main.py

 import yaml
 import julia
 from os.path import join
-from numpy import argmax
+from numpy import argmax, ceil
+import argparse
 
 from helpers import read_inputs, init_folder, xarray_to_ndarray, generate_jl_input, \
     get_potential_per_site, capacity_to_cardinality
@@ -12,14 +13,45 @@ import logging
 logging.basicConfig(level=logging.INFO, format=f"%(levelname)s %(asctime)s - %(message)s", datefmt='%Y-%m-%d %H:%M:%S')
 logger = logging.getLogger(__name__)
 
+def parse_args():
+
+    parser = argparse.ArgumentParser(description='Command line arguments.')
+
+    parser.add_argument('--c', type=float)
+    parser.add_argument('--run_DGH', type=bool, default=False)
+    parser.add_argument('--run_SGH', type=bool, default=False)
+    parser.add_argument('--run_LS', type=bool, default=False)
+    parser.add_argument('--LS_init_algorithm', type=str, default=None)
+
+    parsed_args = vars(parser.parse_args())
+
+    return parsed_args
+
+
+def single_true(iterable):
+    i = iter(iterable)
+    return any(i) and not any(i)
+
+
 if __name__ == '__main__':
 
+    args = parse_args()
+
+    if not single_true([args['run_LS'], args['run_DGH'], args['run_SGH']]):
+        raise ValueError(' More than one run selected in the argparser.')
+
     logger.info('Starting data pre-processing.')
 
     model_parameters = read_inputs('../config_model.yml')
     siting_parameters = model_parameters['siting_params']
     tech_parameters = read_inputs('../config_techs.yml')
 
+    siting_parameters['c'] = args['c']
+    siting_parameters['solution_method']['SA']['set'] = args['run_LS']
+    siting_parameters['solution_method']['SA']['algorithm'] = args['LS_init_algorithm']
+    siting_parameters['solution_method']['DGH']['set'] = args['run_DGH']
+    siting_parameters['solution_method']['SGH']['set'] = args['run_SGH']
+
     data_path = model_parameters['data_path']
     spatial_resolution = model_parameters['spatial_resolution']
     time_horizon = model_parameters['time_slice']
@@ -39,6 +71,8 @@ if __name__ == '__main__':
                                                                  deployment_dict, model_parameters))
 
     jl_dict = generate_jl_input(deployment_dict, site_coordinates, site_positions, legacy_coordinates)
+    c = int(ceil(siting_parameters['c'] * sum(deployment_dict[r][t] for r in deployment_dict.keys()
+                                                                    for t in deployment_dict[r].keys())))
 
     logger.info('Data pre-processing finished. Opening Julia instance.')
     j = julia.Julia(compiled_modules=False)
@@ -53,26 +87,39 @@ if __name__ == '__main__':
         jl_sel, jl_obj, _ = Main.main_SA(jl_dict['index_dict'],
                                          jl_dict['deployment_dict'],
                                          jl_dict['legacy_site_list'],
-                                         criticality_data, siting_parameters['c'],
+                                         criticality_data, c,
                                          params['neighborhood'], params['initial_temp'], params['p'],
                                          params['no_iterations'], params['no_epochs'],
                                          params['no_runs'], params['no_runs_init'],
                                          params['algorithm'])
 
-        output_folder = init_folder(model_parameters, siting_parameters['c'], suffix=f"_SA_{params['algorithm']}")
+        output_folder = init_folder(model_parameters, c, suffix=f"_SA_{params['algorithm']}")
 
     elif siting_parameters['solution_method']['SGH']['set']:
 
         params = siting_parameters['solution_method']['SGH']
-        logger.info('SGH chosen to solve the IP.')
+        logger.info('Stochastic greedy heuristic chosen to solve the IP.')
 
         jl_sel, jl_obj = Main.main_SGH(jl_dict['index_dict'],
                                        jl_dict['deployment_dict'],
                                        jl_dict['legacy_site_list'],
-                                       criticality_data, siting_parameters['c'],
+                                       criticality_data, c,
                                        params['p'], params['no_runs'], params['algorithm'])
 
-        output_folder = init_folder(model_parameters, siting_parameters['c'], suffix="_SGH")
+        output_folder = init_folder(model_parameters, c, suffix="_SGH")
+
+    elif siting_parameters['solution_method']['DGH']['set']:
+
+        params = siting_parameters['solution_method']['DGH']
+        logger.info('Deterministic greedy heuristic chosen to solve the IP.')
+
+        jl_sel, jl_obj = Main.main_DGH(jl_dict['index_dict'],
+                                       jl_dict['deployment_dict'],
+                                       jl_dict['legacy_site_list'],
+                                       criticality_data, c,
+                                       params['no_runs'], params['algorithm'])
+
+        output_folder = init_folder(model_parameters, c, suffix="_DGH")
 
     else:
         raise ValueError(' This solution method is not available.')

src/tools.py

@@ -11,7 +11,7 @@ import xarray as xr
 import xarray.ufuncs as xu
 from geopandas import read_file
 from numpy import arange, interp, float32, datetime64, sqrt, asarray, newaxis, sum, max, unique, \
-    radians, cos, sin, arctan2, zeros
+    radians, cos, sin, arctan2, zeros, ceil
 from pandas import read_csv, Series, DataFrame, date_range, concat, MultiIndex, to_datetime
 from shapely.geometry import Point
 from shapely.ops import nearest_points
@@ -715,7 +715,8 @@ def retrieve_index_dict(deployment_vector, coordinate_dict):
 def retrieve_site_data(model_parameters, capacity_factor_data, criticality_data, deployment_dict,
                        location_mapping, comp_site_coordinates, legacy_sites, output_folder, benchmark):
 
-    c = model_parameters['siting_params']['c']
+    c = int(ceil(model_parameters['siting_params']['c'] * sum(deployment_dict[r][t] for r in deployment_dict.keys()
+                                                              for t in deployment_dict[r].keys())))
 
     output_by_tech = collapse_dict_region_level(capacity_factor_data)
     time_dt = date_range(start=model_parameters['time_slice'][0], end=model_parameters['time_slice'][1], freq='H')