Dynamic simulation-related updates

biosteam/_system.py

@@ -47,7 +47,7 @@
 import numpy as np
 import pandas as pd
 from numpy.linalg import solve
-from scipy.integrate import solve_ivp
+from scipy.integrate import solve_ivp, RK23, RK45, DOP853, Radau, BDF, LSODA
 from . import report
 from thermosteam.network import temporary_units_dump, TemporaryUnit
 import os
@@ -3281,19 +3281,6 @@ def _write_state(self):
         for ws in [i for i in self.streams if i not in self.feeds]:
             ws._state2flows()
 
-    def clear_state(self):
-        """Clear all states and dstates (system, units, and streams)."""
-        self._state = None
-        for u in self.units:
-            u._state = None
-            u._dstate = None
-        for ws in self.streams:
-            y = ws.state.copy()
-            ws._state = None
-            ws._dstate = None
-            ws.state = y*0.0
-            ws.dstate = y*0.0
-
     def mock_convergence(
             self, 
             convergence_rate=10, 
@@ -3457,6 +3444,118 @@ def simulate(self, update_configuration: Optional[bool]=None, units=None,
         if save_configuration: self._save_configuration()
         return outputs
 
+    _dynsim_special_keys = frozenset((
+        't_eval', 'state_reset_hook', 'export_state_to', 'sample_id',
+        'print_msg', 'print_t', 'timeout', 'y0',
+    ))
+    _ivp_solvers = {
+        'RK23': RK23, 'RK45': RK45, 'DOP853': DOP853,
+        'Radau': Radau, 'BDF': BDF, 'LSODA': LSODA,
+    }
+
+    @staticmethod
+    def _solver_options(solver):
+        """Return the names of the options accepted by a solver class, or None if any are."""
+        parameters = signature(solver.__init__).parameters
+        # scipy's own solvers swallow unknown options in `extraneous` and warn once
+        if any(i.kind is i.VAR_KEYWORD and i.name != 'extraneous' for i in parameters.values()): return None
+        return set(parameters) - {'self', 'fun', 't0', 'y0', 't_bound', 'vectorized', 'extraneous'}
+
+    def _check_dynsim_kwargs(self, kwargs):
+        """
+        Validate dynamic simulation keyword arguments. Raise an error if a
+        value is invalid or an option does not exist for any solver.
+        Return the set of keys that only apply to other solvers, which should
+        be dropped (a warning is raised).
+        """
+        if 't_span' in kwargs:
+            t_span = kwargs['t_span']
+            try: t_span = tuple(t_span)
+            except TypeError:
+                raise TypeError(
+                    f't_span must be a (t0, tf) pair of numbers; not {t_span!r}'
+                ) from None
+            if len(t_span) != 2:
+                raise ValueError(f't_span must be a (t0, tf) pair of numbers; got {t_span!r}')
+            try: t0, tf = map(float, t_span)
+            except (TypeError, ValueError):
+                raise TypeError(f't_span must be a (t0, tf) pair of numbers; got {t_span!r}') from None
+            if tf < t0:
+                raise ValueError(
+                    f't_span must be increasing; got t0={t0:g} > tf={tf:g}. '
+                    'Backwards integration is not supported for dynamic systems'
+                )
+            t_eval = kwargs.get('t_eval')
+            if t_eval is not None and kwargs.get('export_state_to'):
+                try: t_eval = np.asarray(t_eval, dtype=float)
+                except (TypeError, ValueError):
+                    raise TypeError(f't_eval must be an iterable of numbers; got {t_eval!r}') from None
+                if t_eval.size and not (t0 <= t_eval.min() and t_eval.max() <= tf):
+                    raise ValueError(
+                        f't_eval must be within t_span={t_span!r} to export the states; '
+                        f'got values in [{t_eval.min():g}, {t_eval.max():g}]'
+                    )
+        hook = kwargs.get('state_reset_hook')
+        if hook:
+            if isinstance(hook, str):
+                if not callable(getattr(self, hook, None)):
+                    raise AttributeError(
+                        f"state_reset_hook {hook!r} is not a method of {self!r}; "
+                        "use 'reset_cache' or a callable"
+                    )
+            elif not callable(hook):
+                raise TypeError(
+                    "state_reset_hook must be 'reset_cache', "
+                    f"a callable, or None; not {hook!r}"
+                )
+        rtol = kwargs.get('rtol')
+        if rtol is not None and np.ndim(rtol) == 0:
+            try: valid = float(rtol) > 0
+            except (TypeError, ValueError): valid = False
+            if not valid: raise ValueError(f'rtol must be a positive number; got {rtol!r}')
+        atol = kwargs.get('atol')
+        if atol is not None:
+            try: valid = bool(np.all(np.asarray(atol, dtype=float) > 0))
+            except (TypeError, ValueError): valid = False
+            if not valid: raise ValueError(f'atol must be positive; got {atol!r}')
+        for key, finite in (('first_step', True), ('max_step', False)):
+            value = kwargs.get(key)
+            if value is None: continue
+            try: valid = 0 < float(value) and (not finite or np.isfinite(float(value)))
+            except (TypeError, ValueError): valid = False
+            if not valid:
+                raise ValueError(
+                    f'{key} must be a positive{" finite" if finite else ""} number; '
+                    f'got {value!r}'
+                )
+        # Options that the chosen solver would silently ignore
+        method = kwargs.get('method', 'RK45')
+        solver = self._ivp_solvers.get(method) if isinstance(method, str) else method
+        if solver is None: return set() # solve_ivp reports the invalid method
+        try: options = self._solver_options(solver)
+        except (TypeError, ValueError, AttributeError): return set() # not a solver class
+        if options is None: return set()
+        known = options | self._dynsim_special_keys | {'t_span', 'method', 'dense_output', 'events', 'args', 'vectorized'}
+        others = set().union(*[self._solver_options(i) for i in self._ivp_solvers.values()])
+        ignored = set(); unknown = []
+        for key in kwargs:
+            if key in known: continue
+            if key in others: ignored.add(key)
+            else: unknown.append(key)
+        if unknown:
+            from difflib import get_close_matches
+            msgs = []
+            for key in unknown:
+                close = get_close_matches(key, known, n=1)
+                msgs.append(repr(key) + (f' (did you mean {close[0]!r}?)' if close else ''))
+            raise TypeError(
+                f'unexpected dynamic simulation keyword argument(s): {", ".join(msgs)}'
+            )
+        if ignored:
+            warn(f'option(s) {sorted(ignored)} have no effect for solver method '
+                 f'{method!r} and are ignored', RuntimeWarning, stacklevel=3)
+        return ignored
+
     def dynamic_run(self, **dynsim_kwargs):
         """
         Run system dynamically without accounting for
@@ -3471,15 +3570,25 @@ def dynamic_run(self, **dynsim_kwargs):
                     Interval of integration (t0, tf).
                     The solver starts with t=t0 and integrates until it reaches t=tf.
                 t_eval : iterable(float)
-                    The time points where status will be saved.
+                    The time points at which the tracked states are interpolated
+                    when exporting with ``export_state_to``. It does not change
+                    which points are recorded in ``System.scope`` (every solver
+                    step is recorded) and is not passed to ``solve_ivp``.
                 state_reset_hook: str|Callable
-                    Hook function to reset the cache state between simulations
-                    for dynamic systems).
-                    Can be "reset_cache" or "clear_state" to call `System.reset_cache`
-                    or `System.clear_state`, or None to avoiding resetting.
+                    What to do with the cached state of a dynamic system
+                    before the run starts. One of:
+
+                    - ``"reset_cache"``: call `System.reset_cache`, so the run
+                      starts again from the initial conditions.
+                    - A function that takes no arguments: it is called before
+                      the run, so you can reset the state yourself.
+                    - ``None`` (default): do not reset, so the run continues
+                      from the current state (e.g., to continue after a spin-up).
                 export_state_to: str
                     If provided with a path, will save the simulated states over time to the given path,
-                    supported extensions are ".xlsx", ".xls", "csv", and "tsv".
+                    supported extensions are ".npy", ".xlsx", ".xls", ".csv", and ".tsv".
+                    A path without an extension is saved as ".npy", and so is
+                    any path when ``sample_id`` is given.
                 sample_id : str
                     ID of the samples to run (for results exporting).
                 print_msg : bool
@@ -3488,15 +3597,46 @@ def dynamic_run(self, **dynsim_kwargs):
                     Whether to print integration time in the console,
                     usually used for debugging.
                 solve_ivp_kwargs
-                    All remaining keyword arguments will be passed to ``solve_ivp``.
-
+                    All remaining keyword arguments will be passed to ``solve_ivp``,
+                    which includes the options of the chosen ``method`` (e.g.,
+                    ``rtol``, ``atol``, ``max_step``). Options that do not exist
+                    for any solver raise a ``TypeError``; options that exist only
+                    for other solvers are dropped with a warning.
+                    ``y0`` is not accepted: the initial state is always taken from
+                    the system's current state (a warning is raised if it is passed).
+
+        Notes
+        -----
+        Keyword arguments are validated before they are stored in
+        ``System.dynsim_kwargs``, so an invalid value does not persist
+        to later runs.
+
+        A ``RuntimeError`` is raised if the solver fails (for example, the step
+        size becomes too small), instead of going on to write the states,
+        export, and run design and cost. The partial results remain in
+        ``System.scope``. A run stopped by an event is not a failure.
+
         See Also
         --------
         `scipy.integrate.solve_ivp <https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html>`_
 
         """
+        if 't_span' not in dynsim_kwargs and 't_span' not in self.dynsim_kwargs:
+            # `dynsim_kwargs` persists between runs, so a previous `t_span` is reused
+            raise TypeError(
+                f'{self!r} is dynamic and needs a time span; pass `t_span`, '
+                'for example `system.simulate(t_span=(0, 10))`. For a '
+                'steady-state simulation, set `system.isdynamic = False`.'
+            )
+        if 'y0' in dynsim_kwargs:
+            warn("`y0` is ignored; the initial state is taken from the system's "
+                 "current state", RuntimeWarning, stacklevel=2)
+            dynsim_kwargs = {k: v for k, v in dynsim_kwargs.items() if k != 'y0'}
+        # Validate before storing, so that an invalid value does not persist
+        ignored = self._check_dynsim_kwargs({**self.dynsim_kwargs, **dynsim_kwargs})
         dk = self.dynsim_kwargs
         dk.update(dynsim_kwargs)
+        for key in ignored: dk.pop(key) # do not keep options that this solver does not use
         dk_cp = dk.copy()
         t_eval = dk_cp.pop('t_eval', None)
         state_reset_hook = dk_cp.pop('state_reset_hook', None)
@@ -3517,7 +3657,10 @@ def dynamic_run(self, **dynsim_kwargs):
             for u in self.units:
                 if not hasattr(u, '_state'): u._init_dynamic()
         # Load initial states
-        self.converge()
+        # Streams only matter for the initial state; once a state is loaded
+        # (continuing a previous run), the static solve is not needed and can
+        # fail from the stream compositions left by the previous run
+        if self._state is None: self.converge()
         y0, idx, nr = self._load_state()
         self.dynsim_kwargs['y0'] = y0.copy()
         self.dynsim_kwargs['print_t'] = print_t # self.dynsim_kwargs might be reset by `state_reset_hook`
@@ -3550,12 +3693,19 @@ def DAE(t, y, _f=self.DAE, _deadline=deadline, _timeout=timeout):
                     )
                 return _f(t, y)
         with np.errstate(invalid='ignore', divide='ignore'):
-            self.scope.sol = sol = solve_ivp(fun=DAE, y0=y0, **dk_cp)
+            # Pass a copy: `y0` is `self._state`, which `DAE` overwrites in place.
+            self.scope.sol = sol = solve_ivp(fun=DAE, y0=y0.copy(), **dk_cp)
         if print_msg:
             if sol.status == 0:
                 print('Simulation completed.')
             else: print(sol.message)
-        self._write_state()          
+        if sol.status < 0:
+            # Do not write the state or go on to design and cost with a failed solve
+            raise RuntimeError(
+                f'dynamic simulation of {self!r} failed at t={sol.t[-1]:g} '
+                f'(solver status {sol.status}): {sol.message}'
+            )
+        self._write_state()
         # Write states to file
         if export_state_to:
             try: file, ext = export_state_to.rsplit('.', 1)

tests/test_model.py

@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 # BioSTEAM: The Biorefinery Simulation and Techno-Economic Analysis Modules
 # Copyright (C) 2020-, Yoel Cortes-Pena <yoelcortes@gmail.com>, 
-#                      Yalin Li <zoe.yalin.li@gmail.com>,
+#                      Yalin Li <mailto.yalin.li@gmail.com>,
 #                      Sarang Bhagwat <sarangb2@gmail.com>
 # 
 # This module is under the UIUC open-source license. See 

tests/test_prediction.py

@@ -1,6 +1,7 @@
 # -*- coding: utf-8 -*-
 # BioSTEAM: The Biorefinery Simulation and Techno-Economic Analysis Modules
-# Copyright (C) 2020-2023, Yoel Cortes-Pena <yoelcortes@gmail.com>, Yalin Li <zoe.yalin.li@gmail.com>
+# Copyright (C) 2020-2023, Yoel Cortes-Pena <yoelcortes@gmail.com>,
+#                          Yalin Li <mailto.yalin.li@gmail.com>
 # 
 # This module is under the UIUC open-source license. See 
 # github.com/BioSTEAMDevelopmentGroup/biosteam/blob/master/LICENSE.txt