Differential-Algebraic Equations
Stan support solving systems of differential-algebraic equations (DAEs) of index 1 (Serban et al. 2021). The solver adaptively refines the solutions in order to satisfy given tolerances.
One can think a differential-algebraic system of equations as ODEs with additional algebraic constraints applied to some of the variables. In such a system, the variable derivatives may not be expressed explicitly with a right-hand-side as in ODEs, but implicitly constrained.
Similar to ODE solvers, the DAE solvers must not only provide the solution to the DAE itself, but also the gradient of the DAE solution with respect to parameters (the sensitivities). Stan’s DAE solver uses the forward sensitivity calculation to expand the base DAE system with additional DAE equations for the gradients of the solution. For each parameter, an additional full set of
Two interfaces are provided for the forward sensitivity solver: one with default tolerances and default max number of steps, and one that allows these controls to be modified. Choosing tolerances is important for making any of the solvers work well – the defaults will not work everywhere. The tolerances should be chosen primarily with consideration to the scales of the solutions, the accuracy needed for the solutions, and how the solutions are used in the model. The same principles in the control parameters section apply here.
Internally Stan’s DAE solver uses a variable-step, variable-order, backward-differentiation formula implementation (Cohen and Hindmarsh 1996; Serban and Hindmarsh 2005).
Notation
A DAE is defined by a set of expressions for the residuals of differential equations and algebraic equations
Example: chemical kinetics
As an example of a system of DAEs, consider following chemical kinetics problem(Robertson 1966). The nondimensionalized DAE consists of two differential equations and one algebraic constraint. The differential equations describe the reactions from reactants
The state equations implicitly defines the state
Unlike solving ODEs, solving DAEs requires a consistent initial condition. That is, one must specify both
Index of DAEs
The index along a DAE solution
Most DAE solvers, including the one in Stan, support only index-1 DAEs. For a high index DAE problem the user must first convert it to a lower index system. This often can be done by carrying out differentiations analytically (Ascher and Petzold 1998).
Coding the DAE system function
The first step in coding an DAE system in Stan is defining the DAE residual function. The system functions require a specific signature so that the solvers know how to use them properly.
The first argument to the residual function is time, passed as a real
; the second argument to the residual function is the system state vector
, the third argument to the residual function is the state derivative vector
. The residual function’s return value is a vector
of the same size as state and state derivatives. Additional arguments can be included in the residual function to pass other information into the solve (these will be passed through the function that starts the DAE solution). These argument can be parameters (in our example, the reaction rate coefficient
The above reaction be coded using the following function in Stan (see the user-defined functions chapter for more information on coding user-defined functions).
vector chem(real t, vector yy, vector yp,
real alpha, real beta, real gamma) {
vector[3] res;
1] = yp[1] + alpha * yy[1] - beta * yy[2] * yy[3];
res[2] = yp[2] - alpha * yy[1] + beta * yy[2] * yy[3] + gamma * yy[2] * yy[2];
res[3] = yy[1] + yy[2] + yy[3] - 1.0;
res[return res;
} }
The function takes in a time t
(a real
), the system state yy
(a vector
), state derivative yp
(a vector
), as well as parameter alpha
(a real
), beta
(a real
), and gamma
(a real
). The function returns a vector
of the residuals at time t
. The DAE coded here does not explicitly depend on t
, however one still needs to specify t
as an argument.
Strict signature
The types in the DAE residual function are strict. The first argument is the time passed as a real
, the second argument is the state passed as a vector
, the third argument is the state derivative passed as a vector
, and the return type is a vector
. A model that does not have this signature will fail to compile. The fourth argument onwards can be any type, granted all the argument types match the types of the respective arguments in the solver call.
All of these are possible DAE signatures:
vector my_dae1(real t, vector y, vector yp, real a0);
vector my_dae2(real t, vector y, vector yp, array[] int a0, vector a1);
vector my_dae3(real t, vector y, vector yp, matrix a0, array[] real a1, row_vector a2);
but these are not allowed:
vector my_dae1(real t, array[] real y, vector yp);
// Second argument is not a vector
array[] real my_dae2(real t, vector y, vector yp);
// Return type is not a vector
vector my_dae3(real t, vector y);
// First argument is not a real and missing the third argument
Solving DAEs
Stan provides a dae
function for solving DAEs, so that the above chemical reaction equation can be solved in the following code.
data {
int N;
vector[3] yy0;
vector[3] yp0;
real t0;
real alpha;
real beta;
array[N] real ts;
array[N] vector[3] y;
}parameters {
real gamma;
}transformed parameters {
vector[3] y_hat[N] = dae(chem, yy0, yp0, t0, ts, alpha, beta, gamma);
}
Since gamma
is a parameter, the DAE solver is called in the transformed parameters block.
Control parameters for DAE solving
Using dae_tol
one can specify the relative_tolerance
, absolute_tolerance
, and max_num_steps
parameters in order to control the DAE solution.
vector[3] y_hat[N] = dae_tol(chem, yy0, yp0, t0, ts,
relative_tolerance,
absolute_tolerance,
max_num_steps, alpha, beta, gamma);
relative_tolerance
and absolute_tolerance
control accuracy the solver tries to achieve, and max_num_steps
specifies the maximum number of steps the solver will take between output time points before throwing an error.
The control parameters must be data variables – they cannot be parameters or expressions that depend on parameters, including local variables in any block other than transformed data and generated quantities. User-defined function arguments may be qualified as only allowing data arguments using the data
qualifier.
The default value of relative and absolute tolerances are
Maximum number of steps
The maximum number of steps can be used to stop a runaway simulation. This can arise in when MCMC moves to a part of parameter space very far from where a differential equation would typically be solved. In particular this can happen during warmup. With the non-stiff solver, this may happen when the sampler moves to stiff regions of parameter space, which will requires small step sizes.