Support for resuming simulations involving delay() and spatialDistribution()

[casella]

In this discussion post for OpenModelica, @mpradovelasco asks:

Is it possible to save and restore the history associated with the delay function so that the simulation of a model can be properly restored from a snapshot?

It would be really interesting to have this capability in pharmacometric and biosystems models, where the delay is widely used to calculate half-life, area under the curve, and other standard indices that support prescription optimization and genomic and personalized precision medicine.
Its implementation on openModelica will add a new issue where OM overcomes other object oriented Modeling and Simulation development environments like EcosimPro :)

I think this is not (only) a tool issue, but also involves at least two issues that are worth discussing here.

Issue 1: initial values for the delay() operator

The spatialDistribution() operator allows to specify the initial content of buffer that stores the internal state of the operator, so a tool could re-use this feature to initialize it with the result of a previous simulation.

Unfortunately, the delay() operator specification has no such provision, as it explicitly states:

Evaluates to 𝑒𝑥𝑝𝑟(time - 𝑑𝑒𝑙𝑎𝑦𝑇𝑖𝑚𝑒) for 𝚝𝚒𝚖𝚎>𝚝𝚒𝚖𝚎.𝚜𝚝𝚊𝚛𝚝+𝑑𝑒𝑙𝑎𝑦𝑇𝑖𝑚𝑒 and 𝑒𝑥𝑝𝑟(time.start) for 𝚝𝚒𝚖𝚎≤𝚝𝚒𝚖𝚎.𝚜𝚝𝚊𝚛𝚝+𝑑𝑒𝑙𝑎𝑦𝑇𝑖𝑚𝑒.

Should we add some additional optional argument to allow intializing the delay buffer memory as well? Of course we should retain the current behaviour when they are not specified. For example, we could add:

  input Real initialRelativeTimes = {0.0, -Modelica.Constants.Inf};
  input Real initialValues = {expr(time.start), expr(time.start)};

I think it would be best to use a relative time axis which is relative to the start time value, to make the operator behaviour more portable in case one want to start from a different point in time while retaining the same content of the delay buffer.

This would also come with the nice side effect that one could initialize models involving delay() in non-trivial initial conditions, e.g. a periodic steady-state.

Issue 2: how to write down these results in a result file

We should probably agree on a standardized way to represent these data in CSV files, so that they can be re-imported to continue a simulation.

[mpradovelasco]

I agree with @casella on the importance of making the operator portable to allow starting from a different point in time, as this could be a valuable feature when restoring a "snapshot".

I have no knowledge of the implementation of the delay() operator, although I feel that to avoid the hassle of changing the delayTime (if it is a parameter, page 24 in MLS3.6), the use of the optional third argument, delayMax, including the annotation (Evaluate = true) could be required to manage the init values of delay().

Thanks for considering the suggestion.

[HansOlsson]

For resuming simulations there are a number of issues:

• In practice we need to modify models, in order to ignore initial equations and some start-values (in some cases that can be done by modifiers but start-values may be final or in protected components, and there is no standardized way to disable initial equations). In practice we have found that this is a major obstacle. Tools such as Dymola can, of course, continue a simulation - but that is tool-specific and doesn't allow resuming it in a more general setting.
• Then there are these extra states - but compared to the above they are less problematic:
  • External objects would need to a save/restore state function; Introduce external object save/restore #1950
  • State machine would need parametrized initial states; Initial state in state machines #2826
  • Synchronous models would need special care for setting ticks in state etc.
  • And then delay/spatialDistribution
  • ...

If we just focus on tool-specific ways of resuming a simulation the blocker in the language is the lack of external object save/restore. And if we add that one could imagine that delay reuses on such a definition, i.e., each delay call is treated as an external object with a specified save/restore function.

[HansOlsson]

I would thus propose the following steps:

1. Allow tools to continue simulations (with tool internals). Required language change:
   • save/restore for External Objects; in particular for FMUs (as the underlying FMUs may have that functionality). Introduce external object save/restore #1950
2. Allow a standard way of continuing simulations. Required language changes:
   • some way of overriding start-values and initial equations even if the model wasn't prepared for it
   • including:
     • standard format for save/restore for External Objects
     • standard format for delay/spatialDistribution (similarly as External Object)
     • parametrized initial states in state-machines
     • ticks in state for synchronous models

For the second step we should likely ensure that you can also have this as a sub-model.

The idea in both of these steps is to be able to do it without changing the internals of the model, as the intended use case is that you have a model and want to be able to re-run it from some point in its trajectory.

Users just want it to work; and they don't want to spend time modifying the internals and investigating if that causes problems, especially as the model may be large, may contain parts made by others, and may contain encrypted parts.

That implies that #2826 or an extended delay-function as in this issue wouldn't help - as it requires modifying the model. Obviously we can still add such extensions, but we shouldn't consider it as a solution to this problem.

[mpradovelasco]

I fully agree on the difficulty of reaching such a milestone in a general case (and certainly in the steps described).

However, in the meantime, the possibility to access and somehow restore the history of relevant functions like delay could be an opportunity to develop specific libraries (it is not a general solution) that implement some particular functionalities. For example, the execution of parametric simulations or even Monte Carlo simulations can be implemented without the need for external code by developing a library of components prepared to save and restore their state under the control of an orchestrator within the same model. In this example, the lack of access to delay history prevents the inclusion of delay-based components.

[HansOlsson]

I fully agree on the difficulty of reaching such a milestone in a general case (and certainly in the steps described).

However, in the meantime, the possibility to access and somehow restore the history of relevant functions like delay could be an opportunity to develop specific libraries (it is not a general solution) that implement some particular functionalities. For example, the execution of parametric simulations or even Monte Carlo simulations can be implemented without the need for external code by developing a library of components prepared to save and restore their state under the control of an orchestrator within the same model. In this example, the lack of access to delay history prevents the inclusion of delay-based components.

I don't understand the link between parametric simulations/Monte Carlo simulations and delay-based components, can you clarify?

If you mean a form of Monte Carlo simulations with stochastic differential equations I would say that the lack of storing delay-based components is the least of the problems.

[mpradovelasco]

Of course @HansOlsson. I am finishing a functionality for a library that mainly deals with physiological - pharmacokinetic systems (it is called Cybim), where the delay function is important to manage several metrics (and other stories). For example, the area under the curve (AUC) along a time interval is a key measure of the activity of a molecule in a particular pharmacotherapy, and AUC is based on the delay().

The execution of both parametric and Monte Carlo analysis is a relevant scenario where this type of models are applied (Monte Carlo simulation is commonly used in population models that utilize a hierarchically based structure well known in pharmacometrics). Although both types of experiments can be defined in an external framework (e.g. using Python, Matlab), it is possible to implement these experiments through a kind of metamodel in OpenModelica, where one part of the model acts as an orchestrator and requests the reinit of the rest of the components. However, delay() based components need a method to control their history, otherwise these components cannot be reinit.

I know that there are additional details that need to be solved in this application, but these can be addressed in the context of a particular library. In fact, I have a solution that is nearing completion for CybSim. I hope to release it as open source this year.

Thanks for your interest.