...
 
Commits (67)
  • Tjerk Vreeken's avatar
    Fix history extrapolation for delayed expressions · 1f64e529
    Tjerk Vreeken authored
    The supposed block interpolation that was in place does not work as
    suggested. It would only properly extrapolate values if the corresponding
    output x-coordinate was outside of the range of input x-coordinates. It
    would however not extrapolate or interpolate like that with NaNs present.
    For example (y-coordinates of input and output only):
              [4,0, 4.0] --> [4.0, 4.0, 4.0, 4.0]
    [nan, nan, 4.0, 4.0] --> [nan, nan, 4.0, 4.0]
    1f64e529
  • Tjerk Vreeken's avatar
    Fix unwanted conversion to MX in substitute calls · 80c8656b
    Tjerk Vreeken authored
    It could happen that a DM expression would be converted to an MX
    expression, which would then later confuse the solver (e.g. when the
    values are the lower/upper bounds of the constraints).
    
    One example is that when passing as expr a 15x1 DM array of all 293.55
    values, the resulting value would be "MX(all_293.65(15x1))". This in
    turn would confuse the solver call, resulting in an exception.
    80c8656b
  • Tjerk Vreeken's avatar
    Remove linux tags from CI file · 0bf28c2a
    Tjerk Vreeken authored
    This allows us to run using GitLab's runners. If we ever end up running
    Windows docker images, we can then make Windows the exception with a
    "windows" tag.
    0bf28c2a
  • Tjerk Vreeken's avatar
    Fix goal consistency and monotonicity checks · 4b82a3ec
    Tjerk Vreeken authored
    If a goal had both a target min and a target max, the latter would not
    get checked for consistency due to the "elif" instead of "if".
    
    The check for target min/max exceedence of the function range compared
    the True/False value returned by has_target_* instead of the actual
    target min/max value.
    
    Furthermore, the monotonicity checks path goals were behind the option
    flag, but those of (normal) goals were not.
    4b82a3ec
  • Tjerk Vreeken's avatar
    Fix goal constraint min/max bounds consistency · 40567dba
    Tjerk Vreeken authored
    It could happen that lower (later) priority goals for an already
    constrained state would result in an answer that would be slightly
    outside the hard constraints due to feasibilty tolerances.
    When combining the soft constraint and existing hard constraint in a new
    hard constraint, the existing hard constraint values could be shifted,
    basically enforcing a solution that may or may not be feasible.
    
    An example for a typical state with goals at multiple priorities:
    
    +----------+------------+------------+
    | Priority | Target Min | Target Max |
    +----------+------------+------------+
    |        1 |        0.0 |        0.7 |
    |        2 |        0.7 |        0.7 |
    +----------+------------+------------+
    
    The first priority could end up with a hard constraint exactly matching,
    i.e. (0.0, 0.7). Exaggerating the numerical differences that can occur,
    the second priority could end up with a value of 0.701 for the state in
    question.
    What would happen before this commit is that we would calculate a new
    hard constraint (0.701, 0.7). This would then be corrected to (0.701,
    0.701), which would then replace the original hard constraint.
    
    It is safer to keep enforcing the old hard constraint (which led to the
    solution of 0.701) making the new constraint (0.7, 0.7). This is the
    behavior that this commit makes happen.
    40567dba
  • Tjerk Vreeken's avatar
    Fix reading csv files in documentation of examples · 90bfb70b
    Tjerk Vreeken authored
    By default, NumPy does not convert bytes to strings. Specifying
    "encoding=None" will use the system default to produce strings.
    
    See also 962659ce, which fixed a similar issue for CSVMixin.
    90bfb70b
  • Tjerk Vreeken's avatar
    Documentation fix for reference to LookupTable · 6bb630e3
    Tjerk Vreeken authored
    Regression in 91af5370.
    6bb630e3
  • Tjerk Vreeken's avatar
    Documentation fix for lookup tables · 6ed9fd84
    Tjerk Vreeken authored
    6ed9fd84
  • Tjerk Vreeken's avatar
    Documentation fix for SimulationProblem · 77c90eb0
    Tjerk Vreeken authored
    Both "finalize" and "get_options" were removed in the conversion from
    JModelica/FMU-based simulation to Pymoca/CasADi-based simulation.
    77c90eb0
  • Tjerk Vreeken's avatar
    Documentation fix for LinearizationMixin · 060e0247
    Tjerk Vreeken authored
    The documentation page was not referenced, and therefore not listed or
    reachable.
    060e0247
  • Tjerk Vreeken's avatar
    b0ba8172
  • Tjerk Vreeken's avatar
    1dfcbfb6
  • Tjerk Vreeken's avatar
    Documentation restore tight x-axis layout · d9d28d8e
    Tjerk Vreeken authored
    Making plots with Matplotlib 2+ has whitespace on both ends of the x-axis,
    which is not what we want. To restore the previous behavior, we explicitly
    tightly scale the x-axis.
    d9d28d8e
  • Jorn Baayen's avatar
    Add 'initial_state' method to documentation. · 4681c30c
    Jorn Baayen authored
    4681c30c
  • Tjerk Vreeken's avatar
    Fix big-M in goal programming example · 281ae441
    Tjerk Vreeken authored
    Not sure why a very large value was chosen, as it is well established
    that big-M values should typically be as small as possible. Very high
    values can result in bad relaxations and make branching very hard.
    
    In this case, a more reasonable M of 2 is both faster (~2x), and also
    results in a better solution. The first two priorities remain the same,
    but the third priority goes from 12.23 to 10.85.
    281ae441
  • Tjerk Vreeken's avatar
    Doc update of getting started instructions · aeea1515
    Tjerk Vreeken authored
    Many basic installation and running instructions were still aimed at
    RTC-Tools 2.0 and the way that was packaged. The most important changes
    to 2.2 are:
    
    - Using pip instead of pre-packaged installer
    - Having to explicitly download examples (if installing from PyPI)
    - Copying Modelica libraries to a more convenient location to load them
      in OMEdit.
    aeea1515
  • Tjerk Vreeken's avatar
    efc13d64
  • Tjerk Vreeken's avatar
    Doc fixes for relative references · 313ee467
    Tjerk Vreeken authored
    313ee467
  • Tjerk Vreeken's avatar
    Doc typo and formatting fixes · 1a4be698
    Tjerk Vreeken authored
    1a4be698
  • Tjerk Vreeken's avatar
    Use delay() operator in tests · c3043141
    Tjerk Vreeken authored
    Some tests were already addressed in 14c1702e, but not all.
    c3043141
  • Tjerk Vreeken's avatar
    Fix parameter inlining for delayed expressions · fa49e368
    Tjerk Vreeken authored
    If we inline a parameter, we have to make sure to also inline it in the
    delay expression and duration. The latter was already taken care of
    somewhere else in the code, but has been moved near the newly added
    processing of the delay expressions.
    
    Note that evaluating the delayed_feedback_function() could still
    give non-NaN results even if inputs were missing, which made this bug
    somewhat hard to encounter and track down.
    The fact that this happens is due to the expand() call; if we would call
    the original MX function, it would complain about free variables. Why
    CasADi does not always return NaN for outputs of which not all inputs
    were specified is unknown though.
    fa49e368
  • Tjerk Vreeken's avatar
    Fix parameter inlining for path objective/constraints · 0b6c34f7
    Tjerk Vreeken authored
    When using goal programming, parameters would only get inlined in the
    first priority (i.e. the first call to the transcribe() method).
    
    Note that regular objective and constraints cannot/should not refer to
    parameters with the variable() or state() methods (which return an MX
    symbol), and instead always rely on state_at() and parameters(). The
    latter two always return a numberic value, and we therefore we do not
    need to perform any inlining on the regular objective and constraints.
    0b6c34f7
  • Tjerk Vreeken's avatar
    Fix theta parameter availability · 4df6d5bb
    Tjerk Vreeken authored
    When using the HomotopyMixin, a call to parameters() would fail in
    pre(), because the parameter __theta is only initialized in the
    optimize() call.
    
    Also add a comment explaining why we do not initialize theta right away.
    4df6d5bb
  • Teresa Piovesan's avatar
    399a7f7c
  • Tjerk Vreeken's avatar
    Fix discrepancy between goal and path goal scaling · cd33c423
    Tjerk Vreeken authored
    Path goals were not divided by the number of time steps, meaning that
    formulting a problem with a bunch of normal goals using state_at() (one
    for each time step) and a path goal would lead to different objectives and
    results.
    
    To fix the discrepancy, we disable the division by the number of
    goals/time steps by default. This is conservative in the sense that we
    then guarantee at least as accurate a result as before this commit. For
    very tricky/badly scaled problems it may however lead to convergence
    issues.
    
    Should scaling with the problem size be desired, the new
    "scale_by_problem_size" goal programming option can be set to True.
    
    Closes #1046
    cd33c423
  • Tjerk Vreeken's avatar
    bde691b4
  • Teresa Piovesan's avatar
    Doc fix type of scale_by_problem_size · 9ce64c09
    Teresa Piovesan authored
    9ce64c09
  • Olav van Duin's avatar
    Fix flooring of dates/times in PI-XML import · a2f00aff
    Olav van Duin authored
    Before, we would floor based on the base time (midnight 1-1-1970). This
    would mean that having a start date of e.g. 09:35 and a time step of
    00:15 did not work well together. The start date would get rounded down
    to 09:30 instead of the desired 09:35. This commit makes sure that such
    cases are handled correctly, but that rounding is still performed (e.g.
    09:51 will get rounded down to 09:50).
    
    Closes #1020
    a2f00aff
  • Tjerk Vreeken's avatar
    Disallow function_range for minimization goals · 5f3e17df
    Tjerk Vreeken authored
    Minimization goals only use the function nominal for scaling. Requiring
    the user to specify a function range is misleading.
    
    Closes #1051
    5f3e17df
  • Tjerk Vreeken's avatar
    Improve goal function range vs. targets check · f97a0ab4
    Tjerk Vreeken authored
    We now check that the target min/max are not equal to the function range
    lower bound/upper bound respectively as well.
    
    A target minimum equal to the lower bound of the function range does not
    make any sense, as the corresponding epsilon is then free (and will
    therefore always become zero). Effectively, the goal will behave as a
    critical goal, and users should use those instead if such behavior is
    desired.
    f97a0ab4
  • Olav van Duin's avatar
    Only do duplicate parameter check for PI-XML · 094cdac4
    Olav van Duin authored
    Before, we would (optionally) check if a parameter read from the PI-XML
    files (rtcParameterConfig*.xml) was already set before. This would lead
    to a warning if a parameter occurred in multiple PI-XML files, but also
    if a value was already present in the model file(s). The latter type of
    warning is unnecessary as we expect parameters to be updated at runtime.
    This commit makes sure that only warnings are raised for duplicates
    found in the PI-XML files.
    094cdac4
  • Teresa Piovesan's avatar
    Fix examples to comply to new funct_range check · bb84fbfd
    Teresa Piovesan authored
    User can not specify (anymore) the function range for minimization goals.
    Change the examples to comply with the check introduced in commit
    e9f500df.
    bb84fbfd
  • Jorn Baayen's avatar
    CollInt: Handle MX objective with 0x0 shape · 3bb2217b
    Jorn Baayen authored
    If self.objective() returns a 0x0 MX, do not use it as a base to add
    path objectives to, as the result would still be 0x0.
    
    This is due to the fact that `ca.MX() + 1 = ca.MX()`.
    3bb2217b
  • Jorn Baayen's avatar
    Return ca.MX(0) by default in objectives/residuals · db12302f
    Jorn Baayen authored
    Because `ca.MX() + 1 = ca.MX()`, returning an empty symbol could be
    confusing if people override e.g. self.objective().
    db12302f
  • Teresa Piovesan's avatar
    Doc fix goal function range/nominals · f413630d
    Teresa Piovesan authored
    Add explanation on how the function range and function nominal are used
    for minimization goals vs target goals. This also fixes consistency with
    the check introduced in commit e9f500d7.
    f413630d
  • Tjerk Vreeken's avatar
    Simulation: opt-in partial workaround for delay() · 1ca2c04d
    Tjerk Vreeken authored
    A workaround is added to support a delay of zero (or make delay expression
    behave as such). By default, this workaround is disabled. The user has to
    explicitly set "_force_zero_delay" to True to enable it.
    1ca2c04d
  • Teresa Piovesan's avatar
    Doc fix goal_programming example · 72053fbb
    Teresa Piovesan authored
    Delete in doc example reference to function range for
    minimization goal.
    72053fbb
  • Teresa Piovesan's avatar
    Fix error in scale_by_problem_size implementation · eb3a953d
    Teresa Piovesan authored
    Typo in the scaling of the accumulated objective in the path_objective
    definition.
    eb3a953d
  • Ivo's avatar
    Fix function_range in initial state estimator goals · bef18ad3
    Ivo authored
    Since e9f500d7 it is no longer allowed to specify a function_range for a
    minimization goal. Also, the scaling is only done with the
    function_nominal since 5ac6be99, and function_nominal should therefore
    be `max_deviation`.
    
    Closes #1062
    bef18ad3
  • Jesse VanderWees's avatar
    Doc note in MI example about solver support · 8045046a
    Jesse VanderWees authored
    Solver support is a FAQ not adequately addressed in the documentation.
    We add a note to the mixed-integer example directing people to
    solver_options().
    8045046a
  • Tjerk Vreeken's avatar
    ad7366d0
  • Teresa Piovesan's avatar
    Fix goal relaxation with minimization goals · f7d3b695
    Teresa Piovesan authored
    The default value for fix_minimized_values is True, and therefore a goal
    relaxation on minimization goals would result in two constraints, one of
    which would be possibly non-convex.
    
    The idea behind fix_minimized_values is to impose a possibly non-convex
    _equality_ constraint for interior point methods, to avoid numerical
    issues with small slack variables. We therefore now check if the goal
    has any relaxation, in which case we impose a single inequality
    constraint instead of a single equality constraint.
    f7d3b695
  • Teresa Piovesan's avatar
    Handle fix_minimized_values for normal goals · 8ece605d
    Teresa Piovesan authored
    The fix_minimized_values goal programming option only affected path
    goals. The same logic should however apply to normal goals.
    8ece605d
  • Teresa Piovesan's avatar
    Make fix_minimized_values False for non-default solvers · e9e4ecfa
    Teresa Piovesan authored
    Solvers like Gurobi and CPLEX cannot handle quadratic equality
    constraints. The default value for fix_minimized_values being True might
    therefore lead to issues if more than one minimization goal is present.
    
    To provide backwards compatibility, we make sure that the default value
    of True is preserverd for IPOPT and BONMIN. Only in the case that the
    user overrides the solver will the default value be False.
    e9e4ecfa
  • Tjerk Vreeken's avatar
    Fix big-M in mixed integer example · c19ab862
    Tjerk Vreeken authored
    Similar to the big-M fix for the goal proramming example done in
    commit a1101002.
    c19ab862
  • Tjerk Vreeken's avatar
    Fix scale_by_problem_size for path target goals · c103e18f
    Tjerk Vreeken authored
    If a path goal has a target min and/or max, the epsilon will be handled in
    the normal objective (not path_objective()). We therefore have to make
    sure to also divide by the number of time steps after summing an epsilon
    over time if the user has enabled the "scale_by_problem_size" options.
    c103e18f
  • Tjerk Vreeken's avatar
    Improve goal scale by problem size · e7b7c863
    Tjerk Vreeken authored
    The weighing of goals and path goals should not depend on whether they are
    part of the objective or path objective. To make sure it does not, we have
    to divide both objectives by the _total_ number of subproblem objectives.
    e7b7c863
  • Tjerk Vreeken's avatar
    Inlude function_range in repr() of StateGoal · efbe3e71
    Tjerk Vreeken authored
    The function_range was already added to the repr() of Goal, but not to
    StateGoal yet. We add it there as well to provide helpful information
    when e.g. the function_range derived from the bounds and the target
    min/max do not agree.
    efbe3e71
  • Tjerk Vreeken's avatar
    Rename constraint_relaxation to violation_relaxation · d3e1d3a8
    Tjerk Vreeken authored
    The "constraint_relaxation" name was confusing, as the relaxation is
    only indirectly applied to the eventual hard constraint (i.e. it is
    scaled by the goal's function nominal as well). The term
    violation_relaxation much better describes that this relaxation is
    applied to the violation variables (a.k.a. epsilons).
    d3e1d3a8
  • Tjerk Vreeken's avatar
    Change solver_output from 2-D to 1-D array · b4141be4
    Tjerk Vreeken authored
    Most arrays like lbg/ubg/x0 are already 1-D NumPy arrays if they are not
    a CasADi matrix (MX/DM). It therefore makes sense if the output of the
    solver is also a 1-D array, requiring only a single index to be
    specified to get an element (instead of X[i, 0]).
    
    We ensure that the results dictionary does not present views onto the
    raw solver_output by making a full copy of it first. This is just a
    safeguard though, as it is undefined behavior what should happen if
    someone changes the arrays in the results like that.
    
    Contrary to the previous type annotation, the solver output was already
    a numpy array (a 2-D one, see commit 40875b0a).
    
    Closes #1064
    b4141be4
  • Tjerk Vreeken's avatar
    Pin CasADi and Pymoca versions · 89ffc704
    Tjerk Vreeken authored
    These two packages are relatively unstable or known to make
    substantial (backwards-incompatible) changes. We therefore pin them to
    release series known/expected to work, while allowing bug fixes using a
    wildcard specifier.
    89ffc704
  • Tjerk Vreeken's avatar
    Fix formatting of default parameters with types · a93ccfdb
    Tjerk Vreeken authored
    According to Python's formatting rules, default parameters without types
    should have no spaces around the equals operator. However, if the type of
    a parameter _is_ specified, spaces around the operator are required. See
    also PEP-3107 for examples of such code.
    
    Flake8 only started reporting this formatting error starting with version
    3.6 (error code E252), causing our test pipeline to fail even though there
    were no changes to the codebase. This commit makes sure that the
    formatting checks pass again.
    a93ccfdb
  • Tjerk Vreeken's avatar
    Fix formatting of broken escape sequences · 0a119b4c
    Tjerk Vreeken authored
    Flake8 only started reporting this formatting error starting with
    version 3.6 (error code W605), causing our test pipeline to fail even
    though there were no changes to the codebase. This commit makes sure
    that the formatting checks pass again.
    0a119b4c
  • Teresa Piovesan's avatar
    Fix seed derivatives · 25856836
    Teresa Piovesan authored
    The value of the derivatives of the initial time step and of the extra
    variables is included the seed.
    
    Closes issue #1066
    25856836
  • Jesse VanderWees's avatar
    Fix saving function values of critical goals · a1a3e66f
    Jesse VanderWees authored
    Function values of critical goals were not being written to timeseries,
    even when the variable function_value_timeseries_id was set.
    a1a3e66f
  • Jorn Baayen's avatar
    CSVMixin: Output CSV files with time stamps from times() method · 968c5456
    Jorn Baayen authored
    The data interpolation assumes that self.times() dictates the time stamps
    for the output data. The time stamps in the CSV output file, however, were
    taken from the input CSV file and not from the times() method. This
    implicitly assumed that self.times() equates the time steps from the CSV
    input file. This is not always the case, as self.times() may be overridden
    by a subclass.
    968c5456
  • Jorn Baayen's avatar
    PIMixin: Output PI files with time stamps from times() method · 9ceda84d
    Jorn Baayen authored
    Similar fix to previous commit.
    9ceda84d
  • Jesse VanderWees's avatar
    Properly handle LookupTable __call__() parameters · 0b4c7cf9
    Jesse VanderWees authored
    The handling of LookupTable __call__() parameters was somewhat sloppy,
    without clear error codes if parameters were not supported. Furthermore,
    calling with np.nan would expose undefined behaviour. This commit
    improves error messages and does explicit handling of np.nan parameters.
    0b4c7cf9
  • Tjerk Vreeken's avatar
    Fix reading symbolic nominals in ModelicaMixin · 5b1d8e4e
    Tjerk Vreeken authored
    Earlier on in the affected method we cast to make sure that we always deal
    with a ca.MX instance. However, the substitute call is likely to return a
    ca.DM, which would then fail later on when calling ".to_DM()" on it.
    
    Closes #1081
    5b1d8e4e
  • Tjerk Vreeken's avatar
    Fix failing constraint() call before optimize() · b9a30ca4
    Tjerk Vreeken authored
    When using GoalProgrammingMixin, it was not possible to call constraint()
    or path_constraints() in e.g. pre(). This is because the subproblem
    constraint variables were not fully initialized yet; they were not
    indexable with "ensemble_member". Only when optimize() was called were
    the variables properly initialized.
    b9a30ca4
  • Jesse VanderWees's avatar
    Don't set nan constant inputs in modelica mixin · ad65d4af
    Jesse VanderWees authored
    Modelica mixin was masking missing timeseries detection by setting
    timeseries full of NaN. This was causing hard-to-debug
    InvalidNumberDetected solver errors.
    ad65d4af
  • Tjerk Vreeken's avatar
    Fix interpolation call in map_path_expression · 2d4a1c73
    Tjerk Vreeken authored
    This could lead to an exception being thrown by CasADi, because the
    integer for "mode" would be passed as the boolean for "equidistant". We
    are lucky that there was no implicit casting of int->bool, or this would
    have turned into an obscure bug.
    2d4a1c73
  • Jesse VanderWees's avatar
    Fix log_solver_failure_as_error being ignored · 5d8228d9
    Jesse VanderWees authored
    Commit b2dc65d6 caused the log_solver_failure_as_error parameter to
    optimize() to be ignored.
    5d8228d9
  • Tjerk Vreeken's avatar
    SimulationProblem: Fix scaling of equations with nominals · b4ccb638
    Tjerk Vreeken authored
    SimulationProblem was only ever using the nominals defined in the Modelica
    file, even if e.g. the user was overriding get_variable_nominal().
    
    Also note that variable nominals can be symbolic, but that we would like
    to handle numeric values/floats as much as possible. We therefore
    evaluate the nominals to their numeric values as soon as we know all
    parameter values.
    b4ccb638
  • Tjerk Vreeken's avatar
    SimulationProblem: Fix scaling of bounds with nominals · 2d4c12fb
    Tjerk Vreeken authored
    The bounds were not scaled at all with the nominals. With the equations
    being scaled, this would lead to many infeasible problems.
    2d4c12fb
  • Teresa Piovesan's avatar
    Fix timeseries import PIMixin for ensembles · 4c3af937
    Teresa Piovesan authored
    The looping in PIMixin was done correctly, but the underlying Timeseries
    class in pi.py only ever returned the values of the first ensemble
    member.
    4c3af937
  • Teresa Piovesan's avatar
    Fix timeseries export PIMixin for ensembles · 196d1fb7
    Teresa Piovesan authored
    The underlying Timeseries class in pi.py raised an AssertionError when
    trying to set the "ensemble_size" without first setting
    "contains_ensemble" to True.
    196d1fb7
......@@ -46,7 +46,6 @@ py35:linux:
TOXENV: py35
tags:
- docker
- linux
py36:linux:
<<: *unittests
......@@ -55,7 +54,6 @@ py36:linux:
TOXENV: py36
tags:
- docker
- linux
# coverage
coverage:
......
......@@ -165,7 +165,7 @@ html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
# html_static_path = ['_static']
# Add any extra paths that contain custom files (such as robots.txt or
# .htaccess) here, relative to this directory. These files are copied
......
......@@ -20,9 +20,14 @@ does not want to end up with too much water at the end of the six days. They
have chosen to use RTC-Tools to calculate how much water to release and when
to release it.
The folder ``<installation directory>\RTCTools2\examples\basic``
contains a complete RTC-Tools optimization problem. An RTC-Tools
directory has the following structure:
If you installed using source, the library and examples directory are
available in the git repositories. If you installed using pip directly, you
first need to download/copy the examples and libraries to a convenient
location. See :ref:`getting-started-download-examples` and
:ref:`getting-started-copy-libraries` for detailed instructions.
The folder ``<examples directory>\basic`` contains a complete RTC-Tools
optimization problem. An RTC-Tools directory has the following structure:
* ``input``: This folder contains the model input data. These are several files
in comma separated value format, ``csv``.
......@@ -40,17 +45,15 @@ 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 :ref:`getting-started-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``
* Select ``<library directory>\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``
* Select ``<examples directory>\basic\model\Example.mo``
Once loaded, we have an OpenModelica Connection Editor window that looks like
this:
......
......@@ -4,12 +4,12 @@ Using an Ensemble Forecast
.. note::
This example is an extension of :doc:`lookup_table`. It
assumes prior knowledge of goal programming and the lookuptables components
of RTC-Tools. If you are a first-time user of RTC-Tools, see
assumes prior knowledge of goal programming and the lookup tables components
of RTC-Tools. If you are a first-time user of RTC-Tools, see
:doc:`basic`.
Then biggest change to RTC-Tools when using an ensemble is the structure of the
directory. The folder ``RTCTools2\examples\ensemble``
directory. The folder ``<examples directory>\ensemble``
contains a complete RTC-Tools ensemble optimization problem. An RTC-Tools
ensemble directory has the following structure:
......@@ -23,12 +23,12 @@ ensemble directory has the following structure:
* ``ensemble.csv``: a file where the names and probabilities of the ensemble
members are defined
* ``forcast1``
* ``forecast1``
* the file ``timeseries_import.csv``
* the file ``initial_state.csv``
* ``forcast2``
* ``forecast2``
* ``timeseries_import.csv``
* ``initial_state.csv``
......@@ -37,11 +37,11 @@ ensemble directory has the following structure:
* ``output``: The folder where the output is saved:
* ``forcast1``
* ``forecast1``
* ``timeseries_export.csv``
* ``forcast2``
* ``forecast2``
* ``timeseries_export.csv``
......@@ -246,8 +246,8 @@ 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 CSV-reading software can
The results from the run are found in ``output/forecast1/timeseries_export.csv``
and ``output/forecast2/timeseries_export.csv``. Any CSV-reading software can
import it, but this is how results can be plotted using the python library
matplotlib:
......
......@@ -91,7 +91,7 @@ Declaring Goals
'''''''''''''''
Goals are defined as classes that inherit the ``Goal`` parent class. The
components of goals can be found in :doc:`../optimization/multi_objective`. In
components of goals can be found in :doc:`../../optimization/multi_objective`. In
this example, we demonstrate three ways to define a goal in RTC-Tools.
First, we have a high priority goal to keep the water level within a minimum and
......@@ -105,9 +105,8 @@ this goal, called a ``StateGoal``:
:lineno-match:
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
non-path 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:
``Q_pump``. This goal has a lower priority than the water level range goal.
This goal does not use a helper class:
.. literalinclude:: ../../../examples/goal_programming/src/example.py
:language: python
......
......@@ -65,7 +65,7 @@ Declaring Goals
'''''''''''''''
Goals are defined as classes that inherit the ``Goal`` parent class. The
components of goals can be found in :doc:`../optimization/multi_objective`. In
components of goals can be found in :doc:`../../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
......@@ -126,7 +126,7 @@ 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 :doc:`../optimization/basics`.
file. You can find more information on timeseries here :doc:`../../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
......
......@@ -12,6 +12,14 @@ Mixed Integer Optimization: Pumps and Orifices
hydraulic model, and assumes basic exposure to RTC-Tools. To start with
basics, see :doc:`basic`.
.. note::
By default, if you define any integer or boolean variables in the model, RTC-Tools
will switch from IPOPT to BONMIN. You can modify solver options by overriding
the
:meth:`solver_options()<rtctools.optimization.optimization_problem.OptimizationProblem.solver_options>`
method. Refer to CasADi's nlpsol interface for a list of supported solvers.
The Model
---------
......@@ -34,7 +42,7 @@ 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``
``<examples directory>\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
......
......@@ -11,7 +11,7 @@ delimiter = ','
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:]
dtype='object' + ',float' * (ncols - 1), names=True, encoding=None)[1:]
# Generate Plot
n_subplots = 2
......@@ -42,5 +42,7 @@ for i in range(n_subplots):
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)
plt.autoscale(enable=True, axis='x', tight=True)
# Output Plot
plt.show()
......@@ -14,7 +14,7 @@ def get_results(forecast_name):
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)
dtype='object' + ',float' * (ncols - 1), names=True, encoding=None)
# Generate Plot
n_subplots = 2
......@@ -54,5 +54,7 @@ for i in range(len(axarr)):
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)
plt.autoscale(enable=True, axis='x', tight=True)
# Output Plot
plt.show()
......@@ -11,7 +11,7 @@ delimiter = ','
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:]
dtype='object' + ',float' * (ncols - 1), names=True, encoding=None)[1:]
# Generate Plot
n_subplots = 3
......@@ -49,5 +49,7 @@ for i in range(n_subplots):
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)
plt.autoscale(enable=True, axis='x', tight=True)
# Output Plot
plt.show()
......@@ -11,7 +11,7 @@ delimiter = ','
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:]
dtype='object' + ',float' * (ncols - 1), names=True, encoding=None)[1:]
# Generate Plot
n_subplots = 2
......@@ -42,5 +42,7 @@ for i in range(n_subplots):
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)
plt.autoscale(enable=True, axis='x', tight=True)
# Output Plot
plt.show()
......@@ -11,7 +11,7 @@ delimiter = ','
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:]
dtype='object' + ',float' * (ncols - 1), names=True, encoding=None)[1:]
# Generate Plot
f, axarr = plt.subplots(2, sharex=True)
......@@ -41,5 +41,7 @@ for i in range(len(axarr)):
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)
plt.autoscale(enable=True, axis='x', tight=True)
# Output Plot
plt.show()
......@@ -11,7 +11,7 @@ delimiter = ','
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)
dtype='object' + ',float' * (ncols - 1), names=True, encoding=None)
# Generate Plot
n_subplots = 2
......@@ -42,5 +42,7 @@ for i in range(n_subplots):
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)
plt.autoscale(enable=True, axis='x', tight=True)
# Output Plot
plt.show()
......@@ -23,7 +23,7 @@ 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``
The folder ``<examples directory>\simulation``
contains a complete RTC-Tools simulation problem. An RTC-Tools
directory has the following structure:
......@@ -48,12 +48,12 @@ 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``
* Select ``<library directory>\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``
* Select ``<examples directory\simulation\model\Example.mo``
Once loaded, we have an OpenModelica Connection Editor window that looks like
this:
......
......@@ -9,80 +9,86 @@ For most users, the easiest way to install RTC-Tools using the `pip <https://pip
Using the Pip Package Manager
-----------------------------
Although not required, it is recommended to install RTC-Tools in a virtual
environment. See the `official Python tutorial
<https://docs.python.org/3/tutorial/venv.html>`_ for more information on how
to set up and activate a virtual environment.
RTC-Tools, including its dependencies, can be installed using the `pip <https://pip.pypa.io/>`_ package manager::
# Install RTC-Tools using pip package manager
python3 -m pip install rtctools
# Install RTC-Tools and Channel Flow using pip package manager
pip install rtc-tools rtc-tools-channel-flow
From Source
-----------
Although not required, it is recommended to build and install RTC-Tools in a `virtual environment
<https://virtualenv.pypa.io/en/stable/>`_.
Dependencies
~~~~~~~~~~~~
RTC-Tools 2 requires `Python <https://www.python.org>`_ >= 3.6 (*not Python 2*) with the following packages:
The latest RTC-Tools and Channel Flow source can be downloaded using git::
- `numpy <https://www.numpy.org/>`_ >= 1.11.10
# Get RTC-Tools source
git clone https://gitlab.com/deltares/rtc-tools.git
- `pymoca <https://github.com/pymoca/pymoca/>`_
# Get RTC-Tools's Modelica library
git clone https://gitlab.com/deltares/rtc-tools-channel-flow.git
- `CasADi <https://github.com/casadi/casadi/>`_ >= 3.2.0
Then you can install this latest version as follows::
The dependencies can be installed using the `pip <https://pip.pypa.io/>`_ package manager::
pip install ./rtc-tools
pip install ./rtc-tools-channel-flow
# Install dependencies using pip
python3 -m pip install numpy pymoca casadi
Or if you would like to have an editable installation (e.g. as developer)::
Acquiring the source
~~~~~~~~~~~~~~~~~~~~
pip install -e ./rtc-tools
pip install -e ./rtc-tools-channel-flow
The latest RTC-Tools source can be downloaded using git::
# Get RTC-Tools source
git clone https://gitlab.com/deltares/rtc-tools.git
.. _getting-started-download-examples:
# Get RTC-Tools's Modelica library
git clone https://gitlab.com/deltares/rtc-tools-channel-flow.git
Downloading and running examples
================================
Linux/MacOS-X
~~~~~~~~~~~~~
To check whether the installation was succesful, the basic example can be
used. If RTC-Tools was not installed from source, the examples need to be
downloaded first::
To actually build and install RTC-Tools, run::
# Download the examples to the current folder (.)
rtc-tools-download-examples .
# Build and install RTC-Tools
python setup.py install
# Navigate to the basic example
cd rtc-tools-examples/basic/src
To check whether the installation was succesful, the basic example can be
used. Only the environment variable pointing to the
Deltares Modelica library remains for the user to set::
# Run the example
python example.py
export DELTARES_LIBRARY_PATH=\`readlink -f ../rtc-tools-channel-flow\`
If the installation was succesful, you should see that the solver succeeds:
cd examples/basic/src
.. image:: images/basic_example_console.png
# Set the correct environment variables, and run the example
python3 example.py
Elsewhere in this documentation we refer to the folder containing the examples
as ``<examples directory>``. Depending on the method of installation this can
then either be:
Windows
~~~~~~~
* ``\path\to\rtc-tools-examples``, when having downloaded the examples
* ``\path\to\source\of\rtc-tools\examples``, when having installed RTC-Tools from source
To check whether the installation was succesful, the basic example can be
used. Only the environment variable pointing to the
Deltares Modelica library remains for the user to set::
.. _getting-started-copy-libraries:
set DELTARES_LIBRARY_PATH=C:\path\to\rtc-tools-channel-flow
Copying Modelica libraries
==========================
cd /D C:\path\to\rtc-tools\basic\src
Because the Modelica libraries are distributed as pip packages, their location
inside Python's site-packages can be somewhat inconvient. To copy the Modelica
libraries to a more convenient location, you can use the ``rtc-tools-copy-libraries``
command::
# Set the correct environment variables, and run the example
python3 example.py
# Copy all Modelica libraries of RTC-Tools to the current folder (.)
rtc-tools-copy-libraries .
If the installation was succesful, you should see that the solver succeeds:
You should now have a folder ``Deltares``, containing amongst others a
``package.mo`` file, a ``ChannelFlow`` folder and folders of any other RTC-
Tools extensions you installed.
.. image:: images/basic_example_console.png
Elsewhere in this documentation we refer to the library folder containing the
``Deltares`` folder as ``<library directory>``.
.. _getting-started-omedit:
......@@ -101,7 +107,7 @@ 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, :doc:`examples/basic`.
example, :doc:`examples/optimization/basic`.
.. _running-rtc-tools:
......@@ -109,17 +115,15 @@ example, :doc:`examples/basic`.
Running RTC-Tools
=================
RTC-Tools is run from a command line shell. On Windows, both ``PowerShell``
and ``cmd`` can be used. On Linux/MacOS you could use the terminal application
with a shell of your liking.
RTC-Tools is run from a command line shell. If you installed using the Windows
executable, the RTC-Tools Shell can be started by clicking::
Start -> All Programs -> RTC-Tools -> Shell
Once you have started the shell, navigate to the ``src`` directory of the case
you wish to optimize, e.g.::
Once you have started the shell and loaded the correct virtual environment (if
applicable), navigate to the ``src`` directory of the case you wish to
optimize, e.g.::
cd \path\to\RTCTools2\examples\basic\src
cd \path\to\rtc-tools-examples\basic\src
Then, to run the case with RTC-Tools, run the ``src`` python script, e.g.::
......
......@@ -12,7 +12,7 @@
RTC-Tools documentation
=======================
This is the documentation for RTC-Tools 2.0, the `Deltares
This is the documentation for RTC-Tools, the `Deltares
<https://www.deltares.nl/>`_ toolbox for control and optimization of
environmental systems.
......
......@@ -11,6 +11,7 @@ Contents:
optimization/modelica_models
optimization/csv_io
optimization/fews_io
optimization/linearization
optimization/lookup_tables
optimization/homotopy
optimization/initial_state_estimation
......
......@@ -6,7 +6,7 @@ Basics
:show-inheritance:
.. autoclass:: rtctools.optimization.optimization_problem.OptimizationProblem
:members: bounds, constant_inputs, constraints, control, control_at, delayed_feedback, der, der_at, ensemble_member_probability, ensemble_size, get_timeseries, history, initial_time, integral, interpolate, lookup_tables, objective, optimize, parameters, path_constraints, path_objective, post, pre, seed, set_timeseries, solver_options, solver_success, state, state_at, states_in, timeseries_at
:members: bounds, constant_inputs, constraints, control, control_at, delayed_feedback, der, der_at, ensemble_member_probability, ensemble_size, get_timeseries, history, initial_state, initial_time, integral, interpolate, lookup_tables, objective, optimize, parameters, path_constraints, path_objective, post, pre, seed, set_timeseries, solver_options, solver_success, state, state_at, states_in, timeseries_at
:show-inheritance:
.. autofunction:: rtctools.util.run_optimization_problem
\ No newline at end of file
.. autofunction:: rtctools.util.run_optimization_problem
Lookup tables
=============
.. autoclass:: rtctools.optimization.optimization_problem.LookupTable
.. autoclass:: rtctools.optimization.csv_lookup_table_mixin.LookupTable
:members: __call__
:show-inheritance:
.. autoclass:: rtctools.optimization.csv_lookup_table_mixin.CSVLookupTableMixin
:members: -
:members: lookup_tables
:show-inheritance:
\ No newline at end of file
......@@ -2,7 +2,7 @@ Basics
======
.. autoclass:: rtctools.simulation.simulation_problem.SimulationProblem
:members: initialize, pre, post, setup_experiment, finalize, update, simulate, reset, get_start_time, get_end_time, get_current_time, get_options, get_var, get_var_count, get_var_name, get_var_type, get_var_rank, get_var_shape, get_variables, get_parameter_variables, get_input_variables, get_output_variables, compiler_options
:members: initialize, pre, post, setup_experiment, update, simulate, reset, get_start_time, get_end_time, get_current_time, get_var, get_var_count, get_var_name, get_var_type, get_var_rank, get_var_shape, get_variables, get_parameter_variables, get_input_variables, get_output_variables, compiler_options
:show-inheritance:
.. autofunction:: rtctools.util.run_simulation_problem
\ No newline at end of file
......@@ -13,11 +13,10 @@ from rtctools.util import run_optimization_problem
class WaterVolumeRangeGoal(StateGoal):
def __init__(self, optimization_problem):
# Call super class first, and pass in the optimization problem
super().__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')
super().__init__(optimization_problem)
state = 'storage.V'
priority = 1
......@@ -31,8 +30,8 @@ class MinimizeQreleaseGoal(StateGoal):
order = 1
class Example(GoalProgrammingMixin, ControlTreeMixin, CSVLookupTableMixin,
CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):
class Example(GoalProgrammingMixin, CSVMixin, CSVLookupTableMixin, ModelicaMixin,
ControlTreeMixin, CollocatedIntegratedOptimizationProblem):
"""
An extention of the goal programming and lookuptable examples that
demonstrates how to work with ensembles.
......
time,Q_orifice,Q_pump,is_downhill,sea_level,storage_level
2013-05-19 22:00:00,0.000000,0.000000,1.000000,0.000000,0.400000
2013-05-19 23:00:00,0.000000,0.000000,1.000000,0.100000,0.418000
2013-05-20 00:00:00,1.671714,0.000000,1.000000,0.200000,0.429982
2013-05-20 01:00:00,3.889201,0.000000,1.000000,0.300000,0.433981
2013-05-20 02:00:00,2.125117,1.202160,1.000000,0.400000,0.440003
2013-05-20 03:00:00,0.000000,6.046332,0.000000,0.500000,0.436236
2013-05-20 04:00:00,0.000000,5.318681,0.000000,0.600000,0.435088
2013-05-20 05:00:00,0.000000,5.062000,0.000000,0.700000,0.434865
2013-05-20 06:00:00,0.000000,5.009117,0.000000,0.800000,0.434832
2013-05-20 07:00:00,0.000000,5.000202,0.000000,0.900000,0.434832
2013-05-20 08:00:00,0.000000,4.999319,0.000000,1.000000,0.434834
2013-05-20 09:00:00,0.000000,5.000742,0.000000,0.900000,0.434831
2013-05-20 10:00:00,0.000000,5.000108,0.000000,0.800000,0.434831
2013-05-20 11:00:00,0.000000,4.992524,0.000000,0.700000,0.434858
2013-05-20 12:00:00,0.000000,4.944152,0.000000,0.600000,0.435059
2013-05-20 13:00:00,0.000000,4.622809,0.000000,0.500000,0.436417
2013-05-20 14:00:00,2.125114,1.878932,1.000000,0.400000,0.440002
2013-05-20 15:00:00,3.975640,1.024350,1.000000,0.300000,0.440002
2013-05-20 16:00:00,5.105053,0.000000,1.000000,0.200000,0.439624
2013-05-20 17:00:00,5.543917,0.000000,1.000000,0.100000,0.437666
2013-05-20 18:00:00,5.286357,0.000000,1.000000,0.000000,0.436635
2013-05-20 00:00:00,1.669620,0.000000,1.000000,0.200000,0.429989
2013-05-20 01:00:00,3.889305,0.000000,1.000000,0.300000,0.433988
2013-05-20 02:00:00,2.125114,1.204191,1.000000,0.400000,0.440002
2013-05-20 03:00:00,0.000000,5.040872,0.000000,0.500000,0.439855
2013-05-20 04:00:00,0.000000,5.422016,0.000000,0.600000,0.438336
2013-05-20 05:00:00,0.000000,5.426202,0.000000,0.700000,0.436802
2013-05-20 06:00:00,0.000000,5.341535,0.000000,0.800000,0.435572
2013-05-20 07:00:00,0.000000,5.264185,0.000000,0.900000,0.434621
2013-05-20 08:00:00,0.000000,5.221975,0.000000,1.000000,0.433822
2013-05-20 09:00:00,0.000000,5.202363,0.000000,0.900000,0.433093
2013-05-20 10:00:00,0.000000,5.156304,0.000000,0.800000,0.432531
2013-05-20 11:00:00,0.000000,4.994088,0.000000,0.700000,0.432552
2013-05-20 12:00:00,0.000000,4.576840,0.000000,0.600000,0.434075
2013-05-20 13:00:00,0.000000,3.755049,0.000000,0.500000,0.438557
2013-05-20 14:00:00,2.125106,2.473546,1.000000,0.400000,0.440002
2013-05-20 15:00:00,3.975632,1.024423,1.000000,0.300000,0.440002
2013-05-20 16:00:00,5.105254,0.000000,1.000000,0.200000,0.439623
2013-05-20 17:00:00,5.544216,0.000000,1.000000,0.100000,0.437664
2013-05-20 18:00:00,5.286836,0.000000,1.000000,0.000000,0.436631
......@@ -34,10 +34,8 @@ class MinimizeQpumpGoal(Goal):
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 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
......@@ -52,7 +50,6 @@ class MinimizeChangeInQpumpGoal(Goal):
# 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
......@@ -77,7 +74,7 @@ class Example(GoalProgrammingMixin, CSVMixin, ModelicaMixin,
# 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
M = 2 # The so-called "big-M"
constraints.append((self.state('H_sea') - self.state('storage.HQ.H') -
(1 - self.state('is_downhill')) * M, -np.inf, 0.0))
constraints.append((self.state('H_sea') - self.state('storage.HQ.H') +
......
......@@ -17,11 +17,10 @@ class WaterVolumeRangeGoal(StateGoal):
# 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().__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')
super().__init__(optimization_problem)
state = 'storage.V'
priority = 1
......
......@@ -4,19 +4,19 @@ time,Q_orifice,Q_pump,is_downhill,sea_level,storage_level
2013-05-20 00:00:00,4.729408,0.000000,1.000000,0.200000,0.398123
2013-05-20 01:00:00,3.423235,0.000000,1.000000,0.300000,0.403799
2013-05-20 02:00:00,1.378672,0.000000,1.000000,0.400000,0.416836
2013-05-20 03:00:00,0.000000,4.476889,0.000000,0.500000,0.418719
2013-05-20 04:00:00,0.000000,4.103700,0.000000,0.600000,0.421946
2013-05-20 05:00:00,0.000000,3.817639,0.000000,0.700000,0.426203
2013-05-20 06:00:00,0.000000,3.544285,0.000000,0.800000,0.431443
2013-05-20 07:00:00,0.000000,3.265951,0.000000,0.900000,0.437686
2013-05-20 08:00:00,0.000000,2.976000,0.000000,1.000000,0.444972
2013-05-20 09:00:00,0.000000,2.671708,0.000000,0.900000,0.453354
2013-05-20 10:00:00,0.000000,2.340619,0.000000,0.800000,0.462928
2013-05-20 11:00:00,0.000000,1.983423,0.000000,0.700000,0.473787
2013-05-20 12:00:00,0.000000,1.592914,0.000000,0.600000,0.486053
2013-05-20 13:00:00,0.000315,1.125490,1.000000,0.500000,0.500000
2013-05-20 03:00:00,0.000000,4.403457,0.000000,0.500000,0.418984
2013-05-20 04:00:00,0.000000,4.039766,0.000000,0.600000,0.422441
2013-05-20 05:00:00,0.000000,3.770278,0.000000,0.700000,0.426868
2013-05-20 06:00:00,0.000000,3.517466,0.000000,0.800000,0.432205
2013-05-20 07:00:00,0.000000,3.256937,0.000000,0.900000,0.438480
2013-05-20 08:00:00,0.000000,2.994531,0.000000,1.000000,0.445699
2013-05-20 09:00:00,0.000000,2.712544,0.000000,0.900000,0.453934
2013-05-20 10:00:00,0.000000,2.390271,0.000000,0.800000,0.463329
2013-05-20 11:00:00,0.000000,2.029730,0.000000,0.700000,0.474022
2013-05-20 12:00:00,0.000000,1.633157,0.000000,0.600000,0.486143
2013-05-20 13:00:00,0.001503,1.149293,1.000000,0.500000,0.500000
2013-05-20 14:00:00,3.360000,1.640000,1.000000,0.400000,0.500000
2013-05-20 15:00:00,4.751758,0.248242,1.000000,0.300000,0.500000
2013-05-20 16:00:00,5.475820,0.000000,1.000000,0.200000,0.498287
2013-05-20 17:00:00,5.815266,0.000000,1.000000,0.100000,0.495352
2013-05-20 18:00:00,5.494218,0.000000,1.000000,0.000000,0.493573
2013-05-20 16:00:00,5.485497,0.000000,1.000000,0.200000,0.498252
2013-05-20 17:00:00,5.819139,0.000000,1.000000,0.100000,0.495303
2013-05-20 18:00:00,5.505331,0.000000,1.000000,0.000000,0.493484
......@@ -26,8 +26,7 @@ class Example(CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):
def path_constraints(self, ensemble_member):
# Call super to get default constraints
constraints = super().path_constraints(ensemble_member)
# M is a handy big number
M = 1e10
M = 2 # The so-called "big-M"
# Release through orifice downhill only. This constraint enforces the
# fact that water only flows downhill.
......
......@@ -41,10 +41,10 @@ setup(
platforms = ['Windows', 'Linux', 'Mac OS-X', 'Unix'],
packages = find_packages("src"),
package_dir = {"": "src"},
install_requires = ["casadi >= 3.2.0",
install_requires = ["casadi == 3.4.*",
"numpy >= 1.14.0",
"scipy >= 1.0.0",
"pymoca >= 0.2.6"],
"pymoca == 0.3.*"],
tests_require = ['nose'],
test_suite = 'nose.collector',
python_requires='>=3.5',
......
......@@ -44,7 +44,7 @@ def reduce_matvec(e, v):
def substitute_in_external(expr, symbols, values):
if len(symbols) == 0:
if len(symbols) == 0 or all(isinstance(x, ca.DM) for x in expr):
return expr
else:
f = ca.Function('f', symbols, expr)
......
......@@ -642,10 +642,13 @@ class Timeseries:
raise Exception('Unsupported unit type: ' + el.get('unit'))
def __floor_date_time(self, dt, tdel):
# Floor a PI date time based on a PI time step
"""
Floor a PI date time to an integer number of PI time steps from the
start date time.
"""
roundTo = tdel.total_seconds()
seconds = (dt - dt.min).seconds
seconds = (dt - self.__start_datetime).seconds
# // is a floor division:
rounding = (seconds + roundTo / 2) // roundTo * roundTo
return dt + datetime.timedelta(0, rounding - seconds, -dt.microsecond)
......@@ -1040,4 +1043,4 @@ class Timeseries:
ensemble member.
"""
for key in self.__values[ensemble_member].keys():
yield key, self.get(key)
yield key, self.get(key, ensemble_member)
......@@ -232,13 +232,13 @@ class ControlTreeMixin(OptimizationProblem):
def extract_controls(self, ensemble_member=0):
# Solver output
X = self.solver_output
X = self.solver_output.copy()
# Extract control inputs
results = {}
for variable in self.controls:
results[variable] = np.array(self.variable_nominal(
variable) * X[self.__control_indices[ensemble_member][variable], 0]).ravel()
results[variable] = self.variable_nominal(
variable) * X[self.__control_indices[ensemble_member][variable]].ravel()
# Done
return results
......
......@@ -3,7 +3,7 @@ import glob
import logging
import os
import pickle
from typing import List, Tuple, Union
from typing import Iterable, List, Tuple, Union
import casadi as ca
......@@ -28,7 +28,7 @@ class LookupTable:
Lookup table.
"""
def __init__(self, inputs: List[ca.MX], function: ca.Function, tck: Tuple=None):
def __init__(self, inputs: List[ca.MX], function: ca.Function, tck: Tuple = None):
"""
Create a new lookup table object.
......@@ -79,12 +79,23 @@ class LookupTable:
"""
return self.__function
def __call__(self, *args: List[Union[float, Timeseries]]) -> Union[float, Timeseries]:
@property
@cached
def __numeric_function_evaluator(self):
return np.vectorize(
lambda *args: np.nan
if np.any(np.isnan(args))
else np.float(self.function(*args))
)
def __call__(
self, *args: Union[float, Iterable, Timeseries]
) -> Union[float, np.ndarray, Timeseries]:
"""
Evaluate the lookup table.
:param args: Input values.
:type args: Float, iterable of floats, or :class:`Timeseries`
:type args: Float, iterable of floats, or :class:`.Timeseries`
:returns: Lookup table evaluated at input values.
Example use::
......@@ -93,42 +104,95 @@ class LookupTable:
[y1, y2] = lookup_table([1.0, 2.0])
"""
if isinstance(args[0], Timeseries):
return Timeseries(args[0].times, self(args[0].values))
evaluator = self.__numeric_function_evaluator
if len(args) == 1:
arg = args[0]
if isinstance(arg, Timeseries):
return Timeseries(arg.times, self(arg.values))
else:
if hasattr(arg, "__iter__"):
arg = np.fromiter(arg, dtype=float)
return evaluator(arg)
else:
arg = float(arg)
return evaluator(arg).item()
else:
if hasattr(args[0], '__iter__'):
evaluator = np.vectorize(
lambda v: float(self.function(v)))
return evaluator(args[0])
if any(isinstance(arg, Timeseries) for arg in args):
raise TypeError(
"Higher-order LookupTable calls do not yet support Timeseries parameters"
)
elif any(hasattr(arg, "__iter__") for arg in args):
raise TypeError(
"Higher-order LookupTable calls do not yet support vector parameters"
)
else:
return float(self.function(*args))
def reverse_call(self, y, domain=(None, None), detect_range_error=True):
"""
use scipy brentq optimizer to do an inverted call to this lookuptable
args = np.fromiter(args, dtype=float)
return evaluator(*args)
def reverse_call(
self,
y: Union[float, Iterable, Timeseries],
domain: Tuple[float, float] = (None, None),
detect_range_error: bool = True,
) -> Union[float, np.ndarray, Timeseries]:
"""Do an inverted call on this LookupTable
Uses SciPy brentq optimizer to simulate a reversed call.
Note: Method does not work with higher-order LookupTables
"""
if isinstance(y, Timeseries):
# Recurse and return
return Timeseries(y.times, self.reverse_call(y.values))
# Get domain information
l_d, u_d = domain
if l_d is None:
l_d = self.domain[0]
if u_d is None:
u_d = self.domain[1]
if detect_range_error:
l_r, u_r = self.range
# Cast y to array of float
if hasattr(y, "__iter__"):
y_array = np.fromiter(y, dtype=float)
else:
y_array = np.array([y], dtype=float)
def function(y_target):
if detect_range_error and (y_target < l_r or y_target > u_r):
raise ValueError('Value {} is not in lookup table range ({}, {})'.format(y_target, l_r, u_r))
return brentq(lambda x: self(x) - y_target, l_d, u_d)
# Find not np.nan
is_not_nan = ~np.isnan(y_array)
y_array_not_nan = y_array[is_not_nan]
if isinstance(y, Timeseries):
return Timeseries(y.times, self.reverse_call(y.values))
# Detect if there is a range violation
if detect_range_error:
l_r, u_r = self.range
lb_viol = y_array_not_nan < l_r
ub_viol = y_array_not_nan > u_r
all_viol = y_array_not_nan[lb_viol | ub_viol]
if all_viol:
raise ValueError(
"Values {} are not in lookup table range ({}, {})".format(
all_viol, l_r, u_r
)
)
# Construct function to do inverse evaluation
evaluator = self.__numeric_function_evaluator
def inv_evaluator(y_target):
"""inverse evaluator function"""
return brentq(lambda x: evaluator(x) - y_target, l_d, u_d)
inv_evaluator = np.vectorize(inv_evaluator)
# Calculate x_array
x_array = np.full_like(y_array, np.nan, dtype=float)
if y_array_not_nan.size != 0:
x_array[is_not_nan] = inv_evaluator(y_array_not_nan)
# Return x
if hasattr(y, "__iter__"):
return x_array
else:
if hasattr(y, '__iter__'):
evaluator = np.vectorize(function)
return evaluator(y)
else:
return function(y)
return x_array.item()
class CSVLookupTableMixin(OptimizationProblem):
......
......@@ -169,7 +169,7 @@ class CSVMixin(OptimizationProblem):
if self.__timeseries_times_sec[i + 1] - self.__timeseries_times_sec[i] != dt:
raise Exception(
'CSVMixin: Expecting equidistant timeseries, the time step towards '
'{} is not the same as the time step(s) before. Set equidistant=False '
'{} is not the same as the time step(s) before. Set csv_equidistant = False '
'if this is intended.'.format(self.__timeseries_times[i + 1]))
def times(self, variable=None):
......@@ -324,8 +324,8 @@ class CSVMixin(OptimizationProblem):
names = ['time'] + sorted({sym.name() for sym in self.output_variables})
formats = ['O'] + (len(names) - 1) * ['f8']
dtype = {'names': names, 'formats': formats}
data = np.zeros(len(self.__timeseries_times), dtype=dtype)
data['time'] = self.__timeseries_times
data = np.zeros(len(times), dtype=dtype)
data['time'] = [self.__timeseries_times[0] + timedelta(seconds=s) for s in times]
for output_variable in self.output_variables:
output_variable = output_variable.name()
try:
......
......@@ -30,17 +30,28 @@ class HomotopyMixin(OptimizationProblem):
# Do not override any previously seeded values, such as goal programming results.
for key, result in self.__results[ensemble_member].items():
times = self.times(key)
if key not in seed and len(result) == len(times):
if key in seed:
continue
elif len(result) == len(times):
# Only include seed timeseries which are consistent
# with the specified time stamps.
seed[key] = Timeseries(times, result)
elif len(result) == 1:
seed[key] = result
return seed
def parameters(self, ensemble_member):
parameters = super().parameters(ensemble_member)
options = self.homotopy_options()
parameters[options['homotopy_parameter']] = self.__theta
try:
# Only set the theta if we are in the optimization loop. We want
# to avoid accidental usage of the parameter value in e.g. pre().
# Note that we use a try-except here instead of hasattr, to avoid
# explicit name mangling.
parameters[options['homotopy_parameter']] = self.__theta
except AttributeError:
pass
return parameters
......
......@@ -8,8 +8,7 @@ class _MeasurementGoal(Goal):
self.__state = state
self.__measurement_id = measurement_id
self.function_range = -max_deviation, max_deviation
self.function_nominal = 1.0
self.function_nominal = max_deviation
def function(self, optimization_problem, ensemble_member):
op = optimization_problem
......@@ -26,8 +25,7 @@ class _SmoothingGoal(Goal):
self.__state1 = state1
self.__state2 = state2
self.function_range = -max_deviation, max_deviation
self.function_nominal = 1.0
self.function_nominal = max_deviation
def function(self, optimization_problem, ensemble_member):
op = optimization_problem
......
......@@ -215,11 +215,19 @@ class ModelicaMixin(OptimizationProblem):
constant_input_names = {sym.name() for sym in self.__mx['constant_inputs']}
for v in self.__pymoca_model.inputs:
if v.symbol.name() in constant_input_names:
constant_inputs[v.symbol.name()] = Timeseries(
times, np.full_like(times, v.value))
if logger.getEffectiveLevel() == logging.DEBUG:
logger.debug("Read constant input {} = {} from Modelica model".format(
v.symbol.name(), v.value))
# NOTE: v.value can be a DM. Cast to float to make sure the
# isnan call succeeds.
value = float(v.value)
if not np.isnan(value):
constant_inputs[v.symbol.name()] = Timeseries(
times, np.full_like(times, value))
if logger.getEffectiveLevel() == logging.DEBUG:
logger.debug(
"Read constant input {} = {} from Modelica model".format(
v.symbol.name(), value
)
)
return constant_inputs
......@@ -354,6 +362,7 @@ class ModelicaMixin(OptimizationProblem):
# If nominal contains parameter symbols, substitute them
if not nominal.is_constant():
[nominal] = substitute_in_external([nominal], self.__mx['parameters'], parameter_values)
nominal = ca.MX(nominal)
if nominal.is_constant():
# Take absolute value (nominal sign is meaningless- a nominal is a magnitude)
......
......@@ -36,7 +36,7 @@ class PIMixin(OptimizationProblem):
:cvar pi_parameter_config_numerical_basename:
Numerical config file basename to read. Default is ``rtcParameterConfig_Numerical``.
:cvar pi_check_for_duplicate_parameters:
Check if duplicate parameters are read. Default is ``False``.
Check if duplicate parameters are read. Default is ``True``.
:cvar pi_validate_timeseries:
Check consistency of timeseries. Default is ``True``.
"""
......@@ -171,6 +171,9 @@ class PIMixin(OptimizationProblem):
# Call parent class first for default values.
parameters = super().parameters(ensemble_member)
# Make a new dictionary for the new parameters
new_parameters = AliasDict(self.alias_relation)
# Load parameters from parameter config
for parameter_config in self.__parameter_config:
for location_id, model_id, parameter_id, value in parameter_config:
......@@ -180,13 +183,16 @@ class PIMixin(OptimizationProblem):
parameter = parameter_id
if self.pi_check_for_duplicate_parameters:
if parameter in parameters.keys():
if parameter in new_parameters.keys():
logger.warning(
'PIMixin: parameter {} defined in file {} was already '
'present. Using value {}.'.format(
'present in another or this parameterConfig file. Using value {}.'.format(
parameter, parameter_config.path, value))
parameters[parameter] = value
new_parameters[parameter] = value
# Update parameters
parameters.update(new_parameters)
# Done
return parameters
......@@ -317,21 +323,31 @@ class PIMixin(OptimizationProblem):
# Call parent class first for default behaviour.
super().post()
# Get time stamps
times = self.times()
if len(set(times[1:] - times[:-1])) == 1:
dt = timedelta(seconds=times[1] - times[0])
else:
dt = None
# Start of write output
# Write the time range for the export file.
self.__timeseries_export.times = self.__timeseries_import.times[self.__timeseries_import.forecast_index:]
self.__timeseries_export.times = [
self.__timeseries_import.times[self.__timeseries_import.forecast_index]
+ timedelta(seconds=s) for s in times]
# Write other time settings
self.__timeseries_export.forecast_datetime = self.__timeseries_import.forecast_datetime
self.__timeseries_export.dt = self.__timeseries_import.dt
self.__timeseries_export.dt = dt
self.__timeseries_export.timezone = self.__timeseries_import.timezone
# Write the ensemble properties for the export file.
if self.ensemble_size > 1:
self.__timeseries_export.contains_ensemble = True
self.__timeseries_export.ensemble_size = self.ensemble_size
self.__timeseries_export.contains_ensemble = self.__timeseries_import.contains_ensemble
# Start looping over the ensembles for extraction of the output values.
times = self.times()
for ensemble_member in range(self.ensemble_size):
results = self.extract_results(ensemble_member)
......
......@@ -29,7 +29,7 @@ class PIMixin(SimulationProblem):
:cvar pi_binary_timeseries: Whether to use PI binary timeseries format. Default is ``False``.
:cvar pi_parameter_config_basenames:
List of parameter config file basenames to read. Default is [``rtcParameterConfig``].
:cvar pi_check_for_duplicate_parameters: Check if duplicate parameters are read. Default is ``False``.
:cvar pi_check_for_duplicate_parameters: Check if duplicate parameters are read. Default is ``True``.
:cvar pi_validate_timeseries: Check consistency of timeseries. Default is ``True``.
"""
......
......@@ -32,6 +32,10 @@ class SimulationProblem:
# Folders in which the referenced Modelica libraries are found
modelica_library_folders = []
# Force workaround for delay support by assuming zero delay. This flag
# will be removed when proper delay support is added.
_force_zero_delay = False
def __init__(self, **kwargs):
# Check arguments
assert('model_folder' in kwargs)
......@@ -62,11 +66,8 @@ class SimulationProblem:
self.__mx['constant_inputs'] = []
self.__mx['lookup_tables'] = []
# TODO: implement delayed feedback
delayed_feedback_variables = []
for v in self.__pymoca_model.inputs:
if v.symbol.name() in delayed_feedback_variables:
if v.symbol.name() in self.__pymoca_model.delay_states:
# Delayed feedback variables are local to each ensemble, and
# therefore belong to the collection of algebraic variables,
# rather than to the control inputs.
......@@ -147,8 +148,19 @@ class SimulationProblem:
for index, derivative_state in enumerate(self.__mx['derivatives']):
derivative_approximation_residuals.append(derivative_state - (X[index] - X_prev[index]) / dt)
if self.__pymoca_model.delay_states and not self._force_zero_delay:
raise NotImplementedError("Delayed states are not supported")
# Delayed feedback (assuming zero delay)
# TODO: implement delayed feedback support for delay != 0
delayed_feedback_equations = []
for delay_state, delay_argument in zip(self.__pymoca_model.delay_states,
self.__pymoca_model.delay_arguments):
logger.warning("Assuming zero delay for delay state '{}'".format(delay_state))
delayed_feedback_equations.append(delay_argument.expr - self.__sym_dict[delay_state])
# Append residuals for derivative approximations
dae_residual = ca.vertcat(self.__dae_residual, *derivative_approximation_residuals)
dae_residual = ca.vertcat(self.__dae_residual, *derivative_approximation_residuals, *delayed_feedback_equations)
# TODO: implement lookup_tables
......@@ -224,12 +236,30 @@ class SimulationProblem:
logger.debug('SimulationProblem: Setting parameter {} = {}'.format(var.symbol.name(), val))
self.set_var(var.symbol.name(), val)
# Nominals can be symbolic, written in terms of parameters. After all
# parameter values are known, we evaluate the numeric values of the
# nominals.
nominal_vars = list(itertools.chain(
self.__pymoca_model.states, self.__pymoca_model.alg_states, self.__pymoca_model.der_states))
symbolic_nominals = ca.vertcat(*[self.get_variable_nominal(v.symbol.name()) for v in nominal_vars])
nominal_evaluator = ca.Function('nominal_evaluator', self.__mx['parameters'], [symbolic_nominals])
n_parameters = len(self.__mx['parameters'])
if n_parameters > 0:
[evaluated_nominals] = nominal_evaluator.call(self.__state_vector[-n_parameters:])
else:
[evaluated_nominals] = nominal_evaluator.call([])
evaluated_nominals = np.array(evaluated_nominals).ravel()
nominal_dict = {v.symbol.name(): n for v, n in zip(nominal_vars, evaluated_nominals)}
# Assemble initial residuals and set values from start attributes into the state vector
constrained_residuals = []
minimized_residuals = []
for var in itertools.chain(self.__pymoca_model.states, self.__pymoca_model.alg_states):
var_name = var.symbol.name()
var_nominal = self.get_variable_nominal(var_name)
var_nominal = nominal_dict[var_name]
# Attempt to cast var.start to python type
mx_start = ca.MX(var.start)
......@@ -298,7 +328,11 @@ class SimulationProblem:
# Make a list of unscaled symbols and a list of their scaled equivalent
unscaled_symbols = []
scaled_symbols = []
for sym_name, nominal in self.__nominals.items():
for sym_name, nominal in nominal_dict.items():
if nominal == 1.0:
# Nothing to scale
continue
# Add the symbol to the lists
symbol = self.__sym_dict[sym_name]
unscaled_symbols.append(symbol)
......@@ -317,8 +351,9 @@ class SimulationProblem:
# State bounds can be symbolic, written in terms of parameters. After all
# parameter values are known, we evaluate the numeric values of bounds.
symbolic_bounds = ca.vertcat(*[ca.horzcat(v.min, v.max) for v in itertools.chain(
self.__pymoca_model.states, self.__pymoca_model.alg_states, self.__pymoca_model.der_states)])
bound_vars = list(itertools.chain(
self.__pymoca_model.states, self.__pymoca_model.alg_states, self.__pymoca_model.der_states))
symbolic_bounds = ca.vertcat(*[ca.horzcat(v.min, v.max) for v in bound_vars])
bound_evaluator = ca.Function('bound_evaluator', self.__mx['parameters'], [symbolic_bounds])
# Evaluate bounds using values of parameters
......@@ -328,6 +363,21 @@ class SimulationProblem:
else:
[evaluated_bounds] = bound_evaluator.call([])
# Scale the bounds with the nominals
nominals = []
for var in bound_vars:
nominals.append(nominal_dict[var.symbol.name()])
evaluated_bounds = np.array(evaluated_bounds) / np.array(nominals)[:,