Process simulation solver

One of the most common tasks is to simply solve a process model for its independent variables at constant thermodynamic and process parameters. The model is solved, if all residuals evaluate to absolute values lower their tolerance.

The process model has the form

\[r(x, p, c) = 0\]

Here, \(x\) is the complete set of independent variables, \(p\) and \(c\) the process model and thermodynamic parameters respectively, and \(r\) the set of residuals. For the model to be solvable, the system must be square (\(\mathrm{dim}\ r = \mathrm{dim}\ x\)), and the system matrix \(J_{rx}\) must not be singular.

Iteratively, we can determine a raw Newton-Raphson update via

\[\Delta x^0\ \leftarrow\ J_{rx}\cdot \Delta x^0 = -r\]

Additionally, the model domain \(\mathbb D \ni x\) is limited, and its bounds described by

\[b(x, p, c) > 0\]

These are not to be confused with inequality constraints. In particular, the case \(b = 0\) is generally not feasible, and, as we solve an equation system, defining such bounds as inequality constraints would make the system non-square and thus not solvable.

Instead, we use \(b\) to limit the step size of the current iteration step, assuming and hoping that the solver will eventually find a solution within the domain. We define

\[A = \left \{-\frac{b_i}{\Delta b_i} \bigg | \Delta b_i < 0 \right \}\quad\text{with}\quad \Delta b = J_{bx}\cdot \Delta x^0\]

Based on a margin parameter \(\gamma \approx 0.9\) that describes how much we allow to approach the bounds, we determine the relaxation factor \(\alpha\) as

\[\alpha = \min(\{1, \gamma\cdot \min(A)\})\]

The applied variable update is then

\[\Delta x = x^\mathrm{(k+1)} - x^\mathrm{(k)} = \alpha\cdot \Delta x^0\]

Example: A pressure is 2 bar, and the raw update vector suggests a step of -3 bar. We receive an element in \(A\) of value 2/3. If this is the most limiting variable, the allowed step-size considering \(\gamma\) is 0.6. As this value is less than unity, it will be applied, and the pressure will be updated from 2 to 0.2 bar - we allow to approach the margin by 90 % of the possible step size.

Note

If we would allow the pressure to become truly zero, we would face two issues: Firstly, the model might divide by pressure or utilise the logarithm of it, and zero is already outside the domain. Secondly, numerical noise could still create negative values and thus consequences, even if zero still was in the domain.

With above scheme, the model is solved iteratively until the residual values drop below their tolerance.

SimulationSolver

class simu.SimulationSolver

Bases: Configurable

The simulation solver assumes both thermodynamic and model parameters to be constant, aiming to find the state variable values such that all residuals evaluate to zero within their tolerance.

__init__(
model: NumericHandler,
*,
max_iter: int = 30,
gamma: float = 0.9,
wall: float = 1e-20,
output: TextIOBase | str = 'stdout',
call_back_iter: Callable[[int, SimulationSolverIterationReport, Sequence[float], Callable[[Sequence[float]], NestedMap[Quantity]]], bool] = None,
retain_solution: bool = True,
)

On construction, the solver object requires a NumericHandler object. The solver object can then be reused for multiple solver runs, for instance with variable parameter values (sensitivity study).

Parameters:

  • model – The numeric handler of a model, in most cases obtained by the expression NumericHandler(ModelClass.top()).

  • max_iter –

    The maximum number of iterations (default 30).

    Note

    Normally, 30 iterations should be sufficient. In other words, if the model is not converged after 30 iterations, chances are quite low that it still will converge at all. The advice would be to try to improve the starting values and to investigate whether the model is properly posed.

  • gamma – As described above, \(\gamma\) (default 0.9) is the fraction of the step-length applied by the solver before hitting the domain boundary. Normally, changing the value is not required. Generally, a lower value makes the model more robust against non-linear domain boundaries (and thus linearisation errors causing the state to exit the domain). A higher value yields slightly faster convergence, if the solution is in comparison with the initial values very close to the domain boundary.

  • wall – Either if there is no solution within the domain of the model (for instance: The material balance forces some of the species flows in a stream to be negative), or if the solver for other reasons is forced to try to leave the model domain, the state will move closer and closer to the domain boundary and not revert. At some point, \(\gamma\) becomes ridiculously small, and we need to give up. This threshold value is defined by wall (default 1e-20).

  • output –

    The io-stream to direct the solver output to, or a descriptive string (case-insensitive):

    • "stdout": The output will be written to standard out (default)

    • "none": No output will be printed.

    Note

    Instead of printing, one might either analyse the returned SimulationSolverReport project after the run, or utilise the call_back_iter callback and process the iteration progress from there.

  • call_back_iter – A callback function (default None), see SimulationSolverCallback, to intercept the solving process. The returned boolean variable determines whether the solver iteration is continued or not.

  • retain_solution – Whether the solver shall, on success, retain the obtained state in the model, such that it can be exported via export_state() and be reused as the initial values for the next solving process.

