Getting Started

The including integrators used by constructing a relaxation problem. Then by adding attributes to the relaxation problem and building the integrator from the relaxation problem. The integrator may then be modified by using the base API. Lastly, the relax!/integrate! are called to compute information about the trajectories and this information is retrieved using API functions.

Setting up a relaxation problem

We begin by constructing a relaxation problem. Currently, the ODERelaxProb is used for all parametric ODEs but specialized variants of this will become available in the future. This manner of relaxation problem is specified using the right-hand side function f! (evaluated in place) along with a function of p defining the initial conditions. A time span and box constraints on p are also expected.

using DynamicBounds

function f!(dx, x, p, t)
    dx[1] = -p[1]*x[1]*x[2] + p[2]*x[3] + p[6]*x[6]
    dx[2] = -p[1]*x[1]*x[2] + p[2]*x[3] + p[3]*x[3]
    dx[3] =  p[1]*x[1]*x[2] - p[2]*x[3] - p[3]*x[3]
    dx[4] =  p[3]*x[3] - p[4]*x[4]*x[5] + p[5]*x[6]
    dx[5] = -p[4]*x[4]*x[5] + p[5]*x[6] + p[6]*x[6]
    dx[6] =  p[4]*x[4]*x[5] - p[5]*x[6] - p[6]*x[6]
    return
end
x0(p) = [34.0, 20.0, 0.0, 0.0, 16.0, 0.0]

tspan = (0.0, 18.0e-5*250)
pL = [0.1; 0.033; 16.0; 5.0; 0.5; 0.3]
pU = 10.0*pL

ode_prob = ODERelaxProb(f!, x0, tspan, pL, pU)

We can often further refine bounds/relaxation on the state variables by making use of invariants and apriori bounds/relaxations. We'll add a polyhedral invariant and apriori bounds state bounds.

M = [0.0 -1.0 -1.0  0.0  0.0  0.0;
     0.0  0.0  0.0  0.0 -1.0 -1.0;
     1.0 -1.0  0.0  1.0 -1.0  0.0]
b = [-20.0; -16.0; -2.0]
polyhedral_invariant = PolyhedralConstraint(M, :==, b)
set!(ode_prob, polyhedral_invariant)

u_lo = zeros(6)
u_hi = [34.0; 20.0; 20.0; 34.0; 16.0; 16.0]
set!(prob, ConstantBounds(u_lo, u_hi))

We now can specify the points at which we want to retrieve information about state variables by adding a QueriedSupport attribute to the problem

set!(prob, QueriedSupport(range(tspan[1], tspan[2], length = 50)))

Creating an integrator from the problem

The integrator is built from the relaxation problem along with keyword variables.

integrator = DifferentialInequality(ode_prob, calculate_relax = true,
                                              calculate_subgradient = false)

Provided that we're interested in computing relaxations as well as constant bounds (w.r.t p) we'll need to specify parameter values within the box constraints.

setall!(integrator, ParameterValue(), p)

Relaxation problems also support the ParameterValue attributes. However, it's often preferable to set the value in the integrator to avoid rebuilding the integrator each time.

Computing relaxations and trajectories

We now compute the relevant information about the state variables and the solution of the pODEs.

relax!(integrator)      # compute relaxations
integrate!(integrator)  # compute trajectories

Querying the integrator/problem for information

We can now query functions using the get API function.

val = get(integrator, Value(3))
lo = get(integrator, Bound{Lower}(3))
hi = get(integrator, Bound{Upper}(3))
cv = get(integrator, Relaxation{Lower}(3))
cc = get(integrator, Relaxation{Upper}(3))

Use the sample problem library

DE_keys = keys(LibRelaxODE.STANDARD_LIBRARY)
for i = 1:length(keys)
    prob = LibRelaxODE.STANDARD_LIBRARY(DE_keys[i])
    integrator = DifferentialInequality(prob)
    relax!(integrator)
end