Skip to content
Snippets Groups Projects
Commit 15ad8b76 authored by Yu-Jin Schröer's avatar Yu-Jin Schröer
Browse files

Upload New File

parent 7226febb
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Optimization of control rules/Optimization_control_rules_load.ipynb 0 → 100644
+ 837
− 0
View file @ 15ad8b76
%% Cell type:markdown id:be61f6efdeea355d tags:
#### Import of libraries
%% Cell type:markdown id:e20bdc73 tags:
General libraries
%% Cell type:code id:e120f145 tags:
``` python
import typing
import datetime
import shutil
import subprocess
from multiprocessing.connection import Client
from time import sleep
from warnings import warn
from pathlib import Path
import functools
import threading
from multiprocessing.pool import ThreadPool
import traceback
import numpy as np
import pandas as pd
import nevergrad as ng
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d
from IPython.utils import io as ipyio
import os
```
%% Cell type:markdown id:13df68a6 tags:
Libraries for the interface with MIKE+
%% Cell type:code id:initial_id tags:
``` python
from mikeplus import DataTableAccess
from mikeplus.engines import Engine1D
from mikeio1d import Res1D
from mikeio1d.query import QueryDataStructure
from mikeio1d.query import QueryDataNode
```
%% Cell type:markdown id:d8ba7af66f30210 tags:
#### Setting variables
%% Cell type:code id:7f4a620fa3107d75 tags:
``` python
# the model_folder is the folder where the necessary input files are saved
# note this should be at the same folder depth as the original folder of this script, so that relative paths lead to the same files
model_folder = Path("statische Fracht Drosseloptimierung/")
# the name of the modelfile
db_path_rel = "240502_Modell Rodenkirchen.mupp"
# set the number of optimization steps
optimization_budget = 60
```
%% Cell type:markdown id:24d2e0b8 tags:
Names of the input Excel files
%% Cell type:code id:5275f8e7 tags:
``` python
event_xlsx = "config_events.xlsx"
parameter_xlsx = "config_parameter.xlsx"
structure_xlsx = "structures.xlsx" # The pollution of the discharges are specified here as well.
```
%% Cell type:markdown id:e8f8ec05dab0280 tags:
Further variables
%% Cell type:code id:2512aea6d68be310 tags:
``` python
# set the id of the sensors and the timezone
timezone = "UTC+01:00"
#Number of seats for the Engine, limits amount of simulation we can run at once
engine_seats = 4
# chose optimizer from nevergrade package
optimizers = ["BayesOptimBO"]
# delete result files(as specified by event config) before a sim, assures that no old file is reread in the case that a simulation does not run
DELETE_RESFILE_BEFORE_SIM = True
# output heaps of additional data, may be useful for debugging, Default = False
VERBOSE_OUTPUT = False
# attempt simulation at most so many times, used when a simulation fails
MAX_SIMULATION_TRY = 1
# Consider a simulation as failed and retry if so many values in the result compared to the reference are missing
MAX_MISSING_RATIO = 0.6
# delete the temp folder at end, Default = True
DELETE_TEMP = True
# time string for output folder
time_path_part = datetime.datetime.now().strftime("%Y-%m-%dT%H.%M.%S")
# set path of temp and output folder
temp_folder = Path("tmp_" + time_path_part)
# output folder
out_folder = Path("out_" + time_path_part)
# wheter output from mike simulation get shown on error
show_sim_out_on_error = True
```
%% Cell type:markdown id:c621349d1a2515d3 tags:
#### Prepartion of the database
%% Cell type:markdown id:4b52f384 tags:
file_to_df: Reads the file as a pandas Dataframe for excel, csv, res1d files
%% Cell type:code id:f8156368dd44c1bd tags:
``` python
def file_to_df(path: Path) -> pd.DataFrame:
if path.suffix == ".res1d":
res1d = Res1D(str(path)) # Create a Res1D instance with the file
df = res1d.read_all() # Read all data from the .res1d file
return df
elif path.suffix == ".xlsx":
return pd.read_excel(path)
elif path.suffix == ".csv":
return pd.read_csv(path)
else:
raise NotImplementedError(f"No method for {path.suffix} implemented")
```
%% Cell type:markdown id:a19c5fb4f02bff3e tags:
Creating output folder for results and temp folder where the model_folder is dublicated, to enshure that the origianl model is not damaged
%% Cell type:code id:460b9199eb10f580 tags:
``` python
out_folder.mkdir(exist_ok=True,parents=True)
assert not temp_folder.exists(),"Make sure not to mess with existing folders"
shutil.copytree(model_folder,temp_folder)
db_path = temp_folder / db_path_rel
# print the path of the temp and out folder
temp_folder.absolute()
out_folder.absolute()
```
%% Output
WindowsPath('C:/Users/Schroeer/Python Skripte/out_2025-02-21T10.40.32')
%% Cell type:markdown id:503240fbd162ad3f tags:
Read the excel config_events
%% Cell type:code id:185e7a9554f35dc0 tags:
``` python
events = file_to_df(temp_folder / event_xlsx).set_index('Event', drop=False)
# check if start time of events is before the end time.
assert (events["Start"] <= events["End"]).all(
axis=None
), "Event end needs to be after start."
# add timezone. Use same as sensor_date
events["Start"]=events["Start"].dt.tz_localize(timezone)
events["End"]=events["End"].dt.tz_localize(timezone)
# check if there are duplicated events
assert events.index.is_unique, "Need unique event"
assert events.drop_duplicates(subset="Event")["SimulationName"].is_unique, "Need exactly one simulation name per event"
# Conversion of the SimulationName column to the String data type, if it is not already there
events["SimulationName"] = events["SimulationName"].astype(str)
events
```
%% Cell type:markdown id:4edfb049ae1aa9bc tags:
Read the excel config_parameter
%% Cell type:code id:2825065f30d0835 tags:
``` python
param_specs = file_to_df(temp_folder / parameter_xlsx)
# give params_specs a sensible string index, this is later used to pass optimization parameters as a dictionary
param_specs.index = pd.Index(
param_specs[["Table", "Column", "Muid"]].agg(";".join, axis="columns"),
name="ParamKey",
)
assert (
param_specs.index.is_unique
), "Index needs to be unique. Otherwise indicates dupplicates in the parameters."
param_specs
```
%% Cell type:markdown id:7fc5cf01 tags:
Read the excel structures
%% Cell type:code id:9f364649 tags:
``` python
structures = file_to_df(temp_folder / structure_xlsx)
structures_list = structures['structures'].tolist()
structures_list
```
%% Cell type:code id:a6365d92 tags:
``` python
fracht_list = structures['fracht'].tolist()
fracht_list
```
%% Cell type:markdown id:e084bb2d61486c76 tags:
Open Mike+ database and start engine
%% Cell type:code id:506b4cc9696676fe tags:
``` python
dta = DataTableAccess(db_path)
# open the database from MIKE+
dta.open_database()
# check if the database is open.
dta.is_database_open()
# Create the Engine1D object that will be used to run MIKE 1D simulations
engine = Engine1D(dta.datatables)
```
%% Cell type:markdown id:9defc99104ef7ae9 tags:
Threading is used to run multiple simulations with differnt engines
%% Cell type:code id:2c7455cfdf9a4d80 tags:
``` python
engines = {}
def init_engine():
engines[threading.current_thread()] = Engine1D(dta.datatables)
print(f'initialized engine for {threading.current_thread()}')
# Use threadpool from the multiprocessing package, this has a rather simple interface
sim_executor = ThreadPool(processes=engine_seats,initializer=init_engine)
```
%% Cell type:markdown id:873214d0391ede5b tags:
#### Functions for the Simulation
%% Cell type:markdown id:679383a7d5714e0c tags:
calculate_discharge_for_structure: Calculate the discharge for a given structure form a given res1d object. Returns the total discharge in l.
%% Cell type:code id:9b28a84a6e9b2a97 tags:
``` python
def calculate_discharge_for_structure(res1d, structure):
try:
values = res1d.structures[structure].Discharge.read()
#convert ot liters
total_discharge = np.sum(values) * 60 * 1000
return total_discharge
except KeyError:
print(f'Structure {structure} not found in the Res1D data.')
return None
except AttributeError:
print (f'Discharge data not availave for strucutre {structure}.')
return None
except Exception as e:
print(f'An error occured while calculating discharge for structure {structure}: {e}')
return None
```
%% Cell type:markdown id:9f76cc017567a78b tags:
check_flood: It is checked for each node whether the value of waterVolumeAbouveGround is greater than zero, i.e. whether flooding is present.
To do this, the additional item ‘Water volume above ground’ must be selected in the MIKE+ result specification for the hydrodinamic result file.
%% Cell type:code id:1d2b04e24d558655 tags:
``` python
def check_flood(res1d) -> bool:
Flood = False
node_names = res1d.nodes
flooded_nodes = []
for node in node_names:
query_node = QueryDataNode(quantity='WaterVolumeAboveGround', name = node)
try:
water_level_values = query_node.get_values(res1d)
if water_level_values is None or len(water_level_values) == 0:
raise ValueError(f'No WaterVolumeAboveGround found for node {node}')
if (water_level_values > 0).any():
Flood = True
flooded_nodes.append(node)
print(f'Flood detected for node {node}')
except Exception as e:
print(f'Error for returning value for Node {node}: {e}')
if Flood:
print('Flood detected in the system!')
print('Flooded nodes:', flooded_nodes)
else:
print('No flood detected in the system!')
return Flood
```
%% Cell type:markdown id:e4df0cdeab58794e tags:
single_sim_run: Runs a single simulation of an event. Returns a sting massage on successfull finish.
%% Cell type:code id:245c0c41a582d710 tags:
``` python
def single_sim_run(ev):
sim_name = events.loc[events['Event'] == ev, 'SimulationName'].iloc[0]
print(f'Running sumlation {sim_name} in {threading.current_thread()}.')
engines[threading.current_thread()].run(sim_name)
return f'Completed simulation fo {ev} in {threading.current_thread()}.'
```
%% Cell type:markdown id:8bc732ab8dce7aa9 tags:
delet_result_file: deletes existing result file before a simulation corresponding to event.
%% Cell type:code id:64ecc39a9f9b792e tags:
``` python
def delete_result_file(event):
file_path = str(temp_folder / events.loc[event, "ResultFile"])
try:
Path(file_path).unlink(missing_ok=True)
except FileNotFoundError:
pass
```
%% Cell type:markdown id:612b3e8081cacd16 tags:
siumlation run: Takes parameters as arguments, runs simulations as defined by events with them and returns the results as a series with Multiindex
%% Cell type:code id:66215a7ad8996ecb tags:
``` python
def simulation_run(**kwparams):
assert set(kwparams.keys()) == set(param_specs.index)
timestamp = datetime.datetime.now().strftime("%Y-%m-%dT%H.%M.%S")
output_filename = f"simulation_results_{timestamp}.xlsx"
result_output_path = Path(out_folder / output_filename)
print(f"Results will be saved to: {result_output_path}")
for k in kwparams.keys():
tab = param_specs.loc[k, 'Table']
col = param_specs.loc[k, 'Column']
muid = param_specs.loc[k, 'Muid']
val = kwparams[k]
dta.set_value(tab, muid, col, val)
print(f'Set value for {tab}.{muid}.{col} to value: {val}')
if DELETE_RESFILE_BEFORE_SIM:
for l in events.index:
delete_result_file(l)
finish_msgs = sim_executor.map(single_sim_run,events["Event"].unique())
print(finish_msgs)
results_table = pd.DataFrame(columns=['Event', 'Structure', 'Discharge', 'Load'])
total_result = 0
for l in events.index:
result_file_path = temp_folder / events.loc[l, "ResultFile"]
result_file_path_str = str(Path(result_file_path))
# Check if the result file exists
if not Path(result_file_path).exists():
print(f"Result file does not exist: {result_file_path}")
continue
try:
# Create Res1D instance for each file
res1d = Res1D(result_file_path_str)
flood = check_flood(res1d)
#If there is flooding on the roads then a penalty will be issued
if flood == True:
results_ev = 10e50
else:
discharge_values = []
for structure, fracht in zip(structures_list, fracht_list):
discharge = calculate_discharge_for_structure(res1d, structure)
if discharge is not None:
weighted_discharge = discharge * fracht
discharge_values.append(weighted_discharge)
results_table = pd.concat([results_table, pd.DataFrame([{
'Event': events.loc[l, "Event"],
'Structure': structure,
'Discharge': discharge,
'Load': weighted_discharge
}])], ignore_index=True)
results_ev = np.sum(discharge_values) if discharge_values else 0
print(f"Result for result file {result_file_path}: {results_ev}")
total_result += results_ev
except Exception as e:
print(f"Error processing result file {result_file_path}: {e}")
traceback.print_exc()
results_ev = 0
print(f"Total result: {total_result}")
print(results_table)
results_table.to_excel(result_output_path, index=False)
print(f"Results saved to {result_output_path}")
return total_result
```
%% Cell type:markdown id:c6c65209016d9c72 tags:
#### Test simulation run of all events with the default parameters
%% Cell type:code id:7ba902b33d9088f9 tags:
``` python
param_specs.dtypes
_test_res = simulation_run(**param_specs["Default"])
display("_test_res:")
display(_test_res)
```
%% Cell type:markdown id:597bc0c3 tags:
#### preparations for optimisation
%% Cell type:markdown id:20ac8de37b5b66af tags:
Define optimization parameters as needed by nevergrade.
%% Cell type:code id:b5506c04386e1a93 tags:
``` python
ng_params = {}
for i in param_specs.index:
vmin = param_specs.at[i, "Min"]
vmax = param_specs.at[i, "Max"]
vdef = param_specs.at[i, "Default"]
ngp = ng.p.Scalar(init=vdef, lower=vmin, upper=vmax)
ng_params[i] = ngp
instrumentation = ng.p.Instrumentation(**ng_params)
instrumentation
```
%% Output
Instrumentation(Tuple(),Dict(msm_Pump;DutyPoint;MWP_0112=Scalar{Cl(0.001,0.12,b)}[sigma=Scalar{exp=2.03}],msm_Pump;DutyPoint;MWP_0142=Scalar{Cl(0.001,0.04,b)}[sigma=Scalar{exp=2.03}],msm_Pump;DutyPoint;MWP_0314_I=Scalar{Cl(0.001,0.9500000000000001,b)}[sigma=Scalar{exp=2.03}],msm_Pump;DutyPoint;MWP_0932_I=Scalar{Cl(0.001,0.083,b)}[sigma=Scalar{exp=2.03}],msm_Pump;DutyPoint;MWP_0944_I=Scalar{Cl(0.001,0.098,b)}[sigma=Scalar{exp=2.03}],msm_Pump;DutyPoint;RWP_0302=Scalar{Cl(0.001,0.5,b)}[sigma=Scalar{exp=2.03}],msm_Pump;DutyPoint;RWP_0335=Scalar{Cl(0.001,0.065,b)}[sigma=Scalar{exp=2.03}],msm_RTCAction;PIDSetPoint;HUB0201_rain=Scalar{Cl(0.062,0.518,b)}[sigma=Scalar{exp=2.03}],msm_RTCAction;PIDSetPoint;HUB0301_rain=Scalar{Cl(0.185,0.405,b)}[sigma=Scalar{exp=2.03}],msm_RTCAction;PIDSetPoint;S_SKU0216=Scalar{Cl(0.005,0.112,b)}[sigma=Scalar{exp=2.03}])):((), {'msm_Pump;DutyPoint;MWP_0142': 0.003, 'msm_Pump;DutyPoint;MWP_0112': 0.009000000000000001, 'msm_RTCAction;PIDSetPoint;HUB0301_rain': 0.261, 'msm_Pump;DutyPoint;MWP_0944_I': 0.073, 'msm_Pump;DutyPoint;MWP_0314_I': 0.004, 'msm_RTCAction;PIDSetPoint;HUB0201_rain': 0.371, 'msm_Pump;DutyPoint;RWP_0335': 0.059000000000000004, 'msm_RTCAction;PIDSetPoint;S_SKU0216': 0.012, 'msm_Pump;DutyPoint;MWP_0932_I': 0.001, 'msm_Pump;DutyPoint;RWP_0302': 0.001})
%% Cell type:markdown id:4523f7b62a1787ff tags:
params_losses_callback: returns two lists and a callback function. When the callback function is executed it puts the parameters and loss into the lists.!!Do not register this function itself as a callback, but the third returned parameter!!
%% Cell type:code id:6ef898b0e0627565 tags:
``` python
def params_losses_callback():
params = []
losses = []
def callback(optim, par, loss):
params.append(par)
losses.append(loss)
# returns two lists and the function
return params, losses, callback
```
%% Cell type:markdown id:cb86d6a2f0ea8c61 tags:
optim_run: Start the optimizer input = name of the optimizer output = dataframe with details from optimization model
%% Cell type:code id:a0e2e505ea5df65a tags:
``` python
def optim_run(optim_name):
# set the optimizer based on the optim_name and instrumentation from params_spec
optimizer = ng.optimizers.registry[optim_name](
parametrization=instrumentation, budget=optimization_budget, num_workers=1
)
# show progressbar. updated when the optimizer gets the result of a step
optimizer.register_callback("tell", ng.callbacks.ProgressBar())
# put params and loss in list
param_list, loss_list, plcb = params_losses_callback()
optimizer.register_callback("tell", plcb)
# get recommendation of parameters after run is finished
reccomendation = optimizer.minimize(simulation_run, verbosity=2)
optim_run_df = pd.concat(
[
pd.DataFrame(data=[p.value[1] for p in param_list]),
pd.DataFrame({"Optimizer": optimizer.name, "Loss": loss_list}),
],
axis="columns",
)
# also add a line for the reccomendation
rec_params = reccomendation.value[1]
optim_run_df.loc["reccomendation", rec_params.keys()] = rec_params
optim_run_df.loc["reccomendation", "Loss"] = reccomendation.loss
optim_run_df.loc["reccomendation", "Optimizer"] = optimizer.name
# save the results of the optimization as a csv and excel file
optim_run_df.to_csv(out_folder / f"optim_run{optimizer.name}.csv")
optim_run_df.to_excel(out_folder / f"optim_run{optimizer.name}.xlsx")
return optim_run_df
```
%% Cell type:markdown id:28e37afd tags:
#### Run optimisation
%% Cell type:markdown id:58dd27759c113b7e tags:
choose optimizers from nevergrad
https://en.wikipedia.org/wiki/CMA-ES
create a list with optimizers
it is also possible to use multiple optimizers and compare the results
create dictionary to store the rundata for each optimizer
%% Cell type:code id:1d8d2932f77d3de2 tags:
``` python
optim_run_dfs = {}
for o in optimizers:
# start the optimization
optim_run_dfs[o] = optim_run(o)
# display results
display(optim_run_dfs[o])
optim_run_dfs[o].plot(subplots=True, figsize=(5, 2 * len(param_specs)))
# optim_run_df already has optimizer included, use .values instead of passing the dict, so that the optimizer is not added again.
optim_run_df = pd.concat(optim_run_dfs.values())
optim_run_df.to_csv(out_folder / "optim_run.csv")
# set the index if the optim_run_df
optim_run_df.index.name = "Step"
optim_run_df = optim_run_df.reset_index().set_index(["Optimizer", "Step"])
```
%% Cell type:markdown id:0bb04f20 tags:
Print the Results of the Optimization
%% Cell type:code id:5a16015edf840d52 tags:
``` python
for key, df in optim_run_df.groupby(level="Optimizer"):
display(key,df.loc[(key, "reccomendation"), param_specs.index])
```
%% Cell type:markdown id:1875e60c tags:
Close the database
%% Cell type:code id:f4ba189120eb70e3 tags:
``` python
dta.close_database()
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment