Develop

readme.md

@@ -10,7 +10,7 @@ git clone https://github.com/UVAdMIST/swmm_mpc.git
 pip install swmm_mpc/
 ```
 ## 2. Install pyswmm
-swmm\_mpc requires a special version of OWA's pyswmm. To install this do:
+swmm\_mpc requires a special version of OWA's pyswmm. This version of pyswmm has an additional feature that allows saving a hotstart file at any time step in the simulation run. To install this version of pyswmm do:
 
 ```
 pip install git+https://github.com/UVAdMIST/pyswmm.git@feature_save_hotstart
@@ -46,51 +46,85 @@ export PATH="/path/to/swmm5/src:$PATH"
 
 # Usage
 The `run_swmm_mpc` function is the main function (maybe the only function) that 
-you will likely use. Here is an example of how I use it. First I write a small 
+you will likely use. Here is an example of how it is used. `run_swmm_mpc` takes one and only one argument: the path to your configuration file (see example below).  First I write a small 
 script to indicate the model and the different parameters I want to use, like 
 this (called `my_swmm_mpc.py` in the next step):
 
-```python
-from swmm_mpc.swmm_mpc import run_swmm_mpc
-
-inp_file = "/path/to/my/model.inp"
-control_horizon = 1. #hr
-control_time_step = 900. #sec
-control_str_ids = ["ORIFICE R1", "ORIFICE R2"]
-results_dir = "/path/to/results/"
-work_dir = "/path/to/work/"
-ngen = 4
-nindividuals = 300
 
-target_depth_dict={'Node St1':{'target':1, 'weight':0.1}, 'Node St2':{'target':1.5, 'weight':0.1}}
+## configuration file
+The configuration file is a json file that specifies all of the parameters you will be using in your swmm_mpc run. There are certain parameters that are required to be specified and others that have a default value and are not required.
 
+### Required parameters in configuration file
+1. `inp_file_path`: [string] path to .inp file relative to config file
+2. `ctl_horizon`: [number] ctl horizon in hours
+3. `ctl_time_step`: [number] control time step in seconds
+4. `ctl_str_ids`: [list of strings] ids of control structures for which controls policies will be found. Each should start with one of the key words ORIFICE, PUMP, or WEIR e.g., [ORIFICE R1, ORIFICE R2]
+5. `work_dir`: [string] directory relative to config file where the temporary files will be created
+6. `results_dir`: [string] directory relative to config file where the results will be written
+7. `opt_method`: [string] optimization method. Currently supported methods are 'genetic_algorithm', and 'bayesian_opt'
+8. `optimization_params`: [dict] dictionary with key values that will be passed to the optimization function for GA this includes
+ * `ngen`: [int] number of generations for GA
+ * `nindividuals`: [int] number of individuals for initial generation in GA
+9. `run_suffix`: [string] will be added to the results filename
 
+### Optional parameters in configuration file
+1. `flood_weight`: [number] overall weight for the sum of all flooding relative to the overall weight for the sum of the absolute deviations from target depths (`dev_weight`). Default: 1
+2. `dev_weight`: [number] overall weight for the sum of the absolute deviations from target depths. This weight is relative to the `flood_weight` Default: 0
+3. `target_depth_dict`: [dict] dictionary where the keys are the nodeids and the values are a dictionary. The inner dictionary has two keys, 'target', and 'weight'. These values specify the target depth for the nodeid and the weight given to that in the cost function. e.g., {'St1': {'target': 1, 'weight': 0.1}} Default: None (flooding from all nodes is weighted equally)
+4. `node_flood_weight_dict`: [dict] dictionary where the keys are the node ids and the values are the relative weights for weighting the amount of flooding for a given node. e.g., {'st1': 10, 'J3': 1}. Default: None
 
-def main():
- run_swmm_mpc(inp_file,
- control_horizon,
- control_time_step,
- control_str_ids,
- work_dir,
- results_dir,
- target_depth_dict=target_depth_dict,
- ngen=ngen,
- nindividuals=nindividuals
- )
+## Example 
+### Example configuration file
+
+```
+{
+ "inp_file_path": "my_swmm_model.inp",
+ "ctl_horizon": 1,
+ "ctl_time_step": 900,
+ "ctl_str_ids": ["ORIFICE R1", "ORIFICE R2"],
+ "work_dir": "work/",
+ "results_dir": "results/",
+ "dev_weight":0.5,
+ "flood_weight":1000,
+ "run_suffix": "my_mpc_run",
+ "opt_method": "genetic_algorithm",
+ "optimization_params":
+ {
+ "num_cores":20,
+ "ngen":8,
+ "nindividuals":120
+ },
+ "target_depth_dict":
+ {
+ "Node St1":
+ {
+ "target":1.69,
+ "weight":1
+ },
+ "Node St2":
+ {
+ "target":1.69,
+ "weight":1
+ }
+ }
+} 
+```
+### Example python file
+```python
 
-if __name__ == "__main__":
- main()
+from swmm_mpc.swmm_mpc import run_swmm_mpc
+run_swmm_mpc('my_config_file.json')
 ```
 
 Then to run it, you simply call the script with python:
 ```
 python my_swmm_mpc.py
 ```
 # Dockerized code
-A Docker image with swmm_mpc and all of its dependencies can be found at [https://hub.docker.com/r/jsadler2/swmm_mpc/](https://hub.docker.com/r/jsadler2/swmm_mpc/). You would run it like so:
+A Docker image with swmm_mpc and all of its dependencies can be found at [https://hub.docker.com/r/jsadler2/swmm_mpc/](https://hub.docker.com/r/jsadler2/swmm_mpc/). You would run it like so (**this assumes your results\_dir, your workdir, your .inp file, and your config file (\*.json) are all in the same directory**): 
 
 ```
-docker run -v /path/to/inputfile.inp:/path/to/inputfile.inp -v /path/to/workdir/:/path/to/workdir/ -v /path/to/results/dir/:/path/to/results/dir/ -v /path/to/run_script.py:/run_script.py jsadler2/swmm_mpc:latest python /run_script.py
+docker run -v /path/to/run_dir/:/run_dir/ jsadler2/swmm_mpc:latest python /run_script.py
 ```
 # Example model 
 An example use case model is found on HydroShare: [https://www.hydroshare.org/resource/73b38d6417ac4352b9dae38a78a47d81/](https://www.hydroshare.org/resource/73b38d6417ac4352b9dae38a78a47d81/).

setup.py

@@ -13,6 +13,7 @@
  'matplotlib',
  'numpy',
  'deap',
+ 'GPyOpt',
  ],
  dependency_links=[
  'git+https://github.com/uva-hydroinformatics/pyswmm.git@feature_save_hotstart#egg=pyswmm-0',

swmm_mpc/evaluate.py

@@ -1,11 +1,13 @@
 import string
+import numpy as np
 import sys
 import random
 import os
 from shutil import copyfile
 import subprocess
 from rpt_ele import rpt_ele
 import update_process_model_input_file as up
+import swmm_mpc as sm
 
 
 def get_flood_cost_from_dict(rpt, node_flood_weight_dict):
@@ -49,6 +51,22 @@ def get_deviation_cost(rpt, target_depth_dict):
  return sum(node_deviation_costs)
 
 
+def get_cost(rpt_file, node_flood_weight_dict, flood_weight, target_depth_dict,
+ dev_weight):
+ # read the output file
+ rpt = rpt_ele('{}'.format(rpt_file))
+
+ # get flooding costs
+ node_fld_cost = get_flood_cost(rpt, node_flood_weight_dict)
+
+ # get deviation costs
+ deviation_cost = get_deviation_cost(rpt, target_depth_dict)
+
+ # convert the contents of the output file into a cost
+ cost = flood_weight*node_fld_cost + dev_weight*deviation_cost
+ return cost
+
+
 def bits_to_decimal(bits):
  bits_as_string = "".join(str(i) for i in bits)
  return float(int(bits_as_string, 2))
@@ -68,7 +86,7 @@ def bits_to_perc(bits):
 def bit_to_on_off(bit):
  """
  convert single bit to "ON" or "OFF"
- bit: [int] or [list] 
+ bit: [int] or [list]
  """
  if type(bit) == list:
  if len(bit) > 1:
@@ -107,7 +125,7 @@ def split_gene_by_ctl_ts(gene, control_str_ids, n_steps):
  # get the segment for the control
  gene_seg = gene[:n_bits]
  # split to get the different time steps
- gene_seg_per_ts = split_list(gene_seg, n_steps) 
+ gene_seg_per_ts = split_list(gene_seg, n_steps)
  # add the gene segment to the overall list
  split_gene.append(gene_seg_per_ts)
  # move the beginning of the gene to the end of the current ctl segment
@@ -153,67 +171,94 @@ def gene_to_policy_dict(gene, control_str_ids, n_control_steps):
  return fmted_policies
 
 
-def evaluate(individual, hs_file_path, process_file_path, sim_dt,
- control_time_step, n_control_steps, control_str_ids,
- node_flood_weight_dict, target_depth_dict, flood_weight,
- dev_weight):
- FNULL = open(os.devnull, 'w')
+def list_to_policy(policy, control_str_ids, n_control_steps):
+ """
+ ASSUMPTION: round decimal number to BOOLEAN
+ """
+ split_policies = split_list(policy, len(control_str_ids))
+ fmted_policies = dict()
+ for i, control_id in enumerate(control_str_ids):
+ control_type = control_id.split()[0]
+ if control_type == 'ORIFICE' or control_type == 'WEIR':
+ fmted_policies[control_id] = split_policies[i]
+ elif control_type == 'PUMP':
+ on_off = [bit_to_on_off(round(p)) for p in split_policies[i]]
+ fmted_policies[control_id] = on_off
+ return fmted_policies
+
+
+def format_policies(policy, control_str_ids, n_control_steps, opt_method):
+ if opt_method == 'genetic_algorithm':
+ return gene_to_policy_dict(policy, control_str_ids, n_control_steps)
+ elif opt_method == 'bayesian_opt':
+ return list_to_policy(policy, control_str_ids, n_control_steps)
+
+
+def prep_tmp_files(proc_inp, work_dir):
  # make process model tmp file
  rand_string = ''.join(random.choice(
  string.ascii_lowercase + string.digits) for _ in range(9))
 
  # make a copy of the process model input file
- process_file_dir, process_file_name = os.path.split(process_file_path)
- tmp_process_base = process_file_name.replace('.inp',
- '_tmp_{}_{}'.format(
- sim_dt,
- rand_string
- ))
- tmp_process_inp = os.path.join(process_file_dir,
- tmp_process_base + '.inp')
- tmp_process_rpt = os.path.join(process_file_dir,
- tmp_process_base + '.rpt')
- copyfile(process_file_path, tmp_process_inp)
+ tmp_proc_base = proc_inp.replace('.inp',
+ '_tmp_{}'.format(rand_string))
+ tmp_proc_inp = tmp_proc_base + '.inp'
+ tmp_proc_rpt = tmp_proc_base + '.rpt'
+ copyfile(proc_inp, tmp_proc_inp)
 
  # make copy of hs file
- hs_file_name = os.path.split(hs_file_path)[1]
+ hs_file_path = up.read_hs_filename(proc_inp)
+ hs_file_name = os.path.split(hs_file_path)[-1]
  tmp_hs_file_name = hs_file_name.replace('.hsf',
  '_{}.hsf'.format(rand_string))
+ tmp_hs_file_path = os.path.join(sm.run.work_dir, tmp_hs_file_name)
+ copyfile(hs_file_path, tmp_hs_file_path)
+ return tmp_proc_inp, tmp_proc_rpt, tmp_hs_file_path
 
- tmp_hs_file = os.path.join(process_file_dir, tmp_hs_file_name)
- copyfile(hs_file_path, tmp_hs_file)
 
- # convert individual to percentages
- fmted_policies = gene_to_policy_dict(individual, control_str_ids,
- n_control_steps)
+def evaluate(*individual):
+ """
+ evaluate the performance of an individual given the inp file of the process
+ model, the individual, the control params (ctl_str_ids, horizon, step),
+ and the cost function params (dev_weight/dict, flood weight/dict)
+ """
+ FNULL = open(os.devnull, 'w')
+ # prep files
+ tmp_inp, tmp_rpt, tmp_hs = prep_tmp_files(sm.run.inp_process_file_path,
+ sm.run.work_dir)
+
+ # format policies
+ if sm.run.opt_method == 'genetic_algorithm':
+ individual = individual[0]
+ elif sm.run.opt_method == 'bayesian_opt':
+ individual = np.squeeze(individual)
+
+ fmted_policies = format_policies(individual, sm.run.ctl_str_ids,
+ sm.run.n_ctl_steps, sm.run.opt_method)
 
  # update controls
- up.update_controls_and_hotstart(tmp_process_inp,
- control_time_step,
+ up.update_controls_and_hotstart(tmp_inp,
+ sm.run.ctl_time_step,
  fmted_policies,
- tmp_hs_file)
+ tmp_hs)
 
  # run the swmm model
  if os.name == 'nt':
  swmm_exe_cmd = 'swmm5.exe'
  elif sys.platform.startswith('linux'):
  swmm_exe_cmd = 'swmm5'
- cmd = '{} {} {}'.format(swmm_exe_cmd, tmp_process_inp,
- tmp_process_rpt)
+ cmd = '{} {} {}'.format(swmm_exe_cmd, tmp_inp,
+ tmp_rpt)
  subprocess.call(cmd, shell=True, stdout=FNULL, stderr=subprocess.STDOUT)
 
- # read the output file
- rpt = rpt_ele('{}'.format(tmp_process_rpt))
-
- # get flooding costs
- node_flood_cost = get_flood_cost(rpt, node_flood_weight_dict)
-
- # get deviation costs
- deviation_cost = get_deviation_cost(rpt, target_depth_dict)
-
- # convert the contents of the output file into a cost
- cost = flood_weight*node_flood_cost + dev_weight*deviation_cost
- os.remove(tmp_process_inp)
- os.remove(tmp_process_rpt)
- os.remove(tmp_hs_file)
- return cost,
+ # get cost
+ cost = get_cost(tmp_rpt,
+ sm.run.node_flood_weight_dict,
+ sm.run.flood_weight,
+ sm.run.target_depth_dict,
+ sm.run.dev_weight)
+
+ os.remove(tmp_inp)
+ os.remove(tmp_rpt)
+ os.remove(tmp_hs)
+ return cost

swmm_mpc/run_baeopt.py

@@ -0,0 +1,56 @@
+import GPyOpt
+import swmm_mpc as sm
+import evaluate as ev
+import numpy as np
+from GPyOpt.methods import BayesianOptimization as BayOpt
+
+
+def get_bounds(ctl_str_ids, nsteps):
+ bounds = []
+ for ctl in ctl_str_ids:
+ var_num = 0
+ for j in range(nsteps):
+ ctl_type = ctl.split()[0]
+ ctl_bounds = {}
+
+ if ctl_type == 'WEIR' or ctl_type == 'ORIFICE':
+ var_type = 'continuous'
+ elif ctl_type == 'PUMP':
+ var_type = 'discrete'
+
+ ctl_bounds['name'] = 'var_{}'.format(var_num)
+ ctl_bounds['type'] = var_type
+ ctl_bounds['domain'] = (0, 1)
+ var_num += 1
+ bounds.append(ctl_bounds)
+ return bounds
+
+
+def run_baeopt(opt_params):
+ # set up opt params
+ bounds = get_bounds(sm.run.ctl_str_ids, sm.run.n_ctl_steps)
+ max_iter = opt_params.get('max_iter', 15)
+ max_time = opt_params.get('max_time', 120)
+ initial_guess = opt_params.get('initial_guess', [])
+ if len(initial_guess) == 0:
+ initial_guess = None
+ else:
+ initial_guess = np.array([initial_guess])
+
+ eps = opt_params.get('eps', 0.01)
+ model_type = opt_params.get('model_type', 'GP')
+ acquisition_type = opt_params.get('acquisition_type', 'EI')
+
+ # instantiate object
+ bae_opt = BayOpt(ev.evaluate,
+ domain=bounds,
+ model_type=model_type,
+ acquisition_type='EI',
+ X=initial_guess,
+ evaluator_type='local_penalization',
+ num_cores=opt_params['num_cores'],
+ batch_size=opt_params['num_cores'],
+ )
+ bae_opt.run_optimization(max_iter, max_time, eps)
+ return bae_opt.x_opt, bae_opt.fx_opt
+