RTCTools documentation¶
This is the documentation for RTCTools 2.0, the Deltares toolbox for control and optimization of environmental systems.
Visit the RTCTools website for a general product description and information on available services.
This first chapter covers getting the software running on your computer. The subsequent two chapters describe the RTCTools Python API. The fourth and final chapter discusses several illustrative examples, including the use of goal programming for multiobjective optimization, as well as the use of forecast ensembles.
Contents¶
Getting Started¶
Installation¶
For most users, especially on Windows, the easiest way to install RTCTools and its dependencies is using the Windows Installer. For users on Linux, it is necessary to build it from source using git.
Windows Installer¶
A complete Windows installer is available from the Deltares Download Portal. This will install a selfcontaining RTCTools 2 configuration including all prerequisites, and will not conflict with other Python or JModelica installations on the system.
 Download the installer from the Deltares Download Portal.
 Start the installation by opening the downloaded executable.
 Choose the desired location for RTCTools 2.
 Run the Basic Example via the Start Menu.
If the installation was succesful, you should see that the solver succeeds:
From Source¶
Although not required, it is recommended to build and install RTCTools and JModelica (see dependencies) in a virtual environment.
Dependencies¶
RTCTools 2 has the following system dependencies:
For most Linux distributions only Python is available from the standard repositories. JModelica and its dependencies will need to be built and installed from source. We refer to their respective installation instructions, and list below the instructions pertaining to RTCTools itself.
Acquiring the source¶
The latest RTCTools source can be downloaded using git:
# get RTCTools source
git clone https://gitlab.com/deltares/rtctools.git
# Get RTCTools's Modelica library
git clone https://gitlab.com/deltares/rtctoolschannelflow.git
Ubuntu / Debian¶
Building RTCTools requires one additional Python package over JMmodelica:
# Change directory to where RTCTools was downloaded
cd rtctools
# Install additional dependencies of RTCTools
pip install mock
Now all that remains is to actually build and install RTCTools:
python setup.py install
To check whether the installation was succesful, the basic example can be used. It is importent first to set the correct environment variables for JModelica and RTCTools. Luckily, JModelica comes with a convenient script which does most of this for you. Only the environment variable pointing to the Deltares Modelica library remains for the user to set:
export DELTARES_LIBRARY_PATH=\`readlink f ../rtctoolschannelflow\`
cd examples/basic/src
# Set the correct environment variables, and run the example
/path/to/JModelica/bin/jm_python.sh example.py
Windows¶
Building RTCTools on Windows is easiest by using the JModelica SDK. Be sure to:
 Build using JModelica’s CasADi 2.4 branch
 Update CasADi to the required version (see dependencies)
A further dependency is on Cython. Instructions for building extensions on Windows are available in the Cython docs.
Using the Visual C++ 2008 32bit Command Prompt, it is then possible to build and install RTCTools by running:
python setup.py install
To check whether the installation was succesful, the basic example can be used. It is importent first to set the correct environment variables for JModelica and RTCTools. Luckily, JModelica comes with a convenient script which does this for you. Only the environment variable pointing to the Deltares Modelica library remains for the user to set:
set DELTARES_LIBRARY_PATH=C:\path\to\rtctoolschannelflow
cd /D C:\path\to\rtctools\basic\src
# Set the correct environment variables, and run the example
C:\path\to\JModelica\Python.bat example.py
Getting OMEdit¶
RTCTools uses the Modelica language to describe the mathematics of the system we wish to optimize. There are several editors for Modelica models, but the OpenModelica Connection Editor, or OMEdit, is a free and opensource graphical connection editor that can be used to construct RTCTools models. To download it for windows, click here: https://www.openmodelica.org/download/downloadwindows
Once installed, you can start OMEdit by clicking:
Start > All Programs > OpenModelica > OpenModelica Connection Editor
With OMEdit installed, you can start using it by following along with the basic example, examples/basic.
Running RTCTools¶
RTCTools is run from a command line shell. If you installed using the Windows executable, the RTCTools Shell can be started by clicking:
Start > All Programs > RTCTools > Shell
Once you have started the shell, navigate to the src
directory of the case
you wish to optimize, e.g.:
cd \path\to\RTCTools2\examples\basic\src
Then, to run the case with RTCTools, run the src
python script, e.g.:
python example.py
You will see the progress of RTCTools in your shell. All your standard shell commands can be used in the RTCTools shell. For example, you can use:
python example.py > log.txt
to pipe RTCTools output to a log file.
Optimization¶
Contents:
Basics¶

class
rtctools.optimization.timeseries.
Timeseries
¶ Bases:
object
Time series object, bundling time stamps with values.

__init__
(self, times, values)¶ Create a new time series object.
Parameters:  times – Iterable of time stamps.
 values – Iterable of values.

times
¶ Array of time stamps.

values
¶ Array of values.


class
rtctools.optimization.optimization_problem.
OptimizationProblem
¶ Bases:
object
Base class for all optimization problems.

bounds
(self)¶ Returns variable bounds as a dictionary mapping variable names to a pair of bounds. A bound may be a constant, or a time series.
Returns: A dictionary of variable names and (upper, lower)
bound pairs. The bounds may be numbers orTimeseries
objects.Example:
def bounds(self): return {'x': (1.0, 2.0), 'y': (2.0, 3.0)}

constant_inputs
(self, ensemble_member)¶ Returns a dictionary of constant inputs.
Parameters: ensemble_member – The ensemble member index. Returns: A dictionary of constant input names and time series.

constraints
(self, ensemble_member)¶ Returns a list of constraints for the given ensemble member.
Call
OptimizationProblem.state_at()
to return a symbol representing a model variable at a given time.Parameters: ensemble_member – The ensemble member index. Returns: A list of triples (f, m, M)
, with anMX
object representing the constraint functionf
, lower boundm
, and upper boundM
. The bounds must be numbers.Example:
def constraints(self, ensemble_member): t = 1.0 constraint1 = (2 * self.state_at('x', t, ensemble_member), 2.0, 4.0) constraint2 = (self.state_at('x', t, ensemble_member) + self.state_at('y', t, ensemble_member), 2.0, 3.0) return [constraint1, constraint2]

control
(self, variable)¶ Returns an
MX
symbol for the given control input, not bound to any time.Parameters: variable – Variable name. Returns: MX
symbol for given control input.Raises: KeyError

control_at
(self, variable, t, ensemble_member=0, scaled=False)¶ Returns an
MX
symbol representing the given control input at the given time.Parameters:  variable – Variable name.
 t – Time.
 ensemble_member – The ensemble member index.
 scaled – True to return the scaled variable.
Returns: MX
symbol representing the control input at the given time.Raises: KeyError

delayed_feedback
(self)¶ Returns the delayed feedback mappings. These are given as a list of triples \((x, y, \tau)\), to indicate that \(y = x(t  \tau)\).
Returns: A list of triples. Example:
def delayed_feedback(self): fb1 = ['x', 'y', 0.1] fb2 = ['x', 'z', 0.2] return [fb1, fb2]

der
(self, variable)¶ Returns an
MX
symbol for the time derivative given state, not bound to any time.Parameters: variable – Variable name. Returns: MX
symbol for given state.Raises: KeyError

der_at
(self, variable, t, ensemble_member=0)¶ Returns an expression for the time derivative of the specified variable at time t.
Parameters:  variable – Variable name.
 t – Time.
 ensemble_member – The ensemble member index.
Returns: MX
object representing the derivative.Raises: KeyError

ensemble_member_probability
(self, ensemble_member)¶ The probability of an ensemble member occurring.
Parameters: ensemble_member – The ensemble member index. Returns: The probability of an ensemble member occurring. Raises: IndexError

ensemble_size
¶ The number of ensemble members.

get_timeseries
(self, variable, ensemble_member=0)¶ Looks up a timeseries from the internal data store.
Parameters:  variable – Variable name.
 ensemble_member – The ensemble member index.
Returns: The requested time series.
Return type: Timeseries
Raises: KeyError

history
(self, ensemble_member)¶ Returns the state history. Uses the initial_state() method by default.
Parameters: ensemble_member – The ensemble member index. Returns: A dictionary of variable names and historical time series (up to and including t0).

initial_time
¶ The initial time in seconds.

integral
(self, variable, t0=None, tf=None, ensemble_member=0)¶ Returns an expression for the integral over the interval [t0, tf].
Parameters:  variable – Variable name.
 t0 – Left bound of interval. If equal to None, the initial time is used.
 tf – Right bound of interval. If equal to None, the final time is used.
 ensemble_member – The ensemble member index.
Returns: MX
object representing the integral.Raises: KeyError

interpolate
(self, t, ndarray ts, fs, f_left=np.nan, f_right=np.nan, mode=INTERPOLATION_LINEAR)¶ Linear interpolation over time.
Parameters:  t (float or vector of floats) – Time at which to evaluate the interpolant.
 ts (numpy array) – Time stamps.
 fs – Function values at time stamps ts.
 f_left – Function value left of leftmost time stamp.
 f_right – Function value right of rightmost time stamp.
 mode – Interpolation mode.
Returns: The interpolated value.

lookup_tables
(self, ensemble_member)¶ Returns a dictionary of lookup tables.
Parameters: ensemble_member – The ensemble member index. Returns: A dictionary of variable names and lookup tables.

objective
(self, ensemble_member)¶ The objective function for the given ensemble member.
Call
OptimizationProblem.state_at()
to return a symbol representing a model variable at a given time.Parameters: ensemble_member – The ensemble member index. Returns: An MX
object representing the objective function.Example:
def objective(self, ensemble_member): # Return value of state 'x' at final time: times = self.times() return self.state_at('x', times[1], ensemble_member)

optimize
(self, preprocessing=True, postprocessing=True, log_solver_failure_as_error=True)¶ Perform one initializetranscribesolvefinalize cycle.
Parameters:  preprocessing – True to enable a call to
pre
preceding the opimization.  postprocessing – True to enable a call to
post
following the optimization.
Returns: True on success.
 preprocessing – True to enable a call to

parameters
(self, ensemble_member)¶ Returns a dictionary of parameters.
Parameters: ensemble_member – The ensemble member index. Returns: A dictionary of parameter names and values.

path_constraints
(self, ensemble_member)¶ Returns a list of path constraints. Path constraints apply to all times and ensemble members simultaneously.
Call
OptimizationProblem.state()
to return a time and ensemblememberindependent symbol representing a model variable.Parameters: ensemble_member – The ensemble member index. This index may only be used to supply memberdependent bounds. Returns: A list of triples (f, m, M)
, with anMX
object representing the path constraint functionf
, lower boundm
, and upper boundM
. The bounds may be numbers orTimeseries
objects.Example:
def path_constraints(self, ensemble_member): # 2 * x must lie between 2 and 4 for every time instance. path_constraint1 = (2 * self.state('x'), 2.0, 4.0) # x + y must lie between 2 and 3 for every time instance path_constraint2 = (self.state('x') + self.state('y'), 2.0, 3.0) return [path_constraint1, path_constraint2]

path_objective
(self, ensemble_member)¶ Returns a path objective the given ensemble member. Path objectives apply to all times and ensemble members simultaneously.
Call
OptimizationProblem.state()
to return a time and ensemblememberindependent symbol representing a model variable.Parameters: ensemble_member – The ensemble member index. This index is currently unused, and here for future use only. Returns: A MX
object representing the path objective.Example:
def path_objective(self, ensemble_member): # Minimize x(t) for all t return self.state('x')

post
(self)¶ Postprocessing logic is performed here.

pre
(self)¶ Preprocessing logic is performed here.

seed
(self, ensemble_member)¶ Seeding data. The optimization algorithm is seeded with the data returned by this method.
Parameters: ensemble_member – The ensemble member index. Returns: A dictionary of variable names and seed time series.

set_timeseries
(self, variable, timeseries, ensemble_member=0, output=True, check_consistency=True)¶ Sets a timeseries in the internal data store.
Parameters:  variable – Variable name.
 timeseries (iterable of floats, or
Timeseries
) – Time series data.  ensemble_member – The ensemble member index.
 output – Whether to include this time series in output data files.
 check_consistency – Whether to check consistency between the time stamps on the new timeseries object and any existing time stamps.

solver_options
(self)¶ Returns a dictionary of CasADi optimization problem solver options.
The default solver for continuous problems is Ipopt. The default solver for mixed integer problems is Bonmin.
Returns: A dictionary of CasADi NlpSolver
options. See the CasADi, Ipopt, and Bonmin documentation for details.

state
(self, variable)¶ Returns an
MX
symbol for the given state, not bound to any time.Parameters: variable – Variable name. Returns: MX
symbol for given state.Raises: KeyError

state_at
(self, variable, t, ensemble_member=0, scaled=False)¶ Returns an
MX
symbol representing the given variable at the given time.Parameters:  variable – Variable name.
 t – Time.
 ensemble_member – The ensemble member index.
 scaled – True to return the scaled variable.
Returns: MX
symbol representing the state at the given time.Raises: KeyError

states_in
(self, variable, t0=None, tf=None, ensemble_member=0)¶ Iterates over symbols for states in the interval [t0, tf].
Parameters:  variable – Variable name.
 t0 – Left bound of interval. If equal to None, the initial time is used.
 tf – Right bound of interval. If equal to None, the final time is used.
 ensemble_member – The ensemble member index.
Raises: KeyError

timeseries_at
(self, variable, t, ensemble_member=0)¶ Return the value of a time series at the given time.
Parameters:  variable – Variable name.
 t – Time.
 ensemble_member – The ensemble member index.
Returns: The interpolated value of the time series.
Raises: KeyError


rtctools.util.
run_optimization_problem
(optimization_problem_class, base_folder='..', log_level=logging.INFO, profile=False, profile_casadi=False)¶ Sets up and solves an optimization problem.
This function makes the following assumptions:
 That the
base_folder
contains subfoldersinput
,output
, andmodel
, containing input data, output data, and the model, respectively.  When using
CSVLookupTableMixin
, that the base folder contains a subfolderlookup_tables
.  When using
ModelicaMixin
, that the base folder contains a subfoldermodel
.  When using
ModelicaMixin
, that the toplevel Modelica model name equals the class name.
Parameters:  optimization_problem_class – Optimization problem class to solve.
 base_folder – Base folder.
 log_level – The log level to use.
 profile – Whether or not to enable profiling.
 profile_casadi – Whether or not to enable CasADi profiling.
 That the
Time discretization¶

class
rtctools.optimization.collocated_integrated_optimization_problem.
CollocatedIntegratedOptimizationProblem
¶ Bases:
rtctools.optimization.optimization_problem.OptimizationProblem
Discretizes your model using a mixed collocation/integration scheme.
Collocation means that the discretized model equations are included as constraints between state variables in the optimization problem.
Note
To ensure that your optimization problem only has globally optimal solutions, any model equations that are collocated must be linear. By default, all model equations are collocated, and linearity of the model equations is verified. Working with nonlinear models is possible, but discouraged.
Variables: check_collocation_linearity – If True
, check whether collocation constraints are linear. Default isTrue
.
integrated_states
¶ A list of states that are integrated rather than collocated.
Warning
This is an experimental feature.

integrator_options
(self)¶ Configures the implicit function used for time step integration.
Returns: A dictionary of CasADi ImplicitFunction
options. See the CasADi documentation for details.

interpolation_method
(self, variable=None)¶ Interpolation method for variable.
Parameters: variable – Variable name. Returns: Interpolation method for the given variable.

theta
¶ RTCTools discretizes differential equations of the form
\[\dot{x} = f(x, u)\]using the \(\theta\)method
\[x_{i+1} = x_i + \Delta t \left[\theta f(x_{i+1}, u_{i+1}) + (1  \theta) f(x_i, u_i)\right]\]The default is \(\theta = 1\), resulting in the implicit or backward Euler method. Note that in this case, the control input at the initial time step is not used.
Set \(\theta = 0\) to use the explicit or forward Euler method. Note that in this case, the control input at the final time step is not used.
Warning
This is an experimental feature for \(0 < \theta < 1\).

times
(self, variable=None)¶ List of time stamps for variable.
Parameters: variable – Variable name. Returns: A list of time stamps for the given variable.

Modelica models¶
To learn the basics of modelling with Modelica, please refer to the online book Modelica by Example.

class
rtctools.optimization.modelica_mixin.
ModelicaMixin
¶ Bases:
rtctools.optimization.optimization_problem.OptimizationProblem
Adds a Modelica model to your optimization problem.
During preprocessing, the Modelica files located inside the
model
subfolder are loaded.Variables: modelica_library_folder – Folder in which any referenced Modelica libraries are to be found. Default is mo
.
CSV I/O¶

class
rtctools.optimization.csv_mixin.
CSVMixin
¶ Bases:
rtctools.optimization.optimization_problem.OptimizationProblem
Adds reading and writing of CSV timeseries and parameters to your optimization problem.
During preprocessing, files named
timeseries_import.csv
,initial_state.csv
, andparameters.csv
are read from theinput
subfolder.During postprocessing, a file named
timeseries_export.csv
is written to theoutput
subfolder.In ensemble mode, a file named
ensemble.csv
is read from theinput
folder. This file contains two columns. The first column gives the name of the ensemble member, and the second column its probability. Furthermore, the other XML files appear one level deeper inside the filesystem hierarchy, inside subfolders with the names of the ensemble members.Variables:  csv_delimiter – Column delimiter used in CSV files. Default is
,
.  csv_equidistant – Whether or not the timeseries data is equidistant. Default is
True
.  csv_ensemble_mode – Whether or not to use ensembles. Default is
False
.  csv_validate_timeseries – Check consistency of timeseries. Default is
True
.

max_timeseries_id
(self, variable)¶ Returns the name of the upper bound timeseries for the specified variable.
Parameters: variable – Variable name.

min_timeseries_id
(self, variable)¶ Returns the name of the lower bound timeseries for the specified variable.
Parameters: variable – Variable name.
 csv_delimiter – Column delimiter used in CSV files. Default is
DelftFEWS I/O¶

class
rtctools.optimization.pi_mixin.
PIMixin
¶ Bases:
rtctools.optimization.optimization_problem.OptimizationProblem
Adds DelftFEWS Published Interface I/O to your optimization problem.
During preprocessing, files named
rtcDataConfig.xml
,timeseries_import.xml
,rtcParameterConfig.xml
, andrtcParameterConfig_Numerical.xml
are read from theinput
subfolder.rtcDataConfig.xml
maps tuples of FEWS identifiers, including location and parameter ID, to RTCTools time series identifiers.During postprocessing, a file named
timeseries_export.xml
is written to theoutput
subfolder.Variables:  pi_binary_timeseries – Whether to use PI binary timeseries format. Default is
False
.  pi_parameter_config_basenames – List of parameter config file basenames to read. Default is [
rtcParameterConfig
].  pi_parameter_config_numerical_basename – Numerical config file basename to read. Default is
rtcParameterConfig_Numerical
.  pi_check_for_duplicate_parameters – Check if duplicate parameters are read. Default is
False
.  pi_validate_timeseries – Check consistency of timeseries. Default is
True
.

max_timeseries_id
(self, variable)¶ Returns the name of the upper bound timeseries for the specified variable.
Parameters: variable – Variable name

min_timeseries_id
(self, variable)¶ Returns the name of the lower bound timeseries for the specified variable.
Parameters: variable – Variable name
 pi_binary_timeseries – Whether to use PI binary timeseries format. Default is
Lookup tables¶

class
rtctools.optimization.optimization_problem.
LookupTable
¶ Bases:
object
Lookup table.

__call__
(self, *args)¶ Evaluate the lookup table.
Parameters: args (Float, iterable of floats, or Timeseries
) – Input values.Returns: Lookup table evaluated at input values. Example use:
y = lookup_table(1.0) [y1, y2] = lookup_table([1.0, 2.0])


class
rtctools.optimization.csv_lookup_table_mixin.
CSVLookupTableMixin
¶ Bases:
rtctools.optimization.optimization_problem.OptimizationProblem
Adds lookup tables to your optimization problem.
During preprocessing, the CSV files located inside the
lookup_tables
subfolder are read. In every CSV file, the first column contains the output of the lookup table. Subsequent columns contain the input variables.Cubic BSplines are used to turn the data points into continuous lookup tables.
Optionally, a file
curvefit_options.ini
may be included inside thelookup_tables
folder. This file contains, grouped per lookup table, the following options: monotonicity:
 is an integer, magnitude is ignored
 if positive, causes spline to be monotonically increasing
 if negative, causes spline to be monotonically decreasing
 if 0, leaves spline monotonicity unconstrained
 curvature:
 is an integer, magnitude is ignored
 if positive, causes spline curvature to be positive (convex)
 if negative, causes spline curvature to be negative (concave)
 if 0, leaves spline curvature unconstrained
Note
Currently only onedimensional lookup tables are fully supported. Support for twodimensional lookup tables is experimental.
Variables:  csv_delimiter – Column delimiter used in CSV files. Default is
,
.  csv_lookup_table_debug – Whether to generate plots of the spline fits. Default is
false
.  csv_lookup_table_debug_points – Number of evaluation points for plots. Default is
100
.
Treatment of nonconvexities using homotopy¶
Using homotopy, a convex optimization problem can be continuously deformed into a nonconvex problem.

class
rtctools.optimization.homotopy_mixin.
HomotopyMixin
¶ Bases:
rtctools.optimization.optimization_problem.OptimizationProblem
Adds homotopy to your optimization problem. A homotopy is a continuous transformation between two optimization problems, parametrized by a single parameter \(\theta \in [0, 1]\).
Homotopy may be used to solve nonconvex optimization problems, by starting with a convex approximation at \(\theta = 0.0\) and ending with the nonconvex problem at \(\theta = 1.0\).
Note
It is advised to look for convex reformulations of your problem, before resorting to a use of the (potentially expensive) homotopy process.

homotopy_options
(self)¶ Returns a dictionary of options controlling the homotopy process.
Option Type Default value delta_theta_0
float
1.0
delta_theta_min
float
0.01
homotopy_parameter
string
theta
The homotopy process is controlled by the homotopy parameter in the model, specified by the option
homotopy_parameter
. The homotopy parameter is initialized to0.0
, and increases to a value of1.0
with a dynamically changing step size. This step size is initialized with the value of the optiondelta_theta_0
. If this step size is too large, i.e., if the problem with the increased homotopy parameter fails to converge, the step size is halved. The process of halving terminates when the step size falls below the minimum value specified by the optiondelta_theta_min
.Returns: A dictionary of homotopy options.

Initial state estimation¶

class
rtctools.optimization.initial_state_estimation_mixin.
InitialStateEstimationMixin
¶ Bases:
rtctools.optimization.goal_programming_mixin.GoalProgrammingMixin
Adds initial state estimation to your optimization problem using goal programming.
Before any other goals are evaluated, first, the deviation between initial state measurements and their respective model states is minimized in the least squares sense (1DVAR, priority 2). Secondly, the distance between pairs of states is minimized, again in the least squares sense, so that “smooth” initial guesses are provided for states without measurements (priority 1).
Note
There are types of problems where, in addition to minimizing differences between states and measurements, it is advisable to perform a steadystate initialization using additional initialtime model equations. For hydraulic models, for instance, it is often helpful to require that the timederivative of the flow variables vanishes at the initial time.

initial_state_measurements
(self)¶ List of pairs
(state, measurement_id)
or triples(state, measurement_id, max_deviation)
, relating states to measurement time series IDs.The default maximum deviation is
1.0
.

initial_state_smoothing_pairs
(self)¶ List of pairs
(state1, state2)
or triples(state1, state2, max_deviation)
, relating states the distance of which is to be minimized.The default maximum deviation is
1.0
.

Multiobjective optimization¶

class
rtctools.optimization.goal_programming_mixin.
Goal
¶ Bases:
object
Base class for lexicographic goal programming goals.
A goal is defined by overriding the
function()
method, and setting at least thefunction_range
class variable.Variables:  function_range – Range of goal function. Required.
 function_nominal – Nominal value of function. Used for scaling. Default is
1
.  target_min – Desired lower bound for goal function. Default is
numpy.nan
.  target_max – Desired upper bound for goal function. Default is
numpy.nan
.  priority – Integer priority of goal. Default is
1
.  weight – Optional weighting applied to the goal. Default is
1.0
.  order – Penalization order of goal violation. Default is
2
.  critical – If
True
, the algorithm will abort if this goal cannot be fully met. Default isFalse
.
The target bounds indicate the range within the function should stay, if possible. Goals are, in that sense, soft, as opposed to standard hard constraints.
Four types of goals can be created:
Minimization goal if no target bounds are set:
\[\min f\]Lower bound goal if
target_min
is set:\[m \leq f\]Upper bound goal if
target_max
is set:\[f \leq M\]Combined lower and upper bound goal if
target_min
andtarget_max
are both set:\[m \leq f \leq M\]
Lower priority goals take precedence over higher priority goals.
Goals with the same priority are weighted off against each other in a single objective function.
The goal violation value is taken to the order’th power in the objective function of the final optimization problem.
Example definition of the point goal \(x(t) \geq 1.1\) for \(t=1.0\) at priority 1:
class MyGoal(Goal): def function(self, optimization_problem, ensemble_member): # State 'x' at time t = 1.0 t = 1.0 return optimization_problem.state_at('x', t, ensemble_member) function_range = (1.0, 2.0) target_min = 1.1 priority = 1
Example definition of the path goal \(x(t) \geq 1.1\) for all \(t\) at priority 2:
class MyPathGoal(Goal): def function(self, optimization_problem, ensemble_member): # State 'x' at any point in time return optimization_problem.state('x') function_range = (1.0, 2.0) target_min = 1.1 priority = 2
Note that for path goals, the ensemble member index is not passed to the call to
OptimizationProblem.state()
. This call returns a timeindependent symbol that is also independent of the active ensemble member. Path goals are applied to all times and all ensemble members simultaneously.
function
(self, optimization_problem, ensemble_member)¶ This method returns a CasADi
MX
object describing the goal function.Returns: A CasADi MX
object.

get_function_key
(self, optimization_problem, ensemble_member)¶ Returns a key string uniquely identifying the goal function. This is used to eliminate linearly dependent constraints from the optimization problem.

has_target_bounds
¶ True
if the user goal has min/max bounds.

has_target_max
¶ True
if the user goal has max bounds.

has_target_min
¶ True
if the user goal has min bounds.

is_empty
¶

class
rtctools.optimization.goal_programming_mixin.
StateGoal
¶ Bases:
rtctools.optimization.goal_programming_mixin.Goal
Base class for lexicographic goal programming path goals that act on a single model state.
A state goal is defined by setting at least the
state
class variable.Variables:  state – State on which the goal acts. Required.
 target_min – Desired lower bound for goal function. Default is
numpy.nan
.  target_max – Desired upper bound for goal function. Default is
numpy.nan
.  priority – Integer priority of goal. Default is
1
.  weight – Optional weighting applied to the goal. Default is
1.0
.  order – Penalization order of goal violation. Default is
2
.  critical – If
True
, the algorithm will abort if this goal cannot be fully met. Default isFalse
.
Example definition of the goal \(x(t) \geq 1.1\) for all \(t\) at priority 2:
class MyStateGoal(StateGoal): state = 'x' target_min = 1.1 priority = 2
Contrary to ordinary
Goal
objects,PathGoal
objects need to be initialized with anOptimizationProblem
instance to allow extraction of state metadata, such as bounds and nominal values. Consequently, state goals must be instantiated as follows:my_state_goal = MyStateGoal(optimization_problem)
Note that
StateGoal
is a helper class. State goals can also be defined usingGoal
as direct base class, by implementing thefunction
method and providing thefunction_range
andfunction_nominal
class variables manually.
__init__
(self, optimization_problem)¶ Initialize the state goal object.
Parameters: optimization_problem – OptimizationProblem
instance.

class
rtctools.optimization.goal_programming_mixin.
GoalProgrammingMixin
¶ Bases:
rtctools.optimization.optimization_problem.OptimizationProblem
Adds lexicographic goal programming to your optimization problem.

goal_programming_options
(self)¶ Returns a dictionary of options controlling the goal programming process.
Option Type Default value constraint_relaxation
float
0.0
violation_tolerance
float
inf
mu_reinit
bool
True
fix_minimized_values
bool
True
check_monotonicity
bool
True
equality_threshold
float
1e8
Constraints generated by the goal programming algorithm are relaxed by applying the specified relaxation. Use of this option is normally not required.
A goal is considered to be violated if the violation, scaled between 0 and 1, is greater than the specified tolerance. Violated goals are fixed. Use of this option is normally not required.
The Ipopt barrier parameter
mu
is normally reinitialized a every iteration of the goal programming algorithm, unless mu_reinit is set toFalse
. Use of this option is normally not required.If
fix_minimized_values
is set toTrue
, goal functions will be set to equal their optimized values in optimization problems generated during subsequent priorities. Otherwise, only as upper bound will be set. Use of this option is normally not required.If
check_monotonicity
is set toTrue
, then it will be checked whether goals with the same function key form a monotonically decreasing sequence with regards to the target interval.The option
equality_threshold
controls when a twosided inequality constraint is folded into an equality constraint.The option
interior_distance
controls the distance from the scaled target bounds, starting from which the function value is considered to lie in the interior of the target space.Returns: A dictionary of goal programming options.

goals
(self)¶ User problem returns list of
Goal
objects.Returns: A list of goals.

path_goals
(self)¶ User problem returns list of path
Goal
objects.Returns: A list of path goals.

priority_completed
(self, priority)¶ Called after optimization for goals of certain priority is completed.
Parameters: priority – The priority level that was completed.

priority_started
(self, priority)¶ Called when optimization for goals of certain priority is started.
Parameters: priority – The priority level that was started.

Forecast uncertainty¶

class
rtctools.optimization.control_tree_mixin.
ControlTreeMixin
¶ Bases:
rtctools.optimization.optimization_problem.OptimizationProblem
Adds a stochastic control tree to your optimization problem.

control_tree_options
(self)¶ Returns a dictionary of options controlling the creation of a kary stochastic tree.
Option Type Default value forecast_variables
list
of stringsAll constant inputs branching_times
list
of floatsself.times()
k
int
2
A
k
ary tree is generated, branching at every interior branching time. Ensemble members are clustered to paths through the tree based on average distance over all forecast variables.Returns: A dictionary of control tree generation options.

Simulation¶
Note
For a simulation example, see Simulation examples
Contents:
Basics¶

class
rtctools.simulation.simulation_problem.
SimulationProblem
¶ Bases:
object
FMU simulation runner.
Implements the BMI Interface.

compiler_options
(self)¶ Subclasses can configure the JModelica.org compiler options here.
Returns: A dictionary of JModelica.org compiler options. See the JModelica.org documentation for details.

finalize
(self)¶ Finalize FMU.

get_current_time
(self)¶ Return current time of simulation.
Returns: The current simulation time.

get_end_time
(self)¶ Return end time of experiment.
Returns: The end time of the experiment.

get_input_variables
(self)¶

get_options
(self)¶ Return the available options of the FMU.
Returns: A dictionary of options supported by the FMU.

get_output_variables
(self)¶

get_parameter_variables
(self)¶

get_start_time
(self)¶ Return start time of experiment.
Returns: The start time of the experiment.

get_var
(self, name)¶ Return a numpy array from FMU.
Parameters: name – Variable name. Returns: The value of the variable.

get_var_count
(self)¶ Return the number of variables (internal FMU and user declared).
Returns: The number of variables supported by the FMU.

get_var_name
(self, i)¶ Returns the name of a variable.
Parameters: i – Index in ordered dictionary returned by FMUmethod get_model_variables. Returns: The name of the variable.

get_var_rank
(self, name)¶ Not implemented

get_var_shape
(self, name)¶ Not implemented

get_var_type
(self, name)¶ Return type string, compatible with numpy.
Parameters: name – Variable name. Returns: The type of the variable.

get_variables
(self)¶ Return all variables of FMU (both internal and user defined)
Returns: A list of all variables supported by the FMU.

initialize
(self, config_file=None)¶ Initialize FMU with default values
Parameters: config_file – Path to an initialization file.

post
(self)¶ Any postprocessing takes place here.

pre
(self)¶ Any preprocessing takes place here.

reset
(self)¶ Reset the FMU.

setup_experiment
(self, start, stop, dt=1, tol=None)¶ Create an experiment.
Parameters:  start – Start time for the simulation.
 stop – Final time for the simulation.
 dt – Time step size.
 tol – Tolerance of the underlying FMU method.

simulate
(self)¶ Run model from start_time to end_time.

update
(self, dt)¶ Performs one timestep.
The method
setup_experiment
must have been called before.Parameters: dt – Time step size.


rtctools.util.
run_simulation_problem
(simulation_problem_class, base_folder='..', log_level=logging.INFO)¶ Sets up and runs a simulation problem.
Parameters:  simulation_problem_class – Optimization problem class to solve.
 base_folder – Folder within which subfolders “input”, “output”, and “model” exist, containing input and output data, and the model, respectively.
 log_level – The log level to use.
CSV I/O¶

class
rtctools.simulation.csv_mixin.
CSVMixin
¶ Bases:
rtctools.simulation.simulation_problem.SimulationProblem
Adds reading and writing of CSV timeseries and parameters to your simulation problem.
During preprocessing, files named
timeseries_import.csv
,initial_state.csv
, andparameters.csv
are read from theinput
subfolder.During postprocessing, a file named
timeseries_export.csv
is written to theoutput
subfolder.In ensemble mode, a file named
ensemble.csv
is read from theinput
folder. This file contains two columns. The first column gives the name of the ensemble member, and the second column its probability. Furthermore, the other XML files appear one level deeper inside the filesystem hierarchy, inside subfolders with the names of the ensemble members.Variables:  csv_delimiter – Column delimiter used in CSV files. Default is
,
.  csv_validate_timeseries – Check consistency of timeseries. Default is
True
.

timeseries_at
(self, variable, t)¶ Return the value of a time series at the given time.
Parameters:  variable – Variable name.
 t – Time.
Returns: The interpolated value of the time series.
Raises: KeyError
 csv_delimiter – Column delimiter used in CSV files. Default is
DelftFEWS I/O¶

class
rtctools.simulation.pi_mixin.
PIMixin
¶ Bases:
rtctools.simulation.simulation_problem.SimulationProblem
Adds DelftFEWS Published Interface I/O to your simulation problem.
During preprocessing, files named
rtcDataConfig.xml
,timeseries_import.xml
, and``rtcParameterConfig.xml`` are read from theinput
subfolder.rtcDataConfig.xml
maps tuples of FEWS identifiers, including location and parameter ID, to RTCTools time series identifiers.During postprocessing, a file named
timeseries_export.xml
is written to theoutput
subfolder.Variables:  pi_binary_timeseries – Whether to use PI binary timeseries format. Default is
False
.  pi_parameter_config_basenames – List of parameter config file basenames to read. Default is [
rtcParameterConfig
].  pi_check_for_duplicate_parameters – Check if duplicate parameters are read. Default is
False
.  pi_validate_timeseries – Check consistency of timeseries. Default is
True
.

timeseries_at
(self, variable, t)¶ Return the value of a time series at the given time.
Parameters:  variable – Variable name.
 t – Time.
Returns: The interpolated value of the time series.
Raises: KeyError
 pi_binary_timeseries – Whether to use PI binary timeseries format. Default is
Examples¶
This section provides examples demonstrating key features of RTCTools.
Optimization examples¶
This section provides examples demonstrating key features of RTCTools optimization.
Filling a Reservoir¶
Overview¶
The purpose of this example is to understand the technical setup of an RTCTools model, how to run the model, and how to interpret the results.
The scenario is the following: A reservoir operator is trying to fill a reservoir. They are given a sixday forecast of inflows given in 12hour increments. The operator wants to save as much of the inflows as possible, but does not want to end up with too much water at the end of the six days. They have chosen to use RTCTools to calculate how much water to release and when to release it.
The folder <installation directory>\RTCTools2\examples\basic
contains a complete RTCTools optimization problem. An RTCTools
directory has the following structure:
input
: This folder contains the model input data. These are several files in comma separated value format,csv
.model
: This folder contains the Modelica model. The Modelica model contains the physics of the RTCTools model.output
: The folder where the output is saved in the filetimeseries_export.csv
.src
: This folder contains a Python file. This file contains the configuration of the model and is used to run the model .
The Model¶
The first step is to develop a physical model of the system. The model can be viewed and edited using the OpenModelica Connection Editor (OMEdit) program. For how to download and start up OMEdit, see Getting OMEdit.
Make sure to load the Deltares library before loading the example:
 Load the Deltares library into OMEdit
 Using the menu bar: File > Open Model/Library File(s)
 Select
<installation directory>\RTCTools2\mo\Deltares\package.mo
 Load the example model into OMEdit
 Using the menu bar: File > Open Model/Library File(s)
 Select
<installation directory>\RTCTools2\examples\basic\model\Example.mo
Once loaded, we have an OpenModelica Connection Editor window that looks like this:
The model Example.mo
represents a simple system with the following
elements:
 a reservoir, modeled as storage element
Deltares.ChannelFlow.SimpleRouting.Storage.Storage
,  an inflow boundary condition
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Inflow
,  an outfall boundary condition
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Terminal
,  connectors (black lines) connecting the elements.
You can use the mouseover feature help to identify the predefined models from the Deltares library. You can also drag the elements around the connectors will move with the elements. Adding new elements is easy just drag them in from the Deltares Library on the sidebar. Connecting the elements is just as easy click and drag between the ports on the elements.
In text mode, the Modelica model looks as follows (with annotation statements removed):
1 2 3 4 5 6 7 8 9 10 11 12 13 14  model Example
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Inflow inflow;
Deltares.ChannelFlow.SimpleRouting.Storage.Storage storage(V(nominal=4e5, min=2e5, max=6e5));
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Terminal outfall;
input Modelica.SIunits.VolumeFlowRate Q_in(fixed = true);
input Modelica.SIunits.VolumeFlowRate Q_release(fixed = false, min = 0.0, max = 6.5);
output Modelica.SIunits.Volume V_storage;
equation
connect(inflow.QOut, storage.QIn);
connect(storage.QOut, outfall.QIn);
storage.Q_release = Q_release;
inflow.Q = Q_in;
V_storage = storage.V;
end Example;

The three water system elements (storage, inflow, and outfall) appear under
the model Example
statement. The equation
part connects these three
elements with the help of connections. Note that storage
extends the partial
model QSISO
which contains the connectors QIn
and QOut
.
With QSISO
, storage
can be connected on two sides. The storage
element also has a variable Q_release
, which is the decision variable the
operator controls.
OpenModelica Connection Editor will automatically generate the element and
connector entries in the text text file. Defining inputs and outputs requires
editing the text file directly. Relationships between the inputs and outputs and
the library elements must also be defined in the equation
section.
In addition to elements, the input variables Q_in
and Q_release
are also
defined. Q_in
is determined by the forecast and the operator cannot control
it, so we set Q_in(fixed = true)
. The actual values of Q_in
are stored
in timeseries_import.csv
. In the equation
section, equations are defined
to relate the inputs to the appropriate water system elements.
Because we want to view the water volume in the storage element in the output
file, we also define an output variable V_storage
.
The Optimization Problem¶
The python script is created and edited in a text editor. In general, the python script consists of the following blocks:
 Import of packages
 Definition of the optimization problem class
 Constructor
 Objective function
 Definition of constraints
 Any additional configuration
 A run statement
Importing Packages¶
Packages are imported using from ... import ...
at the top of the file. In
our script, we import the classes we want the class to inherit, the
package run_optimization_problem
form the rtctools.util
package, and
any extra packages we want to use. For this example, the import block looks
like:
1 2 3 4 5  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.util import run_optimization_problem

Optimization Problem¶
The next step is to define the optimization problem class. We construct the class by declaring the class and inheriting the desired parent classes. The parent classes each perform different tasks related to importing and exporting data and solving the optimization problem. Each imported class makes a set of methods available to the our optimization class.
8  class Example(CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):

Next, we define an objective function. This is a class method that returns the value that needs to be minimized.
12 13 14 15 16  def objective(self, ensemble_member):
# Minimize water pumped. The total water pumped is the integral of the
# water pumped from the starting time until the stoping time. In
# practice, self.integral() is a summation of all the discrete states.
return self.integral('Q_release', ensemble_member)

Constraints can be declared by declaring the path_constraints()
method. Path
constraints are constraints that are applied every timestep. To set a constraint
at an individual timestep, we could define it inside the constraints()
method.
Other parent classes also declare this method, so we call the super()
method
so that we don’t overwrite their behaviour.
18 19 20 21 22 23  def path_constraints(self, ensemble_member):
# Call super() class to not overwrite default behaviour
constraints = super(Example, self).path_constraints(ensemble_member)
# Constrain the volume of storage between 380000 and 420000 m^3
constraints.append((self.state('storage.V'), 380000, 420000))
return constraints

Run the Optimization Problem¶
To make our script run, at the bottom of our file we just have to call
the run_optimization_problem()
method we imported on the optimization
problem class we just created.
26  run_optimization_problem(Example)

The Whole Script¶
All together, the whole example script is as follows:
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  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.util import run_optimization_problem
class Example(CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):
"""
A basic example for introducing users to RTCTools 2
"""
def objective(self, ensemble_member):
# Minimize water pumped. The total water pumped is the integral of the
# water pumped from the starting time until the stoping time. In
# practice, self.integral() is a summation of all the discrete states.
return self.integral('Q_release', ensemble_member)
def path_constraints(self, ensemble_member):
# Call super() class to not overwrite default behaviour
constraints = super(Example, self).path_constraints(ensemble_member)
# Constrain the volume of storage between 380000 and 420000 m^3
constraints.append((self.state('storage.V'), 380000, 420000))
return constraints
# Run
run_optimization_problem(Example)

Running RTCTools¶
To run this basic example in RTCTools, navigate to the basic example src
directory in the RTCTools shell and run the example using python
example.py
. For more details about using RTCTools, see
Running RTCTools.
Extracting Results¶
The results from the run are found in output\timeseries_export.csv
. Any
CSVreading software can import it, but this is what the results look like when
plotted in Microsoft Excel:
This plot shows that the operator is able to keep the water level within the bounds over the entire time horizon and end with a full reservoir.
Feel free to experiment with this example. See what happens if you change the
max of Q_release
(in the Modelica file) or if you make the objective
function negative (in the python script).
Mixed Integer Optimization: Pumps and Orifices¶
Note
This example focuses on how to incorporate mixed integer components into a hydraulic model, and assumes basic exposure to RTCTools. To start with basics, see Filling a Reservoir.
The Model¶
For this example, the model represents a typical setup for the dewatering of lowland areas. Water is routed from the hinterland (modeled as discharge boundary condition, right side) through a canal (modeled as storage element) towards the sea (modeled as water level boundary condition on the left side). Keeping the lowland area dry requires that enough water is discharged to the sea. If the sea water level is lower than the water level in the canal, the water can be discharged to the sea via gradient flow through the orifice (or a weir). If the sea water level is higher than in the canal, water must be pumped.
To discharge water via gradient flow is free, while pumping costs money. The control task is to keep the water level in the canal below a given flood warning level at minimum costs. The expected result is that the model computes a control pattern that makes use of gradient flow whenever possible and activates the pump only when necessary.
The model can be viewed and edited using the OpenModelica Connection Editor
program. First load the Deltares library into OpenModelica Connection Editor,
and then load the example model, located at
RTCTools2\examples\mixed_integer\model\Example.mo
. The model Example.mo
represents a simple water system with the following elements:
 a canal segment, modeled as storage element
Deltares.ChannelFlow.Hydraulic.Storage.Linear
,  a discharge boundary condition
Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Discharge
,  a water level boundary condition
Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Level
,  a pump
Deltares.ChannelFlow.Hydraulic.Structures.Pump
 an orifice modeled as a pump
Deltares.ChannelFlow.Hydraulic.Structures.Pump
In text mode, the Modelica model looks as follows (with annotation statements removed):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  model Example
// Declare Model Elements
Deltares.ChannelFlow.Hydraulic.Storage.Linear storage(A=1.0e6, H_b=0.0, HQ.H(min=0.0, max=0.5));
Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Discharge discharge;
Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Level level;
Deltares.ChannelFlow.Hydraulic.Structures.Pump pump;
Deltares.ChannelFlow.Hydraulic.Structures.Pump orifice;
// Define Input/Output Variables and set them equal to model variables
input Modelica.SIunits.VolumeFlowRate Q_pump(fixed=false, min=0.0, max=7.0) = pump.Q;
input Boolean is_downhill;
input Modelica.SIunits.VolumeFlowRate Q_in(fixed=true) = discharge.Q;
input Modelica.SIunits.Position H_sea(fixed=true) = level.H;
input Modelica.SIunits.VolumeFlowRate Q_orifice(fixed=false, min=0.0, max=10.0) = orifice.Q;
output Modelica.SIunits.Position storage_level = storage.HQ.H;
output Modelica.SIunits.Position sea_level = level.H;
equation
// Connect Model Elements
connect(orifice.HQDown, level.HQ);
connect(storage.HQ, orifice.HQUp);
connect(storage.HQ, pump.HQUp);
connect(discharge.HQ, storage.HQ);
connect(pump.HQDown, level.HQ);
end Example;

The five water system elements (storage, discharge boundary condition, water
level boundary condition, pump, and orifice) appear under the model Example
statement. The equation
part connects these five elements with the help of
connections. Note that Pump
extends the partial model HQTwoPort
which
inherits from the connector HQPort
. With HQTwoPort
, Pump
can be
connected on two sides. level
represents a model boundary condition (model
is meant in a hydraulic sense here), so it can be connected to one other
element only. It extends the HQOnePort
which again inherits from the
connector HQPort
.
In addition to elements, the input variables Q_in
, H_sea
, Q_pump
,
and Q_orifice
are also defined. Because we want to view the water levels in
the storage element in the output file, we also define output
variables storage_level
and sea_level
. It is usually easiest to set input
and output variables equal to their corresponding model variable in the same line.
To maintain the linearity of the model, we input the Boolean is_downhill
as
a way to keep track of whether water can flow by gravity to the sea. This
variable is not used directly in the hydraulics, but we use it later in the
constraints in the python file.
The Optimization Problem¶
The python script consists of the following blocks:
 Import of packages
 Definition of the optimization problem class
 Constructor
 Objective function
 Definition of constraints
 Additional configuration of the solver
 A run statement
Importing Packages¶
For this example, the import block is as follows:
1 2 3 4 5 6  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.util import run_optimization_problem
from numpy import inf

Note that we are also importing inf
from numpy
. We will use this later
in the constraints.
Optimization Problem¶
Next, we construct the class by declaring it and inheriting the desired parent classes.
9  class Example(CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):

Now we define an objective function. This is a class method that returns the value that needs to be minimized. Here we specify that we want to minimize the volume pumped:
17 18 19 20 21  def objective(self, ensemble_member):
# Minimize water pumped. The total water pumped is the integral of the
# water pumped from the starting time until the stoping time. In
# practice, self.integral() is a summation of all the discrete states.
return self.integral('Q_pump', ensemble_member)

Constraints can be declared by declaring the path_constraints()
method.
Path constraints are constraints that are applied every timestep. To set a
constraint at an individual timestep, define it inside the constraints
method.
The orifice BooleanSubmergedOrifice
requires special constraints to be set
in order to work. They are implemented below in the path_constraints()
method. their parent classes also declare this method, so we call the
super()
method so that we don’t overwrite their behaviour.
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58  def path_constraints(self, ensemble_member):
# Call super to get default constraints
constraints = super(Example, self).path_constraints(ensemble_member)
# M is a handy big number
M = 1e10
# Release through orifice downhill only. This constraint enforces the
# fact that water only flows downhill.
constraints.append(
(self.state('Q_orifice') + (1  self.state('is_downhill')) * 10,
0.0, 10.0))
# Make sure is_downhill is true only when the sea is lower than the
# water level in the storage.
constraints.append((self.state('H_sea')  self.state('storage.HQ.H') 
(1  self.state('is_downhill')) * M, inf, 0.0))
constraints.append((self.state('H_sea')  self.state('storage.HQ.H') +
self.state('is_downhill') * M, 0.0, inf))
# Orifice flow constraint. Uses the equation:
# Q(HUp, HDown, d) = width * C * d * (2 * g * (HUp  HDown)) ^ 0.5
# Note that this equation is only valid for orifices that are submerged
# units: description:
w = 3.0 # m width of orifice
d = 0.8 # m hight of orifice
C = 1.0 # none orifice constant
g = 9.8 # m/s^2 gravitational acceleration
constraints.append(
(((self.state('Q_orifice') / (w * C * d)) ** 2) / (2 * g) +
self.state('orifice.HQDown.H')  self.state('orifice.HQUp.H') 
M * (1  self.state('is_downhill')),
inf, 0.0))
return constraints

Finally, we want to apply some additional configuration, reducing the amount of information the solver outputs:
61 62 63 64 65  def solver_options(self):
options = super(Example, self).solver_options()
# Restrict solver output
options['print_level'] = 1
return options

Run the Optimization Problem¶
To make our script run, at the bottom of our file we just have to call
the run_optimization_problem()
method we imported on the optimization
problem class we just created.
68  run_optimization_problem(Example)

The Whole Script¶
All together, the whole example script is as follows:
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 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.util import run_optimization_problem
from numpy import inf
class Example(CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):
"""
This class is the optimization problem for the Example. Within this class,
the objective, constraints and other options are defined.
"""
# This is a method that returns an expression for the objective function.
# RTCTools always minimizes the objective.
def objective(self, ensemble_member):
# Minimize water pumped. The total water pumped is the integral of the
# water pumped from the starting time until the stoping time. In
# practice, self.integral() is a summation of all the discrete states.
return self.integral('Q_pump', ensemble_member)
# A path constraint is a constraint where the values in the constraint are a
# Timeseries rather than a single number.
def path_constraints(self, ensemble_member):
# Call super to get default constraints
constraints = super(Example, self).path_constraints(ensemble_member)
# M is a handy big number
M = 1e10
# Release through orifice downhill only. This constraint enforces the
# fact that water only flows downhill.
constraints.append(
(self.state('Q_orifice') + (1  self.state('is_downhill')) * 10,
0.0, 10.0))
# Make sure is_downhill is true only when the sea is lower than the
# water level in the storage.
constraints.append((self.state('H_sea')  self.state('storage.HQ.H') 
(1  self.state('is_downhill')) * M, inf, 0.0))
constraints.append((self.state('H_sea')  self.state('storage.HQ.H') +
self.state('is_downhill') * M, 0.0, inf))
# Orifice flow constraint. Uses the equation:
# Q(HUp, HDown, d) = width * C * d * (2 * g * (HUp  HDown)) ^ 0.5
# Note that this equation is only valid for orifices that are submerged
# units: description:
w = 3.0 # m width of orifice
d = 0.8 # m hight of orifice
C = 1.0 # none orifice constant
g = 9.8 # m/s^2 gravitational acceleration
constraints.append(
(((self.state('Q_orifice') / (w * C * d)) ** 2) / (2 * g) +
self.state('orifice.HQDown.H')  self.state('orifice.HQUp.H') 
M * (1  self.state('is_downhill')),
inf, 0.0))
return constraints
# Any solver options can be set here
def solver_options(self):
options = super(Example, self).solver_options()
# Restrict solver output
options['print_level'] = 1
return options
# Run
run_optimization_problem(Example)

Running the Optimization Problem¶
Note
An explaination of bonmin behaviour and output goes here.
Extracting Results¶
The results from the run are found in output/timeseries_export.csv
. Any
CSVreading software can import it, but this is how results can be plotted using
the python library matplotlib:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime
data_path = '../../../examples/mixed_integer/output/timeseries_export.csv'
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)[1:]
# Generate Plot
f, axarr = plt.subplots(2, sharex=True)
axarr[0].set_title('Water Level and Discharge')
# 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(u'Flow Rate [m\u00B3/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()
# 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].legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False)
# Output Plot
plt.show()
(Source code, png, hires.png, pdf)
Observations¶
Note that in the results plotted above, the pump runs with a constantly varying throughput. To smooth out the flow through the pump, consider using goal programming to apply a path goal minimizing the derivative of the pump at each timestep. For an example, see the third goal in Declaring Goals.
Goal Programming: Defining Multiple Objectives¶
Note
This example focuses on how to implement multiobjective optimization in RTCTools using Goal Programming. It assumes basic exposure to RTCTools. If you are a firsttime user of RTCTools, see Filling a Reservoir.
Goal programming is a way to satisfy (sometimes conflicting) goals by ranking the goals by priority. The optimization algorithm will attempt to optimize each goal one at a time, starting with the goal with the highest priority and moving down through the list. Even if a goal cannot be satisfied, the goal programming algorithm will move on when it has found the best possible answer. Goals can be roughly divided into two types:
 As long as we satisfy the goal, we do not care by how much. If we cannot satisfy a goal, any lower priority goals are not allowed to increase the amount by which we exceed (which is equivalent to not allowing any change at all to the exceedance).
 We try to achieve as low a value as possible. Any lower priority goals are not allowed to result in an increase of this value (which is equivalent to not allowing any change at all).
In this example, we will be specifying two goals, on for each type. The higher priority goal will be to maintain the water level of the storage element between two levels. The lower priority goal will be to minimize the total volume pumped.
The Model¶
Note
This example uses the same hydraulic model as the MILP example. For a detailed explanation of the hydraulic model, including how to to formulate mixed integers in your model, see Mixed Integer Optimization: Pumps and Orifices.
For this example, the model represents a typical setup for the dewatering of lowland areas. Water is routed from the hinterland (modeled as discharge boundary condition, right side) through a canal (modeled as storage element) towards the sea (modeled as water level boundary condition on the left side). Keeping the lowland area dry requires that enough water is discharged to the sea. If the sea water level is lower than the water level in the canal, the water can be discharged to the sea via gradient flow through the orifice (or a weir). If the sea water level is higher than in the canal, water must be pumped.
In OpenModelica Connection Editor, the model looks like this:
In text mode, the Modelica model looks as follows (with annotation statements removed):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  model Example
// Declare Model Elements
Deltares.ChannelFlow.Hydraulic.Storage.Linear storage(A=1.0e6, H_b=0.0, HQ.H(min=0.0, max=0.5));
Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Discharge discharge;
Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Level level;
Deltares.ChannelFlow.Hydraulic.Structures.Pump pump;
Deltares.ChannelFlow.Hydraulic.Structures.Pump orifice;
// Define Input/Output Variables and set them equal to model variables
input Modelica.SIunits.VolumeFlowRate Q_pump(fixed=false, min=0.0, max=7.0) = pump.Q;
input Boolean is_downhill;
input Modelica.SIunits.VolumeFlowRate Q_in(fixed=true) = discharge.Q;
input Modelica.SIunits.Position H_sea(fixed=true) = level.H;
input Modelica.SIunits.VolumeFlowRate Q_orifice(fixed=false, min=0.0, max=10.0) = orifice.Q;
output Modelica.SIunits.Position storage_level = storage.HQ.H;
output Modelica.SIunits.Position sea_level = level.H;
equation
// Connect Model Elements
connect(orifice.HQDown, level.HQ);
connect(storage.HQ, orifice.HQUp);
connect(storage.HQ, pump.HQUp);
connect(discharge.HQ, storage.HQ);
connect(pump.HQDown, level.HQ);
end Example;

The Optimization Problem¶
When using goal programming, the python script consists of the following blocks:
 Import of packages
 Declaration of Goals
 Declaration of the optimization problem class
 Constructor
 Declaration of constraint methods
 Specification of Goals
 Declaration of a
priority_completed()
method  Declaration of a
pre()
method  Declaration of a
post()
method  Additional configuration of the solver
 A run statement
Importing Packages¶
For this example, the import block is as follows:
1 2 3 4 5 6 7 8  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.goal_programming_mixin \
import GoalProgrammingMixin, Goal, StateGoal
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.util import run_optimization_problem
from numpy import inf, diff, absolute

Declaring Goals¶
Goals are defined as classes that inherit the Goal
parent class. The
components of goals can be found in ../optimization/multi_objective. In
this example, we demonstrate three ways to define a goal in RTCTools.
First, we have a high priority goal to keep the water level within a minimum and
maximum. Since we are applying this goal to a specific state (model variable) in
our model at every time step, we can inherit a special helper class to define
this goal, called a StateGoal
:
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25  class WaterLevelRangeGoal(StateGoal):
# Applying a state goal to every time step is easily done by defining a goal
# that inherits StateGoal. StateGoal is a helper class that uses the state
# to determine the function, function range, and function nominal
# automatically.
state = 'storage.HQ.H'
# One goal can introduce a single or two constraints (min and/or max). Our
# target water level range is 0.43  0.44. We might not always be able to
# realize this, but we want to try.
target_min = 0.43
target_max = 0.44
# Because we want to satisfy our water level target first, this has a
# higher priority (=lower number).
priority = 1

We also want to save energy, so we define a goal to minimize the integral of
Q_pump
. This goal has a lower priority than the water level range goal. With
nonpath goals, the function range must be large enough to enclose the integral
of the variable over all the timesteps. This goal does not use a helper class:
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44  class MinimizeQpumpGoal(Goal):
# This goal does not use a helper class, so we have to define the function
# method, range and nominal explicitly. We do not specify a target_min or
# target_max in this class, so the goal programming mixin will try to
# minimize the expression returned by the function method.
def function(self, optimization_problem, ensemble_member):
return optimization_problem.integral('Q_pump')
# Every goal needs a rough (over)estimate (enclosure) of the range of the
# function defined above. The nominal is used to scale the value returned by
# the function method so that the value is on the order of 1.
function_range = (0, 540000.0)
function_nominal = 100.0
# The lower the number returned by this function, the higher the priority.
priority = 2
# The penalty variable is taken to the order'th power.
order = 1

We add a third goal minimizing the changes in``Q_pump``, and give it the least priority. This goal smooths out the operation of the pump so that it changes state as few times as possible. To get an idea of what the pump would have done without this goal, see Mixed Integer: Observations. The order of this goal must be 2, so that it penalizes both positive and negative derivatives. Order of 2 is the default, but we include it here explicitly for the sake of clarity.
47 48 49 50 51 52 53 54 55 56 57 58  class MinimizeChangeInQpumpGoal(Goal):
# To reduce pump power cycles, we add a third goal to minimize changes in
# Q_pump. This will be passed into the optimization problem as a path goal
# because it is an an individual goal that should be applied at every time
# step.
def function(self, optimization_problem, ensemble_member):
return optimization_problem.der('Q_pump')
function_range = (10.0, 10.0)
function_nominal = 5.0
priority = 3
# Default order is 2, but we want to be explicit
order = 2

Optimization Problem¶
Next, we construct the class by declaring it and inheriting the desired parent classes.
61 62  class Example(GoalProgrammingMixin, CSVMixin, ModelicaMixin,
CollocatedIntegratedOptimizationProblem):

Constraints can be declared by declaring the path_constraints()
method.
Path constraints are constraints that are applied every timestep. To set a
constraint at an individual timestep, define it inside the constraints()
method.
The “orifice” requires special constraints to be set in order to work. They are
implemented below in the path_constraints()
method. Other parent classes
also declare this method, so we call the super()
method so that we don’t
overwrite their behaviour.
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99  def path_constraints(self, ensemble_member):
# We want to add a few hard constraints to our problem. The goal
# programming mixin however also generates constraints (and objectives)
# from on our goals, so we have to call super() here.
constraints = super(Example, self).path_constraints(ensemble_member)
# Release through orifice downhill only. This constraint enforces the
# fact that water only flows downhill
constraints.append((self.state('Q_orifice') +
(1  self.state('is_downhill')) * 10, 0.0, 10.0))
# Make sure is_downhill is true only when the sea is lower than the
# water level in the storage.
M = 1e10 # M is a handy big number
constraints.append((self.state('H_sea')  self.state('storage.HQ.H') 
(1  self.state('is_downhill')) * M, inf, 0.0))
constraints.append((self.state('H_sea')  self.state('storage.HQ.H') +
self.state('is_downhill') * M, 0.0, inf))
# Orifice flow constraint. Uses the equation:
# Q(HUp, HDown, d) = width * C * d * (2 * g * (HUp  HDown)) ^ 0.5
# Note that this equation is only valid for orifices that are submerged
# units: description:
w = 3.0 # m width of orifice
d = 0.8 # m hight of orifice
C = 1.0 # none orifice constant
g = 9.8 # m/s^2 gravitational acceleration
constraints.append(
(((self.state('Q_orifice') / (w * C * d)) ** 2) / (2 * g) +
self.state('orifice.HQDown.H')  self.state('orifice.HQUp.H') 
M * (1  self.state('is_downhill')),
inf, 0.0))
return constraints

Now we pass in the goals. There are path goals and normal goals, so we have to pass them in using separate methods. A path goal is a specific kind of goal that applies to a particular variable at an individual time step, but that we want to set for all the timesteps.
Nonpath goals are more general goals that are not iteratively applied at every
timestep. We use the goals()
method to pass a list of these goals to the
optimizer.
101 102  def goals(self):
return [MinimizeQpumpGoal()]

For the goals that want to apply our goals to every timestep, so we use the
path_goals()
method. This is a method that returns a list of the path goals
we defined above. Note that with path goals, each timestep is implemented as an
independant goal if we cannot satisfy our min/max on time step A, it will not
affect our desire to satisfy the goal at time step B. Goals that inherit
StateGoal
are always path goals and must always be initialized with the
parameter self
.
104 105 106 107  def path_goals(self):
# Sorting goals on priority is done in the goal programming mixin. We
# do not have to worry about order here.
return [WaterLevelRangeGoal(self), MinimizeChangeInQpumpGoal()]

If all we cared about were the results, we could end our class declaration here. However, it is usually helpful to track how the solution changes after optimizing each priority level. To track these changes, we need to add three methods.
The method pre()
is already defined in RTCTools, but we would like to add
a line to it to create a variable for storing intermediate results. To do this,
we declare a new pre()
method, call super(Example, self).pre()
to ensure
that the original method runs unmodified, and add in a variable declaration to
store our list of intermediate results:
109 110 111 112 113 114  def pre(self):
# Call super() class to not overwrite default behaviour
super(Example, self).pre()
# We keep track of our intermediate results, so that we can print some
# information about the progress of goals at the end of our run.
self.intermediate_results = []

Next, we define the priority_completed()
method to inspect and summarize the
results. These are appended to our intermediate results variable after each
priority is completed.
116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133  def priority_completed(self, priority):
# We want to show that the results of our highest priority goal (water
# level) are remembered. The other information we want to see is how our
# lower priority goal (Q_pump) progresses. We can write some code that
# sumerizes the results and stores it.
# A little bit of tolerance when checking for acceptance, because
# strictly speaking 0.4299... is smaller than 0.43.
_min = 0.43  1e4
_max = 0.44 + 1e4
results = self.extract_results()
n_level_satisfied = sum(
[1 for x in results['storage.HQ.H'] if _min <= x <= _max])
q_pump_integral = sum(results['Q_pump'])
q_pump_sum_changes = sum(absolute(diff(results['Q_pump'])))
self.intermediate_results.append(
(priority, n_level_satisfied, q_pump_integral, q_pump_sum_changes))

We want some way to output our intermediate results. This is accomplished using
the post()
method. Again, we nedd to call the super()
method to avoid
overwiting the internal method.
135 136 137 138 139 140 141 142 143 144  def post(self):
# Call super() class to not overwrite default behaviour
super(Example, self).post()
for priority, n_level_satisfied, q_pump_integral, q_pump_sum_changes \
in self.intermediate_results:
print('\nAfter finishing goals of priority {}:'.format(priority))
print('Level goal satisfied at {} of {} time steps'.format(
n_level_satisfied, len(self.times())))
print('Integral of Q_pump = {:.2f}'.format(q_pump_integral))
print('Sum of Changes in Q_pump: {:.2f}'.format(q_pump_sum_changes))

Finally, we want to apply some additional configuration, reducing the amount of information the solver outputs:
147 148 149 150  def solver_options(self):
options = super(Example, self).solver_options()
options['print_level'] = 1
return options

Run the Optimization Problem¶
To make our script run, at the bottom of our file we just have to call
the run_optimization_problem()
method we imported on the optimization
problem class we just created.
153  run_optimization_problem(Example)

The Whole Script¶
All together, the whole example script is as follows:
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 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.goal_programming_mixin \
import GoalProgrammingMixin, Goal, StateGoal
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.util import run_optimization_problem
from numpy import inf, diff, absolute
class WaterLevelRangeGoal(StateGoal):
# Applying a state goal to every time step is easily done by defining a goal
# that inherits StateGoal. StateGoal is a helper class that uses the state
# to determine the function, function range, and function nominal
# automatically.
state = 'storage.HQ.H'
# One goal can introduce a single or two constraints (min and/or max). Our
# target water level range is 0.43  0.44. We might not always be able to
# realize this, but we want to try.
target_min = 0.43
target_max = 0.44
# Because we want to satisfy our water level target first, this has a
# higher priority (=lower number).
priority = 1
class MinimizeQpumpGoal(Goal):
# This goal does not use a helper class, so we have to define the function
# method, range and nominal explicitly. We do not specify a target_min or
# target_max in this class, so the goal programming mixin will try to
# minimize the expression returned by the function method.
def function(self, optimization_problem, ensemble_member):
return optimization_problem.integral('Q_pump')
# Every goal needs a rough (over)estimate (enclosure) of the range of the
# function defined above. The nominal is used to scale the value returned by
# the function method so that the value is on the order of 1.
function_range = (0, 540000.0)
function_nominal = 100.0
# The lower the number returned by this function, the higher the priority.
priority = 2
# The penalty variable is taken to the order'th power.
order = 1
class MinimizeChangeInQpumpGoal(Goal):
# To reduce pump power cycles, we add a third goal to minimize changes in
# Q_pump. This will be passed into the optimization problem as a path goal
# because it is an an individual goal that should be applied at every time
# step.
def function(self, optimization_problem, ensemble_member):
return optimization_problem.der('Q_pump')
function_range = (10.0, 10.0)
function_nominal = 5.0
priority = 3
# Default order is 2, but we want to be explicit
order = 2
class Example(GoalProgrammingMixin, CSVMixin, ModelicaMixin,
CollocatedIntegratedOptimizationProblem):
"""
An introductory example to goal programming in RCTTools
"""
def path_constraints(self, ensemble_member):
# We want to add a few hard constraints to our problem. The goal
# programming mixin however also generates constraints (and objectives)
# from on our goals, so we have to call super() here.
constraints = super(Example, self).path_constraints(ensemble_member)
# Release through orifice downhill only. This constraint enforces the
# fact that water only flows downhill
constraints.append((self.state('Q_orifice') +
(1  self.state('is_downhill')) * 10, 0.0, 10.0))
# Make sure is_downhill is true only when the sea is lower than the
# water level in the storage.
M = 1e10 # M is a handy big number
constraints.append((self.state('H_sea')  self.state('storage.HQ.H') 
(1  self.state('is_downhill')) * M, inf, 0.0))
constraints.append((self.state('H_sea')  self.state('storage.HQ.H') +
self.state('is_downhill') * M, 0.0, inf))
# Orifice flow constraint. Uses the equation:
# Q(HUp, HDown, d) = width * C * d * (2 * g * (HUp  HDown)) ^ 0.5
# Note that this equation is only valid for orifices that are submerged
# units: description:
w = 3.0 # m width of orifice
d = 0.8 # m hight of orifice
C = 1.0 # none orifice constant
g = 9.8 # m/s^2 gravitational acceleration
constraints.append(
(((self.state('Q_orifice') / (w * C * d)) ** 2) / (2 * g) +
self.state('orifice.HQDown.H')  self.state('orifice.HQUp.H') 
M * (1  self.state('is_downhill')),
inf, 0.0))
return constraints
def goals(self):
return [MinimizeQpumpGoal()]
def path_goals(self):
# Sorting goals on priority is done in the goal programming mixin. We
# do not have to worry about order here.
return [WaterLevelRangeGoal(self), MinimizeChangeInQpumpGoal()]
def pre(self):
# Call super() class to not overwrite default behaviour
super(Example, self).pre()
# We keep track of our intermediate results, so that we can print some
# information about the progress of goals at the end of our run.
self.intermediate_results = []
def priority_completed(self, priority):
# We want to show that the results of our highest priority goal (water
# level) are remembered. The other information we want to see is how our
# lower priority goal (Q_pump) progresses. We can write some code that
# sumerizes the results and stores it.
# A little bit of tolerance when checking for acceptance, because
# strictly speaking 0.4299... is smaller than 0.43.
_min = 0.43  1e4
_max = 0.44 + 1e4
results = self.extract_results()
n_level_satisfied = sum(
[1 for x in results['storage.HQ.H'] if _min <= x <= _max])
q_pump_integral = sum(results['Q_pump'])
q_pump_sum_changes = sum(absolute(diff(results['Q_pump'])))
self.intermediate_results.append(
(priority, n_level_satisfied, q_pump_integral, q_pump_sum_changes))
def post(self):
# Call super() class to not overwrite default behaviour
super(Example, self).post()
for priority, n_level_satisfied, q_pump_integral, q_pump_sum_changes \
in self.intermediate_results:
print('\nAfter finishing goals of priority {}:'.format(priority))
print('Level goal satisfied at {} of {} time steps'.format(
n_level_satisfied, len(self.times())))
print('Integral of Q_pump = {:.2f}'.format(q_pump_integral))
print('Sum of Changes in Q_pump: {:.2f}'.format(q_pump_sum_changes))
# Any solver options can be set here
def solver_options(self):
options = super(Example, self).solver_options()
options['print_level'] = 1
return options
# Run
run_optimization_problem(Example)

Running the Optimization Problem¶
Following the execution of the optimization problem, the post()
method
should print out the following lines:
After finishing goals of priority 1:
Level goal satisfied at 19 of 21 time steps
Integral of Q_pump = 74.18
Sum of Changes in Q_pump: 7.83
After finishing goals of priority 2:
Level goal satisfied at 19 of 21 time steps
Integral of Q_pump = 60.10
Sum of Changes in Q_pump: 11.70
After finishing goals of priority 3:
Level goal satisfied at 19 of 21 time steps
Integral of Q_pump = 60.10
Sum of Changes in Q_pump: 10.07
As the output indicates, while optimizing for the priority 1 goal, no attempt
was made to minimize the integral of Q_pump
. The only objective was to
minimize the number of states in violation of the water level goal.
After optimizing for the priority 2 goal, the solver was able to find a solution
that reduced the integral of Q_pump
without increasing the number of
timesteps where the water level exceeded the limit. However, this solution
induced additional variation into the operation of Q_pump
After optimizing the priority 3 goal, the integral of Q_pump
is the same and
the level goal has not improved. Without hurting any higher priority goals,
RTCTools was able to smooth out the operation of the pump.
Extracting Results¶
The results from the run are found in output/timeseries_export.csv
. Any
CSVreading software can import it, but this is how results can be plotted using
the python library matplotlib:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime
data_path = '../../../examples/goal_programming/output/timeseries_export.csv'
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)[1:]
# Generate Plot
n_subplots = 3
f, axarr = plt.subplots(n_subplots, sharex=True, figsize=(8, 3 * n_subplots))
axarr[0].set_title('Water Level and Discharge')
# 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')
# Middle subplot
axarr[1].set_ylabel('Water Level [m]')
axarr[1].plot(results['time'], results['storage_level'], label='Storage',
linewidth=2, color='b')
axarr[1].plot(results['time'], 0.44 * np.ones_like(results['time']), label='Storage Max',
linewidth=2, color='r', linestyle='')
axarr[1].plot(results['time'], 0.43 * np.ones_like(results['time']), label='Storage Min',
linewidth=2, color='g', linestyle='')
# Lower Subplot
axarr[2].set_ylabel(u'Flow Rate [m\u00B3/s]')
axarr[2].plot(results['time'], results['Q_orifice'], label='Orifice',
linewidth=2, color='g')
axarr[2].plot(results['time'], results['Q_pump'], label='Pump',
linewidth=2, color='r')
axarr[2].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
f.autofmt_xdate()
# Shrink each axis by 20% and put a legend to the right of the axis
for i in range(n_subplots):
box = axarr[i].get_position()
axarr[i].set_position([box.x0, box.y0, box.width * 0.8, box.height])
axarr[i].legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False)
# Output Plot
plt.show()
(Source code, png, hires.png, pdf)
Using Lookup Tables¶
Note
This example focuses on how to implement nonlinear storage elements in RTCTools using lookup tables. It assumes basic exposure to RTCTools. If you are a firsttime user of RTCTools, see Filling a Reservoir.
This example also uses goal programming in the formulation. If you are unfamiliar with goal programming, please see Goal Programming: Defining Multiple Objectives.
The Model¶
Note
This example uses the same hydraulic model as the basic example. For a detalied explaination of the hydraulic model, see Filling a Reservoir.
In OpenModelica Connection Editor, the model looks like this:
In text mode, the Modelica model is as follows (with annotation statements removed):
1 2 3 4 5 6 7 8 9 10 11 12  model Example
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Inflow inflow;
Deltares.ChannelFlow.SimpleRouting.Storage.Storage storage(V(nominal=4e5, min=2e5, max=6e5));
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Terminal outfall;
input Modelica.SIunits.VolumeFlowRate Q_in(fixed = true);
input Modelica.SIunits.VolumeFlowRate Q_release(fixed = false, min = 0.0, max = 10.0);
equation
connect(inflow.QOut, storage.QIn);
connect(storage.QOut, outfall.QIn);
storage.Q_release = Q_release;
inflow.Q = Q_in;
end Example;

The Optimization Problem¶
The python script consists of the following blocks:
 Import of packages
 Declaration of Goals
 Declaration of the optimization problem class
 Constructor
 Declaration of a
pre()
method  Specification of Goals
 Declaration of a
priority_completed()
method  Declaration of a
post()
method  Additional configuration of the solver
 A run statement
Importing Packages¶
For this example, the import block is as follows:
1 2 3 4 5 6 7 8 9  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.goal_programming_mixin \
import GoalProgrammingMixin, Goal, StateGoal
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.optimization.csv_lookup_table_mixin import CSVLookupTableMixin
from rtctools.util import run_optimization_problem
import numpy as np

Declaring Goals¶
Goals are defined as classes that inherit the Goal
parent class. The
components of goals can be found in ../optimization/multi_objective. In
this example, we use the helper goal class, StateGoal
.
First, we have a high priority goal to keep the water volume within a minimum and maximum. We use a water volume goal instead of a water level goal when the volumestorage relation of the storage element is nonlinear. The volume of water in the storage element behaves linearly, while the water level does not.
However, goals are usually defined in the form of water level goals. We will
convert the water level goals into volume goals within the optimization
problem class, so we define the __init__()
method so we can pass the
values of the goals in later. We call the super()
method to avoid
overwriting the __init__()
method of the parent class.
11 12 13 14 15 16 17 18 19 20 21 22 23 24  class WaterVolumeRangeGoal(StateGoal):
# We want to add a water volume range goal to our optimization. However, at
# the time of defining this goal we still do not know what the value of the
# min and max are. We add an __init__() method so that the values of these
# goals can be defined when the optimization problem class instantiates
# this goal.
def __init__(self, optimization_problem):
# Call super class first, and pass in the optimization problem
super(WaterVolumeRangeGoal, self).__init__(optimization_problem)
# Assign V_min and V_max the the target range
self.target_min = optimization_problem.get_timeseries('V_min')
self.target_max = optimization_problem.get_timeseries('V_max')
state = 'storage.V'
priority = 1

We also want to save energy, so we define a goal to minimize Q_release
. This
goal has a lower priority.
27 28 29 30 31 32 33  class MinimizeQreleaseGoal(StateGoal):
# GoalProgrammingMixin will try to minimize the following state:
state = 'Q_release'
# The lower the number returned by this function, the higher the priority.
priority = 2
# The penalty variable is taken to the order'th power.
order = 1

Optimization Problem¶
Next, we construct the class by declaring it and inheriting the desired parent classes.
36 37  class Example(GoalProgrammingMixin, CSVLookupTableMixin, CSVMixin,
ModelicaMixin, CollocatedIntegratedOptimizationProblem):

The method pre()
is already defined in RTCTools, but we would like to add
a line to it to create a variable for storing intermediate results. To do this,
we declare a new pre()
method, call super(Example, self).pre()
to ensure
that the original method runs unmodified, and add in a variable declaration to
store our list of intermediate results.
We also want to convert our water level rane goal into a water volume range
goal. We can access the spline function describing the water levelstorage
relation using the lookup_table()
method. We cache the functions for
convenience. The lookup_storage_V()
method can convert timeseries objects,
and we save the water volume goal bounds as timeseries.
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63  def pre(self):
super(Example, self).pre()
# Empty list for storing intermediate_results
self.intermediate_results = []
# Cache lookup tables for convenience and legibility
_lookup_tables = self.lookup_tables(ensemble_member=0)
self.lookup_storage_V = _lookup_tables['storage_V']
# Nonvarrying goals can be implemented as a timeseries like this:
self.set_timeseries('H_min', np.ones_like(self.times()) * 0.44, output=False)
# Q_in is a varying input and is defined in timeseries_import.csv
# However, if we set it again here, it will be added to the output file
self.set_timeseries('Q_in', self.get_timeseries('Q_in'))
# Convert our water level constraints into volume constraints
self.set_timeseries('V_max',
self.lookup_storage_V(self.get_timeseries('H_max')))
self.set_timeseries('V_min',
self.lookup_storage_V(self.get_timeseries('H_min')))

Notice that H_max was not defined in pre(). This is because it was defined as a timeseries import. We access timeseries using get_timeseries() and store them using set_timeseries(). Once a timeseries is set, we can access it later. In addition, all timeseries that are set are automatically included in the output file. You can find more information on timeseries here ../optimization/basics.
Now we pass in the goals. We want to apply our goals to every timestep, so we
use the path_goals()
method. This is a method that returns a list of the
goals we defined above. The WaterVolumeRangeGoal
needs to be instantiated
with the new water volume timeseries we just defined.
65 66 67 68 69  def path_goals(self):
g = []
g.append(WaterVolumeRangeGoal(self))
g.append(MinimizeQreleaseGoal(self))
return g

If all we cared about were the results, we could end our class declaration here. However, it is usually helpful to track how the solution changes after optimizing each priority level. To track these changes, we need to add three methods.
We define the priority_completed()
method to inspect and summerize the
results. These are appended to our intermediate results variable after each
priority is completed.
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90  def priority_completed(self, priority):
results = self.extract_results()
self.set_timeseries('storage_V', results['storage.V'])
_max = self.get_timeseries('V_max').values
_min = self.get_timeseries('V_min').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 <= storage_V, storage_V <= _max))
q_release_integral = sum(results['Q_release'])
self.intermediate_results.append(
(priority, n_level_satisfied, q_release_integral))

We output our intermediate results using the post()
method. Again, we nedd
to call the super()
method to avoid overwiting the internal method.
92 93 94 95 96 97 98 99  def post(self):
# Call super() class to not overwrite default behaviour
super(Example, self).post()
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_release_integral))

Finally, we want to apply some additional configuration, reducing the amount of information the solver outputs:
102 103 104 105  def solver_options(self):
options = super(Example, self).solver_options()
options['print_level'] = 1
return options

Run the Optimization Problem¶
To make our script run, at the bottom of our file we just have to call
the run_optimization_problem()
method we imported on the optimization
problem class we just created.
108  run_optimization_problem(Example)

The Whole Script¶
All together, the whole example script is as follows:
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 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.goal_programming_mixin \
import GoalProgrammingMixin, Goal, StateGoal
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.optimization.csv_lookup_table_mixin import CSVLookupTableMixin
from rtctools.util import run_optimization_problem
import numpy as np
class WaterVolumeRangeGoal(StateGoal):
# We want to add a water volume range goal to our optimization. However, at
# the time of defining this goal we still do not know what the value of the
# min and max are. We add an __init__() method so that the values of these
# goals can be defined when the optimization problem class instantiates
# this goal.
def __init__(self, optimization_problem):
# Call super class first, and pass in the optimization problem
super(WaterVolumeRangeGoal, self).__init__(optimization_problem)
# Assign V_min and V_max the the target range
self.target_min = optimization_problem.get_timeseries('V_min')
self.target_max = optimization_problem.get_timeseries('V_max')
state = 'storage.V'
priority = 1
class MinimizeQreleaseGoal(StateGoal):
# GoalProgrammingMixin will try to minimize the following state:
state = 'Q_release'
# The lower the number returned by this function, the higher the priority.
priority = 2
# The penalty variable is taken to the order'th power.
order = 1
class Example(GoalProgrammingMixin, CSVLookupTableMixin, CSVMixin,
ModelicaMixin, CollocatedIntegratedOptimizationProblem):
"""
An extention of the goal programming example that shows how to incorporate
nonlinear storage elements in the model.
"""
def pre(self):
super(Example, self).pre()
# Empty list for storing intermediate_results
self.intermediate_results = []
# Cache lookup tables for convenience and legibility
_lookup_tables = self.lookup_tables(ensemble_member=0)
self.lookup_storage_V = _lookup_tables['storage_V']
# Nonvarrying goals can be implemented as a timeseries like this:
self.set_timeseries('H_min', np.ones_like(self.times()) * 0.44, output=False)
# Q_in is a varying input and is defined in timeseries_import.csv
# However, if we set it again here, it will be added to the output file
self.set_timeseries('Q_in', self.get_timeseries('Q_in'))
# Convert our water level constraints into volume constraints
self.set_timeseries('V_max',
self.lookup_storage_V(self.get_timeseries('H_max')))
self.set_timeseries('V_min',
self.lookup_storage_V(self.get_timeseries('H_min')))
def path_goals(self):
g = []
g.append(WaterVolumeRangeGoal(self))
g.append(MinimizeQreleaseGoal(self))
return g
# 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 (see post() method below).
def priority_completed(self, priority):
results = self.extract_results()
self.set_timeseries('storage_V', results['storage.V'])
_max = self.get_timeseries('V_max').values
_min = self.get_timeseries('V_min').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 <= storage_V, storage_V <= _max))
q_release_integral = sum(results['Q_release'])
self.intermediate_results.append(
(priority, n_level_satisfied, q_release_integral))
def post(self):
# Call super() class to not overwrite default behaviour
super(Example, self).post()
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_release_integral))
# Any solver options can be set here
def solver_options(self):
options = super(Example, self).solver_options()
options['print_level'] = 1
return options
# Run
run_optimization_problem(Example)

Running the Optimization Problem¶
Following the execution of the optimization problem, the post()
method
should print out the following lines:
After finishing goals of priority 1:
Volume goal satisfied at 12 of 12 time steps
Integral of Q_release = 42.69
After finishing goals of priority 2:
Volume goal satisfied at 12 of 12 time steps
Integral of Q_release = 42.58
As the output indicates, while optimizing for the priority 1 goal, no attempt
was made to minimize the integral of Q_release
. The only objective was to
minimize the number of states in violation of the water level goal.
After optimizing for the priority 2 goal, the solver was able to find a solution
that reduced the integral of Q_release
without increasing the number of
timesteps where the water level exceded the limit.
Extracting Results¶
The results from the run are found in output/timeseries_export.csv
. Any
CSVreading software can import it, but this is how results can be plotted using
the python library matplotlib:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime
data_path = '../../../examples/lookup_table/output/timeseries_export.csv'
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)[1:]
# Generate Plot
n_subplots = 2
f, axarr = plt.subplots(n_subplots, sharex=True, figsize=(8, 3 * n_subplots))
axarr[0].set_title('Water Volume and Discharge')
f.autofmt_xdate()
# Upper subplot
axarr[0].set_ylabel(u'Water Volume [m\u00B3]')
axarr[0].ticklabel_format(style='sci', axis='y', scilimits=(0,0))
axarr[0].plot(results['time'], results['storage_V'], label='Storage',
linewidth=2, color='b')
axarr[0].plot(results['time'], results['V_max'], label='Storage Max',
linewidth=2, color='r', linestyle='')
axarr[0].plot(results['time'], results['V_min'], label='Storage Min',
linewidth=2, color='g', linestyle='')
# Lower Subplot
axarr[1].set_ylabel(u'Flow Rate [m\u00B3/s]')
axarr[1].plot(results['time'], results['Q_in'], label='Inflow',
linewidth=2, color='g')
axarr[1].plot(results['time'], results['Q_release'], label='Release',
linewidth=2, color='r')
# Shrink each axis by 20% and put a legend to the right of the axis
for i in range(n_subplots):
box = axarr[i].get_position()
axarr[i].set_position([box.x0, box.y0, box.width * 0.8, box.height])
axarr[i].legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False)
# Output Plot
plt.show()
(Source code, png, hires.png, pdf)
Using an Ensemble Forecast¶
Note
This example is an extension of Using Lookup Tables. It assumes prior knowledge of goal programming and the lookuptables components of RTCTools. If you are a firsttime user of RTCTools, see Filling a Reservoir.
Then biggest change to RTCTools when using an ensemble is the structure of the
directory. The folder RTCTools2\examples\ensemble
contains a complete RTCTools ensemble optimization problem. An RTCTools
ensemble directory has the following structure:
model
: This folder contains the Modelica model. The Modelica model contains the physics of the RTCTools model.src
: This folder contains a Python file. This file contains the configuration of the model and is used to run the model.input
: This folder contains the model input data pertaining to each ensemble member:ensemble.csv
: a file where the names and probabilities of the ensemble members are definedforcast1
 the file
timeseries_import.csv
 the file
initial_state.csv
 the file
forcast2
timeseries_import.csv
initial_state.csv
 ...
output
: The folder where the output is saved:forcast1
timeseries_export.csv
forcast2
timeseries_export.csv
 ...
The Model¶
Note
This example uses the same hydraulic model as the basic example. For a detailed explanation of the hydraulic model, see Filling a Reservoir.
In OpenModelica Connection Editor, the model looks like this:
In text mode, the Modelica model is as follows (with annotation statements removed):
1 2 3 4 5 6 7 8 9 10 11 12  model Example
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Inflow inflow;
Deltares.ChannelFlow.SimpleRouting.Storage.Storage storage(V(nominal=4e5, min=2e5, max=6e5));
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Terminal outfall;
input Modelica.SIunits.VolumeFlowRate Q_in(fixed = true);
input Modelica.SIunits.VolumeFlowRate Q_release(fixed = false, min = 0.0, max = 6.0);
equation
connect(inflow.QOut, storage.QIn);
connect(storage.QOut, outfall.QIn);
storage.Q_release = Q_release;
inflow.Q = Q_in;
end Example;

The Optimization Problem¶
The python script consists of the following blocks:
 Import of packages
 Declaration of Goals
 Declaration of the optimization problem class
 Constructor
 Set
csv_ensemble_mode = True
 Declaration of a
pre()
method  Specification of Goals
 Declaration of a
priority_completed()
method  Declaration of a
post()
method  Additional configuration of the solver
 A run statement
Importing Packages¶
For this example, the import block is as follows:
1 2 3 4 5 6 7 8 9 10  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.goal_programming_mixin \
import GoalProgrammingMixin, Goal, StateGoal
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.optimization.control_tree_mixin import ControlTreeMixin
from rtctools.optimization.csv_lookup_table_mixin import CSVLookupTableMixin
from rtctools.util import run_optimization_problem
import numpy as np

Declaring Goals¶
First, we have a high priority goal to keep the water volume within a minimum and maximum.
13 14 15 16 17 18 19 20 21  class WaterVolumeRangeGoal(StateGoal):
def __init__(self, optimization_problem):
# Call super class first, and pass in the optimization problem
super(WaterVolumeRangeGoal, self).__init__(optimization_problem)
# Assign V_min and V_max the the target range
self.target_min = optimization_problem.get_timeseries('V_min')
self.target_max = optimization_problem.get_timeseries('V_max')
state = 'storage.V'
priority = 1

We also want to save energy, so we define a goal to minimize Q_release
. This
goal has a lower priority.
24 25 26 27 28 29 30  class MinimizeQreleaseGoal(StateGoal):
# GoalProgrammingMixin will try to minimize the following state
state = 'Q_release'
# The lower the number returned by this function, the higher the priority.
priority = 2
# The penalty variable is taken to the order'th power.
order = 1

Optimization Problem¶
Next, we construct the class by declaring it and inheriting the desired parent classes.
33 34  class Example(GoalProgrammingMixin, ControlTreeMixin, CSVLookupTableMixin,
CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):

We turn on ensemble mode by setting csv_ensemble_mode = True
:
49 50  ensemble_member: [] for ensemble_member in range(self.ensemble_size)}

The method pre()
is already defined in RTCTools, but we would like to add
a line to it to create a variable for storing intermediate results. To do
this, we declare a new pre()
method, call super(Example, self).pre()
to ensure that the original method runs unmodified, and add in a variable
declaration to store our list of intermediate results. This variable is a
dict, reflecting the need to store results from multiple ensemble.
Because the timeseries we set will be the same for both ensemble members, we also make sure that the timeseries we set are set for both ensemble members using for loops.
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74  def pre(self):
# Do the standard preprocessing
super(Example, self).pre()
# Create a dict of empty lists for storing intermediate results from
# each ensemble
self.intermediate_results = {
ensemble_member: [] for ensemble_member in range(self.ensemble_size)}
# Cache lookup tables for convenience and code legibility
_lookup_tables = self.lookup_tables(ensemble_member=0)
self.lookup_storage_V = _lookup_tables['storage_V']
# Nonvarying goals can be implemented as a timeseries
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=e_m)
self.set_timeseries('V_min',
self.lookup_storage_V(self.get_timeseries('H_min')),
ensemble_member=e_m)

Now we pass in the goals:
76 77 78 79 80  def path_goals(self):
g = []
g.append(WaterVolumeRangeGoal(self))
g.append(MinimizeQreleaseGoal(self))
return g

In order to better demonstrate the way that ensembles are handled in RTC
Tools, we modify the control_tree_options()
method. The default setting
allows the control tree to split at every time, but we override this method
and force it to split at a single timestep. See Observations at
the bottom of the page for more information.
82 83 84 85 86 87 88  def control_tree_options(self):
# We want to modify the control tree options, so we override the default
# control_tree_options method. We call super() to get the default options
options = super(Example, self).control_tree_options()
# Change the branching_times list to only contain the fifth timestep
options['branching_times'] = [self.times()[5]]
return options

We define the priority_completed()
method. We ensure that it stores the
results from both ensemble members.
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111  def priority_completed(self, priority):
# 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 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 = 10
_max += tol
_min = tol
n_level_satisfied = sum(
np.logical_and(_min <= V_storage, V_storage <= _max))
q_release_integral = sum(results['Q_release'])
self.intermediate_results[e_m].append((priority, n_level_satisfied,
q_release_integral))

We output our intermediate results using the post()
method:
113 114 115 116 117 118 119 120 121 122  def post(self):
super(Example, self).post()
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[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())))
print("Integral of Q_release = {:.2f}".format(q_release_integral))

Finally, we want to apply some additional configuration, reducing the amount of information the solver outputs:
125 126 127 128 129 130 131 132  def solver_options(self):
options = super(Example, self).solver_options()
# When mumps_scaling is not zero, errors occur. RTCTools does its own
# scaling, so mumps scaling is not critical. Proprietary HSL solvers
# do not exhibit this error.
options['mumps_scaling'] = 0
options['print_level'] = 1
return options

Run the Optimization Problem¶
To make our script run, at the bottom of our file we just have to call
the run_optimization_problem()
method we imported on the optimization
problem class we just created.
135  run_optimization_problem(Example)

The Whole Script¶
All together, the whole example script is as follows:
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 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135  from rtctools.optimization.collocated_integrated_optimization_problem \
import CollocatedIntegratedOptimizationProblem
from rtctools.optimization.goal_programming_mixin \
import GoalProgrammingMixin, Goal, StateGoal
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.optimization.control_tree_mixin import ControlTreeMixin
from rtctools.optimization.csv_lookup_table_mixin import CSVLookupTableMixin
from rtctools.util import run_optimization_problem
import numpy as np
class WaterVolumeRangeGoal(StateGoal):
def __init__(self, optimization_problem):
# Call super class first, and pass in the optimization problem
super(WaterVolumeRangeGoal, self).__init__(optimization_problem)
# Assign V_min and V_max the the target range
self.target_min = optimization_problem.get_timeseries('V_min')
self.target_max = optimization_problem.get_timeseries('V_max')
state = 'storage.V'
priority = 1
class MinimizeQreleaseGoal(StateGoal):
# GoalProgrammingMixin will try to minimize the following state
state = 'Q_release'
# The lower the number returned by this function, the higher the priority.
priority = 2
# The penalty variable is taken to the order'th power.
order = 1
class Example(GoalProgrammingMixin, ControlTreeMixin, CSVLookupTableMixin,
CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):
"""
An extention of the goal programming and lookuptable examples that
demonstrates how to work with ensembles.
"""
# Overide default csv_ensemble_mode = False from CSVMixin before calling pre()
csv_ensemble_mode = True
def pre(self):
# Do the standard preprocessing
super(Example, self).pre()
# Create a dict of empty lists for storing intermediate results from
# each ensemble
self.intermediate_results = {
ensemble_member: [] for ensemble_member in range(self.ensemble_size)}
# Cache lookup tables for convenience and code legibility
_lookup_tables = self.lookup_tables(ensemble_member=0)
self.lookup_storage_V = _lookup_tables['storage_V']
# Nonvarying goals can be implemented as a timeseries
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=e_m)
self.set_timeseries('V_min',
self.lookup_storage_V(self.get_timeseries('H_min')),
ensemble_member=e_m)
def path_goals(self):
g = []
g.append(WaterVolumeRangeGoal(self))
g.append(MinimizeQreleaseGoal(self))
return g
def control_tree_options(self):
# We want to modify the control tree options, so we override the default
# control_tree_options method. We call super() to get the default options
options = super(Example, self).control_tree_options()
# Change the branching_times list to only contain the fifth timestep
options['branching_times'] = [self.times()[5]]
return options
def priority_completed(self, priority):
# 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 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 = 10
_max += tol
_min = tol
n_level_satisfied = sum(
np.logical_and(_min <= V_storage, V_storage <= _max))
q_release_integral = sum(results['Q_release'])
self.intermediate_results[e_m].append((priority, n_level_satisfied,
q_release_integral))
def post(self):
super(Example, self).post()
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[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())))
print("Integral of Q_release = {:.2f}".format(q_release_integral))
# Any solver options can be set here
def solver_options(self):
options = super(Example, self).solver_options()
# When mumps_scaling is not zero, errors occur. RTCTools does its own
# scaling, so mumps scaling is not critical. Proprietary HSL solvers
# do not exhibit this error.
options['mumps_scaling'] = 0
options['print_level'] = 1
return options
# Run
run_optimization_problem(Example)

Running the Optimization Problem¶
Following the execution of the optimization problem, the post()
method
should print out the following lines:
Results for Ensemble Member 0:
After finishing goals of priority 1:
Level goal satisfied at 10 of 12 time steps
Integral of Q_release = 17.34
After finishing goals of priority 2:
Level goal satisfied at 9 of 12 time steps
Integral of Q_release = 17.12
Results for Ensemble Member 1:
After finishing goals of priority 1:
Level goal satisfied at 10 of 12 time steps
Integral of Q_release = 20.82
After finishing goals of priority 2:
Level goal satisfied at 9 of 12 time steps
Integral of Q_release = 20.60
This is the same output as the output for Mixed Integer Optimization: Pumps and Orifices, except now the output for each ensemble is printed.
Extracting Results¶
The results from the run are found in output/forcast1/timeseries_export.csv
and output/forcast2/timeseries_export.csv
. Any CSVreading software can
import it, but this is how results can be plotted using the python library
matplotlib:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime
from pylab import get_cmap
forecast_names = ['forecast1', 'forecast2']
# Import Data
def get_results(forecast_name):
output_dir = '../../../examples/ensemble/output/'
data_path = output_dir + forecast_name + '/timeseries_export.csv'
delimiter = ','
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
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(u'Water Volume in Storage [m\u00B3]')
axarr[0].ticklabel_format(style='sci', axis='y', scilimits=(0,0))
# Lower Subplot
axarr[1].set_ylabel(u'Flow Rate [m\u00B3/s]')
# Plot Ensemble Members
for idx, forecast in enumerate(forecast_names):
# Upper Subplot
results = get_results(forecast)
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 + ':Volume',
linewidth=2, color=get_cmap(cmaps[idx])(shades[1]))
# Lower Subplot
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 30% 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.7, box.height])
axarr[i].legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False)
# Output Plot
plt.show()
(Source code, png, hires.png, pdf)
Observations¶
Note that in the results plotted above, the control tree follows a single path and does not branch until it arrives at the first branching time. Up until the branching point, RTCTools is making decisions that enhance the flexibility of the system so that it can respond as ideally as possible to whichever future emerges. In the case of two forecasts, this means taking a path that falls between the two possible futures. This will cause the water level to diverge from the ideal levels as time progresses. While this appears to be suboptimal, it is preferable to simply gambling on one of the forecasts coming true and ignoring the other. Once the branching time is reached, RTCTools is allowed to optimize for each individual branch separately. Immediately, RTCTools applies the corrective control needed to get the water levels into the acceptable range. If the operator simply picks a forecast to use and guesses wrong, the corrective control will have to be much more drastic and potentially catastrophic.
Simulation examples¶
This section provides examples demonstrating key features of RTCTools simulation.
Tracking a Setpoint¶
Overview¶
The purpose of this example is to understand the technical setup of an RTC Tools simulation model, how to run the model, and how to access the results.
The scenario is the following: A reservoir operator is trying to keep the reservoir’s volume close to a given target volume. They are given a sixday forecast of inflows given in 12hour increments. To keep things simple, we ignore the waterlevelstorage relation of the reservoir and headdischarge relationships in this example. To make things interesting, the reservoir operator is only able to release water at a few discrete flow rates, and only change the discrete flow rate every 12 hours. They have chosen to use the RTC Tools simulator to see if a simple proportional controller will be able to keep the system close enough to the target water volume.
The folder <installation directory>\RTCTools2\examples\simulation
contains a complete RTCTools simulation problem. An RTCTools
directory has the following structure:
input
: This folder contains the model input data. These are several files in comma separated value format,csv
.model
: This folder contains the Modelica model. The Modelica model contains the physics of the RTCTools model.output
: The folder where the output is saved in the filetimeseries_export.csv
.src
: This folder contains a Python file. This file contains the configuration of the model and is used to run the model .
The Model¶
The first step is to develop a physical model of the system. The model can be viewed and edited using the OpenModelica Connection Editor (OMEdit) program. For how to download and start up OMEdit, see Getting OMEdit.
Make sure to load the Deltares library before loading the example:
 Load the Deltares library into OMEdit
 Using the menu bar: File > Open Model/Library File(s)
 Select
<installation directory>\RTCTools2\mo\Deltares\package.mo
 Load the example model into OMEdit
 Using the menu bar: File > Open Model/Library File(s)
 Select
<installation directory>\RTCTools2\examples\simulation\model\Example.mo
Once loaded, we have an OpenModelica Connection Editor window that looks like this:
The model Example.mo
represents a simple system with the following
elements:
 a reservoir, modeled as storage element
Deltares.ChannelFlow.SimpleRouting.Storage.Storage
,  an inflow boundary condition
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Inflow
,  an outfall boundary condition
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Terminal
,  connectors (black lines) connecting the elements.
You can use the mouseover feature help to identify the predefined models from the Deltares library. You can also drag the elements around the connectors will move with the elements. Adding new elements is easy just drag them in from the Deltares Library on the sidebar. Connecting the elements is just as easy click and drag between the ports on the elements.
In text mode, the Modelica model looks as follows (with annotation statements removed):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18  model Example
// Elements
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Inflow inflow(Q = Q_in);
Deltares.ChannelFlow.SimpleRouting.Storage.Storage storage(Q_release = P_control, V(start=storage_V_init, nominal=4e5));
Deltares.ChannelFlow.SimpleRouting.BoundaryConditions.Terminal outfall;
// Initial States
parameter Modelica.SIunits.Volume storage_V_init;
// Inputs
input Modelica.SIunits.VolumeFlowRate P_control(fixed = true);
input Modelica.SIunits.VolumeFlowRate Q_in(fixed = true);
input Modelica.SIunits.VolumeFlowRate storage_V_target(fixed = true);
// Outputs
output Modelica.SIunits.Volume storage_V = storage.V;
output Modelica.SIunits.VolumeFlowRate Q_release = P_control;
equation
connect(inflow.QOut, storage.QIn);
connect(storage.QOut, outfall.QIn);
end Example;

The three water system elements (storage, inflow, and outfall) appear under
the model Example
statement. The equation
part connects these three
elements with the help of connections. Note that storage
extends the partial
model QSISO
which contains the connectors QIn
and QOut
.
With QSISO
, storage
can be connected on two sides. The storage
element also has a variable Q_release
, which is the decision variable the
operator controls.
OpenModelica Connection Editor will automatically generate the element and
connector entries in the text text file. Defining inputs and outputs requires
editing the text file directly and assigning the inputs to the appropriate
element variables. For example, inflow(Q = Q_in)
sets the Q
variable
of the inflow
element equal to Q_in
.
In addition to elements, the input variables Q_in
and P_control
are
also defined. Q_in
is determined by the forecast and the operator cannot
control it, so we set Q_in(fixed = true)
. The actual values of Q_in
are stored in timeseries_import.csv
. P_control
is not defined anywhere
in the model or inputs we will dynamically assign its value every timestep in
the python script, \src\example.py
.
Because we want to view the water volume in the storage element in the output
file, we also define an output variable storage_V
and set it equal to the
corresponding state variable storage.V
. Dito for Q_release = P_control
.
The Simulation Problem¶
The python script is created and edited in a text editor. In general, the python script consists of the following blocks:
 Import of packages
 Definition of the simulation problem class
 Any additional configuration (e.g. overriding methods)
 A run statement
Importing Packages¶
Packages are imported using from ... import ...
at the top of the file. In
our script, we import the classes we want the class to inherit, the
package run_simulation_problem
form the rtctools.util
package, and
any extra packages we want to use. For this example, the import block looks
like:
1 2 3 4 5 6 7 8  from rtctools.simulation.simulation_problem import SimulationProblem
from rtctools.simulation.csv_mixin import CSVMixin
from rtctools.util import run_simulation_problem
import logging
logger = logging.getLogger("rtctools")
class Example(CSVMixin, SimulationProblem):

Simulation Problem¶
The next step is to define the simulation problem class. We construct the class by declaring the class and inheriting the desired parent classes. The parent classes each perform different tasks related to importing and exporting data and running the simulation problem. Each imported class makes a set of methods available to the our simulation class.
8  class Example(CSVMixin, SimulationProblem):

The next, we override any methods where we want to specify nondefault
behaviour. In our simulation problem, we want to define a proportional
controller. In its simplest form, we load the current values of the volume and
target volume variables, calculate their difference, and set P_control
to be
as close as possible to eliminating that difference during the upcoming timestep.
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41  def update(self, dt):
# Get the time step
if dt < 0:
dt = self._dt
# Get relevant model variables
volume = self.get_var('storage.V')
target = self.get_var('storage_V_target')
# Calucate error in storage.V
error = target  volume
# Calculate the desired control
control = error / dt
# Get the closest feasible setting.
bounded_control = min(max(control, self.min_release), self.max_release)
# Set the control variable as the control for the next step of the simulation
self.set_var('P_control', bounded_control)
# Call the super class so that everything else continues as normal
super(Example, self).update(dt)

Run the Simulation Problem¶
To make our script run, at the bottom of our file we just have to call
the run_simulation_problem()
method we imported on the simulation
problem class we just created.
44  run_simulation_problem(Example, log_level=logging.DEBUG)

The Whole Script¶
All together, the whole example script is as follows:
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 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44  from rtctools.simulation.simulation_problem import SimulationProblem
from rtctools.simulation.csv_mixin import CSVMixin
from rtctools.util import run_simulation_problem
import logging
logger = logging.getLogger("rtctools")
class Example(CSVMixin, SimulationProblem):
"""
A basic example for introducing users to RTCTools 2 Simulation
"""
# Min and Max flow rate that the storage is capable of releasing
min_release, max_release = 0., 8. # m^3/s
# Here is an example of overriding the update() method to show how control
# can be build into the python script
def update(self, dt):
# Get the time step
if dt < 0:
dt = self._dt
# Get relevant model variables
volume = self.get_var('storage.V')
target = self.get_var('storage_V_target')
# Calucate error in storage.V
error = target  volume
# Calculate the desired control
control = error / dt
# Get the closest feasible setting.
bounded_control = min(max(control, self.min_release), self.max_release)
# Set the control variable as the control for the next step of the simulation
self.set_var('P_control', bounded_control)
# Call the super class so that everything else continues as normal
super(Example, self).update(dt)
# Run
run_simulation_problem(Example, log_level=logging.DEBUG)

Running RTCTools¶
To run this basic example in RTCTools, navigate to the basic example src
directory in the RTCTools shell and run the example using python
example.py
. For more details about using RTCTools, see
Running RTCTools.
Extracting Results¶
The results from the run are found in output\timeseries_export.csv
. Any
CSVreading software can import it. Here we used matplotlib to generate this plot.
(Source code, png, hires.png, pdf)
Observations¶
This plot shows that the operator is not able to keep the water level within the bounds over the entire time horizon. They may need to increase the controller timestep, use a more complete PID controller, or use model predictive control such as the RTCTools optimization library to get the results they want.