solve(
**kwargs: Any,
) → SimulationSolverReport

This method triggers iterative the solving process. This takes less than 0.1 seconds for small models, and increases with model complexity. For large models, the time per iteration is due to the solving of linear systems cubic in system size, though the model structure might render this a conservative estimate.

With pypardiso installed, the solving of the linear systems is performed on all available CPU cores. However, their solver sometimes chokes and returns a wrong solution. Therefore, the norm of the solution is checked, and scipy.sparse.linalg.spsolve is used in those instances.

For each given keyword argument, the set_option() method is invoked, giving the same effect as if the option was provided with the constructor.

Returns:

The report including the iteration sequence

property initial_state

Freshly extract the initial values from the model. These might have been changed after the solver class was instantiated

property model_parameters: NestedMutMap[Quantity]

A convenience property to access the parameters of the model as a mutable object. The state variables are removed in this instance, as these are rather provided by the solver during the iterative solving process.

SimulationSolverIterationReport

class simu.core.solver.simulation.SimulationSolverIterationReport

This data class object is provided for each iteration during a SimulationSolver run.

max_err: float

For each Residual, the quotient of residual value \(r_i\) and tolerance \(t_i\) is calculated. max_err is the maximum absolute value of these quotients:

\[\mathrm{MET} = \max_i \frac{r_i}{t_i}\]
max_res_name: str

The name of the Residual which causes the value of max_err

relax_factor: float

The applied relaxation factor to stay within the domain of the process model and the thermodynamic models, according to the defined bounds.

min_alpha_name: str

The name of the bound that is most limiting and therefore causing the value of relax_factor

duration: float

The accumulative duration of the solving process inclusive the given iteration

lmet: float

The logarithmic (base 10) value of max_err \(r_i / t_i\), practically defined as

\[\mathrm{LMET} = \log_{10} \left ( \max_i \frac{r_i}{t_i} + 10^{-8} \right )\]

The offset is introduced to not cause NaN values for the lucky case in which all residuals are exactly zero. This can however easily happen for linear systems. Anyhow, \(\mathrm{LMET} < 1\) is already a sufficient condition for convergence.

condition: float

This is a fair condition of the system matrix, calculated after repeatedly scaling the rows and columns to eliminate the effect of variable scaling on the norm. As such, the obtained condition number is more suitable to judge the solvability of the model at hand.

The reported number is \(\log_{10}||\mathbf{J}||\).

SimulationSolverReport

class simu.core.solver.simulation.SimulationSolverReport

The data class object returned from a SimulationSolver run

iterations: Sequence[SimulationSolverIterationReport]

A :class:’SimulationSolverIterationReport` object for each performed iteration

final_state: Sequence[float]

The numerical final state of the model.

Important

This state is not suitable for handling initial values in a robust way, as it can be very sensitive to minor model changes that for instance impact the liquid volumes of equations of state.

Instead, use simu.NumericHandler.export_state(), import_state() and retain_state().

prop_func: Callable[[Sequence[float]], NestedMap[Quantity]]

The function to calculate all properties of the model as function of the state.

property properties: NestedMap[Quantity]

This property returns all properties of the model, evaluated on the final_state attribute. For larger models, this causes noticeable computational effort. For this reason, this evaluation is only done on demand.

SimulationSolverCallback

simu.core.solver.simulation.SimulationSolverCallback

A function (type) to act as a call-back in the SimulationSolver solving process. The arguments are as follows:

  • iteration: The iteration number as integer, incrementing from zero

  • report: The (SimulationSolverIterationReport) object

  • state: The internal state of the model at given iteration as a sequence of floats

  • prop_func: A function to calculate all model properties for the given state.

The last argument is provided instead of the property structure itself, as it might be expensive to calculate the entire property structure in each iteration. This way it can be done on demand.

The callback shall return True, if the solver is to continue, or False otherwise. In the latter case, the solver will raise a IterativeProcessInterrupted exception.

Example:

1 from pprint import pprint
2
3 def my_callback(iteration, report, state, prop_func):
4     # This can be a lot to print
5     all_properties = prop_func(state)
6     pprint(all_properties)
7     return True

alias of Callable[[int, SimulationSolverIterationReport, Sequence[float], Callable[[Sequence[float]], NestedMap[Quantity]]], bool]