2.4.14. Adjoint Sensitivity Analysis

Added in version 6.3.0.

The previous sections discuss using ARKODE for the integration of forward ODE models. This section discusses how to use ARKODE for adjoint sensitivity analysis as introduced in §2.2.19. To use ARKStep for adjoint sensitivity analysis (ASA), users simply setup the forward integration as usual (following §2.4.2) with a few differences. Below we provide an updated version of the ARKODE usage in section §2.4.2 where steps that are unchanged are italicized. The example code examples/arkode/C_serial/ark_lotka_volterra_asa.c demonstrates these steps in detail.

Initialize parallel or multi-threaded environment, if appropriate.
Create the SUNDIALS simulation context object
Set problem dimensions, etc.
Set vector of initial values
Create ARKODE object
Specify a fixed time step size.

Currently the discrete ASA capability only allows a fixed time step size to be used. Call ARKodeSetFixedStep() to set the time step.
Set optional inputs
Specify rootfinding problem
Create a SUNAdjointCheckpointScheme object

Create the SUNAdjointCheckpointScheme object by calling SUNAdjointCheckpointScheme_Create_*. Available SUNAdjointCheckpointScheme implementations are found in section §14.3.
Attach the checkpoint scheme object to ARKODE

Call ARKodeSetAdjointCheckpointScheme().
Advance solution in time
Get optional outputs
Create the sensitivities vector with the terminal condition

The sensitivities vector must be an instance of the ManyVector N_Vector implementation. You will have one subvector for the initial condition sensitivities and an additional subvector if you want sensitivities with respect to parameters. The vectors should contain the terminal conditions for the adjoint problem. The first subvector should contain \(dg(t_f,y(t_f),p)/dy(t_f)\) and the second subvector should contain \(dg(t_f,y(t_f),p)/dp\). The subvectors can be any implementation of the N_Vector class.

For example, in a problem with 10 state variables and 4 parameters using serial computations, the ManyVector can be constructed as follows:
```
sunindextype num_equations = 10;
sunindextype num_params    = 4;
N_Vector sensu0            = N_VNew_Serial(num_equations, sunctx);
N_Vector sensp             = N_VNew_Serial(num_params, sunctx);
N_Vector sens[2]           = {sensu0, sensp};
N_Vector sf                = N_VNew_ManyVector(2, sens, sunctx);
// Set the terminal condition for the adjoint system, which
// should be the the gradient of our cost function at tf.
dgdu(u, sensu0, params);
dgdp(u, sensp, params);
```
Create the SUNAdjointStepper object

Call ERKStepCreateAdjointStepper() or ARKStepCreateAdjointStepper().
Set optional ASA input

Refer to §14.2 for options.
Advance the adjoint sensitivity analysis ODE

Call SUNAdjointStepper_Evolve() or SUNAdjointStepper_OneStep().
Get optional ASA outputs

Refer to §14.2 for options.
Deallocate memory for ASA objects

Deallocate the sensitivities vector, SUNAdjointStepper, and SUNAdjointCheckpointScheme objects.
Deallocate memory for solution vector
Free solver memory

Call SUNStepper_Destroy() and ARKodeFree() to free the memory allocated for the SUNStepper and ARKODE integrator objects.
Free the SUNContext object
Finalize MPI, if used

2.4.14.1. User Callable Functions

This section describes user-callable functions for performing adjoint sensitivity analysis with methods with ERKStep and ARKStep.

int ERKStepCreateAdjointStepper(void *arkode_mem, SUNAdjRhsFn f, sunrealtype tf, N_Vector sf, SUNContext sunctx, SUNAdjointStepper *adj_stepper_ptr)

Creates a SUNAdjointStepper object compatible with the provided ERKStep instance for integrating the adjoint sensitivity system (2.72).

Parameters:

arkode_mem – a pointer to the ERKStep memory block.
f – the adjoint right hand side function which implements \(\Lambda = f_y^*(t, y, p) \lambda\) and, if sensitivities with respect to parameters should be computed, \(\nu = f_p^*(t, y, p) \lambda\).
tf – the terminal time for the adjoint sensitivity system.
sf – the sensitivity vector holding the adjoint system terminal condition. This must be an NVECTOR_MANYVECTOR instance. The first subvector must be \(g^*_y(t_f, y(t_f), p) \in \mathbb{R}^N\). If sensitivities to parameters should be computed, then the second subvector must be \(g^*_p(t_f, y(t_f), p) \in \mathbb{R}^{N_s}\), otherwise only one subvector should be provided.
sunctx – the SUNDIALS simulation context object.
adj_stepper_ptr – the newly created SUNAdjointStepper object.

Return values:

ARK_SUCCESS – if successful.
ARK_MEM_FAIL – if a memory allocation failed.
ARK_ILL_INPUT – if an argument has an illegal value.

Added in version 6.3.0.

Note

Currently fixed time steps must be used. Furthermore, the explicit stability function, inequality constraints, and relaxation features are not yet compatible as they require adaptive time steps.

int ARKStepCreateAdjointStepper(void *arkode_mem, SUNAdjRhsFn fe, SUNAdjRhsFn fi, sunrealtype tf, N_Vector sf, SUNContext sunctx, SUNAdjointStepper *adj_stepper_ptr)

Creates a SUNAdjointStepper object compatible with the provided ARKStep instance for integrating the adjoint sensitivity system (2.72).

Parameters:

arkode_mem – a pointer to the ARKStep memory block.
fe – the adjoint right hand side function which implements \(\Lambda = f_y^{E,*}(t, y, p) \lambda\) and, if sensitivities with respect to parameters should be computed, \(\nu = f_p^*(t, y, p) \lambda\).
fi – not yet support, the user should pass NULL.
tf – the terminal time for the adjoint sensitivity system.
sf – the sensitivity vector holding the adjoint system terminal condition. This must be a NVECTOR_MANYVECTOR instance. The first subvector must be \(g^*_y(t_f, y(t_f), p) \in \mathbb{R}^N\). If sensitivities to parameters should be computed, then the second subvector must be \(g^*_p(t_f, y(t_f), p) \in \mathbb{R}^{N_s}\), otherwise only one subvector should be provided.
sunctx – the SUNDIALS simulation context object.
adj_stepper_ptr – the newly created SUNAdjointStepper object.

Return values:

ARK_SUCCESS – if successful.
ARK_MEM_FAIL – if a memory allocation failed.
ARK_ILL_INPUT – if an argument has an illegal value.

Added in version 6.3.0.

Note

Currently only explicit methods with identity mass matrices are supported for ASA, and fixed time steps must be used. Furthermore, the explicit stability function, inequality constraints, and relaxation features are not yet compatible as they require adaptive time steps.