Commit a24609ed authored by Jesse VanderWees's avatar Jesse VanderWees 🐘

added a plot for ensemble results

parent e32785e3
......@@ -234,12 +234,10 @@ now the output for each ensemble is printed.
Extracting Results
------------------
The results from the run are found in ``output/timeseries_export.csv``. Any
CSV-reading software can import it, but this is what the results look like when
plotted in Microsoft Excel:
The results from the run are found in ``output/forcast1/timeseries_export.csv``
and ``output/forcast2/timeseries_export.csv``. Any CSV-reading software can
import it, but this is how results can be plotted using the python library
matplotlib:
.. note::
TODO: Plot these results
.. image:: ../images/lookuptable_resultplot.png
.. plot:: examples/pyplots/ensemble_results.py
:include-source:
......@@ -2,43 +2,51 @@ import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime
from pylab import get_cmap
data_path = '../../../examples/mixed_integer/output/timeseries_export.csv'
output_dir = '../../../examples/ensemble/output/'
forecast_names = ['forecast1', 'forecast2']
delimiter = ','
# Import Data
ncols = len(np.genfromtxt(data_path, max_rows=1, delimiter=delimiter))
datefunc = lambda x: datetime.strptime(x, '%Y-%m-%d %H:%M:%S')
results = np.genfromtxt(data_path, converters={0: datefunc}, delimiter=delimiter,
dtype='object' + ',float' * (ncols - 1), names=True)
def get_results(forecast_name):
data_path = output_dir + forecast_name + '/timeseries_export.csv'
ncols = len(np.genfromtxt(data_path, max_rows=1, delimiter=delimiter))
datefunc = lambda x: datetime.strptime(x, '%Y-%m-%d %H:%M:%S')
return np.genfromtxt(data_path, converters={0: datefunc}, delimiter=delimiter,
dtype='object' + ',float' * (ncols - 1), names=True)
# Generate Plot
f, axarr = plt.subplots(2, sharex=True)
axarr[0].set_title('Water Level and Discharge')
n_subplots = 2
f, axarr = plt.subplots(n_subplots, sharex=True, figsize=(8, 4 * n_subplots))
axarr[0].set_title('Water Volume and Discharge')
cmaps = ['Blues', 'Greens']
shades = [0.5, 0.8]
f.autofmt_xdate()
# Upper subplot
axarr[0].set_ylabel('Water Level [m]')
axarr[0].plot(results['time'], results['storage_level'], label='Storage',
linewidth=2, color='b')
axarr[0].plot(results['time'], results['sea_level'], label='Sea',
linewidth=2, color='m')
axarr[0].plot(results['time'], 0.5 * np.ones_like(results['time']), label='Storage Max',
linewidth=2, color='r', linestyle='--')
# Lower Subplot
axarr[1].set_ylabel('Flow Rate [m3/s]')
axarr[1].plot(results['time'], results['Q_orifice'], label='Orifice',
linewidth=2, color='g')
axarr[1].plot(results['time'], results['Q_pump'], label='Pump',
linewidth=2, color='r')
axarr[1].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
f.autofmt_xdate()
for idx, forecast in enumerate(forecast_names):
results = get_results(forecast)
axarr[0].set_ylabel('Water Volume in Storage [m3]')
if idx == 0:
axarr[0].plot(results['time'], results['V_max'], label='Max',
linewidth=2, color='r', linestyle='--')
axarr[0].plot(results['time'], results['V_min'], label='Min',
linewidth=2, color='g', linestyle='--')
axarr[0].plot(results['time'], results['V_storage'], label=forecast,
linewidth=2, color=get_cmap(cmaps[idx])(shades[1]))
# Lower Subplot
axarr[1].set_ylabel('Flow Rate [m3/s]')
axarr[1].plot(results['time'], results['Q_in'], label='{}:Inflow'.format(forecast),
linewidth=2, color=get_cmap(cmaps[idx])(shades[0]))
axarr[1].plot(results['time'], results['Q_release'], label='{}:Release'.format(forecast),
linewidth=2, color=get_cmap(cmaps[idx])(shades[1]))
# Shrink each axis by 20% and put a legend to the right of the axis
for i in range(len(axarr)):
box = axarr[i].get_position()
axarr[i].set_position([box.x0, box.y0, box.width * 0.8, box.height])
axarr[i].set_position([box.x0, box.y0, box.width * 0.7, box.height])
axarr[i].legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False)
# Output Plot
......
UTC,Q_in
2013-01-01 01:00:00,1.5
2013-01-02 01:00:00,1.5
2013-01-03 01:00:00,1.5
2013-01-04 01:00:00,1.5
2013-01-05 01:00:00,1.5
2013-01-06 01:00:00,1.5
2013-01-07 01:00:00,1.5
2013-01-08 01:00:00,1.5
2013-01-09 01:00:00,1.5
2013-01-10 01:00:00,1.5
2013-01-11 01:00:00,1.5
2013-01-12 01:00:00,1.5
2013-01-01 01:00:00,1.50
2013-01-02 01:00:00,1.95
2013-01-03 01:00:00,2.40
2013-01-04 01:00:00,2.85
2013-01-05 01:00:00,3.30
2013-01-06 01:00:00,3.75
2013-01-07 01:00:00,3.75
2013-01-08 01:00:00,3.30
2013-01-09 01:00:00,2.85
2013-01-10 01:00:00,2.40
2013-01-11 01:00:00,1.95
2013-01-12 01:00:00,1.50
[storage_V]
monotonicity = 1
[storage_H]
monotonicity = 1
storage_H,storage_V
-0.75,0
-0.5,0
-0.25,0
0,0
0.25,100000
0.5,500000
0.75,1000000
time,H_max,H_min,Q_in,Q_release,V_max,V_min,V_storage
2013-01-01 01:00:00,0.460000,0.440000,1.500000,0.000000,420860.554886,382552.000845,400000.000000
2013-01-02 01:00:00,0.460000,0.440000,1.950000,1.938478,420860.554886,382552.000845,400995.535049
2013-01-03 01:00:00,0.460000,0.440000,2.550000,2.550226,420860.554886,382552.000845,400975.998669
2013-01-04 01:00:00,0.460000,0.440000,3.300000,3.310621,420860.554886,382552.000845,400058.386209
2013-01-05 01:00:00,0.460000,0.440000,4.350000,4.531521,420860.554886,382552.000845,384374.935319
2013-01-06 01:00:00,0.460000,0.440000,5.550000,5.956544,420860.554886,382552.000845,349249.557569
2013-01-07 01:00:00,0.460000,0.440000,7.200000,6.000000,420860.554886,382552.000845,452929.552598
2013-01-08 01:00:00,0.460000,0.440000,5.550000,5.946848,420860.554886,382552.000845,418641.846678
2013-01-09 01:00:00,0.460000,0.440000,4.350000,4.527051,420860.554886,382552.000845,403344.610183
2013-01-10 01:00:00,0.460000,0.440000,3.300000,3.310475,420860.554886,382552.000845,402439.609290
2013-01-11 01:00:00,0.460000,0.440000,2.550000,2.546861,420860.554886,382552.000845,402710.828249
2013-01-12 01:00:00,0.460000,0.440000,1.950000,1.739819,420860.554886,382552.000845,420870.468899
time,H_max,H_min,Q_in,Q_release,V_max,V_min,V_storage
2013-01-01 01:00:00,0.460000,0.440000,1.500000,0.000000,420860.554886,382552.000845,400000.000000
2013-01-02 01:00:00,0.460000,0.440000,1.950000,1.936694,420860.554886,382552.000845,401149.627867
2013-01-03 01:00:00,0.460000,0.440000,2.400000,2.398757,420860.554886,382552.000845,401256.983825
2013-01-04 01:00:00,0.460000,0.440000,2.850000,2.849763,420860.554886,382552.000845,401277.484770
2013-01-05 01:00:00,0.460000,0.440000,3.300000,3.300393,420860.554886,382552.000845,401243.487870
2013-01-06 01:00:00,0.460000,0.440000,3.750000,3.744645,420860.554886,382552.000845,401706.182652
2013-01-07 01:00:00,0.460000,0.440000,3.750000,3.744645,420860.554886,382552.000845,402168.876481
2013-01-08 01:00:00,0.460000,0.440000,3.300000,3.300393,420860.554886,382552.000845,402134.878115
2013-01-09 01:00:00,0.460000,0.440000,2.850000,2.849764,420860.554886,382552.000845,402155.291593
2013-01-10 01:00:00,0.460000,0.440000,2.400000,2.398841,420860.554886,382552.000845,402255.459669
2013-01-11 01:00:00,0.460000,0.440000,1.950000,1.942015,420860.554886,382552.000845,402945.387750
2013-01-12 01:00:00,0.460000,0.440000,1.500000,1.292533,420860.554886,382552.000845,420870.499481
......@@ -62,25 +62,27 @@ class Example(GoalProgrammingMixin, ControlTreeMixin, CSVLookupTableMixin,
# Cache lookup tables for convenience and code legibility
_lookup_tables = self.lookup_tables(ensemble_member=0)
self.lookup_storage_V = _lookup_tables['storage_V']
self.lookup_storage_H = _lookup_tables['storage_H']
# Non-varrying goals can be implemented as a timeseries
for ensemble_member in range(self.ensemble_size):
self.set_timeseries('H_min',
np.ones_like(self.times()) * 0.44,
ensemble_member=ensemble_member)
self.set_timeseries('H_max',
np.ones_like(self.times()) * 0.46,
ensemble_member=ensemble_member)
# Convert our water level goals into volume goals
for ensemble_member in range(self.ensemble_size):
for e_m in range(self.ensemble_size):
self.set_timeseries('H_min', np.ones_like(self.times()) * 0.44,
ensemble_member=e_m)
self.set_timeseries('H_max', np.ones_like(self.times()) * 0.46,
ensemble_member=e_m)
# Q_in is a varying input and is defined in each timeseries_import.csv
# However, if we set it again here, it will be added to the output files
self.set_timeseries('Q_in',
self.get_timeseries('Q_in', ensemble_member=e_m),
ensemble_member=e_m)
# Convert our water level goals into volume goals
self.set_timeseries('V_max',
self.lookup_storage_V(self.get_timeseries('H_max')),
ensemble_member=ensemble_member)
ensemble_member=e_m)
self.set_timeseries('V_min',
self.lookup_storage_V(self.get_timeseries('H_min')),
ensemble_member=ensemble_member)
ensemble_member=e_m)
def path_goals(self):
g = []
......@@ -93,41 +95,31 @@ class Example(GoalProgrammingMixin, ControlTreeMixin, CSVLookupTableMixin,
# We want to print some information about our goal programming problem.
# We store the useful numbers temporarily, and print information at the
# end of our run.
for ensemble_member in range(self.ensemble_size):
results = self.extract_results(ensemble_member)
self.set_timeseries('H_storage',
self.lookup_storage_H(Timeseries(self.times(),
results['storage.V'])),
ensemble_member=ensemble_member)
self.set_timeseries('V_storage',
results['storage.V'],
ensemble_member=ensemble_member)
H_storage = self.get_timeseries('H_storage',
ensemble_member=ensemble_member).values
_max = self.get_timeseries('H_max',
ensemble_member=ensemble_member).values
_min = self.get_timeseries('H_min',
ensemble_member=ensemble_member).values
for e_m in range(self.ensemble_size):
results = self.extract_results(e_m)
self.set_timeseries('V_storage', results['storage.V'], ensemble_member=e_m)
_max = self.get_timeseries('V_max', ensemble_member=e_m).values
_min = self.get_timeseries('V_min', ensemble_member=e_m).values
V_storage = self.get_timeseries('V_storage', ensemble_member=e_m).values
# A little bit of tolerance when checking for acceptance. This
# tolerance must be set greater than the tolerance of the solver.
tol = 1e-3
tol = 10
_max += tol
_min -= tol
n_level_satisfied = sum(
np.logical_and(_min <= H_storage, H_storage <= _max))
np.logical_and(_min <= V_storage, V_storage <= _max))
q_release_integral = sum(results['Q_release'])
self.intermediate_results[ensemble_member].append((priority,
n_level_satisfied,
q_release_integral))
self.intermediate_results[e_m].append((priority, n_level_satisfied,
q_release_integral))
def post(self):
super(Example, self).post()
for ensemble_member in range(self.ensemble_size):
print('\n\nResults for Ensemble Member {}:'.format(ensemble_member))
for e_m in range(self.ensemble_size):
print('\n\nResults for Ensemble Member {}:'.format(e_m))
for priority, n_level_satisfied, q_release_integral in \
self.intermediate_results[ensemble_member]:
self.intermediate_results[e_m]:
print("\nAfter finishing goals of priority {}:".format(priority))
print("Level goal satisfied at {} of {} time steps".format(
n_level_satisfied, len(self.times())))
......
......@@ -89,14 +89,14 @@ class Example(GoalProgrammingMixin, CSVLookupTableMixin, CSVMixin,
_max = self.get_timeseries('V_max').values
_min = self.get_timeseries('V_min').values
values = self.get_timeseries('storage_V').values
storage_V = self.get_timeseries('storage_V').values
# A little bit of tolerance when checking for acceptance.
tol = 10
_max += tol
_min -= tol
n_level_satisfied = sum(
np.logical_and(_min <= values, values <= _max))
np.logical_and(_min <= storage_V, storage_V <= _max))
q_release_integral = sum(results['Q_release'])
self.intermediate_results.append(
(priority, n_level_satisfied, q_release_integral))
......@@ -104,11 +104,11 @@ class Example(GoalProgrammingMixin, CSVLookupTableMixin, CSVMixin,
def post(self):
# Call super() class to not overwrite default behaviour
super(Example, self).post()
for priority, n_level_satisfied, q_pump_integral in self.intermediate_results:
for priority, n_level_satisfied, q_release_integral in self.intermediate_results:
print("\nAfter finishing goals of priority {}:".format(priority))
print("Volume goal satisfied at {} of {} time steps".format(
n_level_satisfied, len(self.times())))
print("Integral of Q_release = {:.2f}".format(q_pump_integral))
print("Integral of Q_release = {:.2f}".format(q_release_integral))
# Any solver options can be set here
def solver_options(self):
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment