diff --git a/biosteam/_system.py b/biosteam/_system.py index df4c12c8..4031f494 100644 --- a/biosteam/_system.py +++ b/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 `_ """ + 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) diff --git a/tests/test_model.py b/tests/test_model.py index 78d366fd..998f9b8b 100644 --- a/tests/test_model.py +++ b/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 , -# Yalin Li , +# Yalin Li , # Sarang Bhagwat # # This module is under the UIUC open-source license. See diff --git a/tests/test_prediction.py b/tests/test_prediction.py index 91e16495..61de959e 100644 --- a/tests/test_prediction.py +++ b/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 , Yalin Li +# Copyright (C) 2020-2023, Yoel Cortes-Pena , +# Yalin Li # # This module is under the UIUC open-source license. See # github.com/BioSTEAMDevelopmentGroup/biosteam/blob/master/LICENSE.txt