DynIbex is a plug-in of Ibex library which is a library for constraint processing over real numbers.

DynIbex offers a set of validated numerical integration methods based on Runge-Kutta schemes to solve initial value problem of ordinary differential equations and for DAE in Hessenberg index 1 form.

• We consider an ordinary differential equation with a given initial condition, that is an initial value problem (IVP) $\dot{\mathbf{y}}(t) = f(t, \mathbf{y}(t)) \quad \text{with}\quad \mathbf{y}(0)= \mathbf{y}_0$ with $$f : \mathbb{R} \times \mathbb{R}^n \to \mathbb{R}^n$$ is assumed to be continuous $$t$$ and globally Lipschitz in $$\mathbf{y}$$. Solving an IVP consists in finding a function $$\mathbf{y}(t)$$ described by the ODE and satisfying the initial condition. The solution of an IVP is denoted by $$\mathbf{y}(t;\mathbf{y}_0)$$. A standard numerical integration method computes a sequence of values $$(t_n, \mathbf{y}_n)_{n\in\mathbb{N}}$$ approximating the solution $$\mathbf{y}(t;\mathbf{y}_0)$$ of an IVP such that $$\mathbf{y}_n \approx \mathbf{y}(t_n;\mathbf{y}_{n-1})$$. The goal of a validated or guaranteed numerical integration method is to compute the sequence of boxes $$(t_n, [\mathbf{y}_n])_{n\in\mathbb{N}}$$ such that $$[\mathbf{y}_n] \supseteq \mathbf{y}(t_n;[\mathbf{y}_{n-1}]) = \{ \mathbf{y}(t_n;\mathbf{y}_{n-1}) : \forall \mathbf{y}_{n-1} \in [\mathbf{y}_{n-1}] \}$$.
• We also consider a differential-algebraic equations in Hessenberg index 1 form with consistent initial conditions, that IVP of the form \begin{align*} \dot{\mathbf{y}}(t) & = f (t, \mathbf{y}(t), \mathbf{x}(t)) \\ 0 & = g(t, \mathbf{y}(t), \mathbf{x}(t)) \end{align*} \quad\text{with}\quad \mathbf{y}(0) = \mathbf{y}_0\quad\text{and}\quad \mathbf{x}(0) = \mathbf{x}_0 with $$f : \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n$$ is assumed to be continuous $$t$$ and globally Lipschitz in $$\mathbf{y}$$, and $$g : \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$$. Solving an IVP consists in finding two functions $$\mathbf{y}(t)$$ and $$\mathbf{x}(t)$$ described by the DAE and satisfying the initial conditions. The solutions of an IVP for DAE is denoted by $$\mathbf{y}(t;\mathbf{y}_0)$$ and $$\mathbf{x}(t;\mathbf{x}_0)$$. Once again, the goal of a validated or guaranteed numerical integration method is to compute two sequences of boxes $$(t_n, [\mathbf{y}_n])_{n\in\mathbb{N}}$$ and $$(t_n, [\mathbf{x}_n])_{n\in\mathbb{N}}$$ such that $$[\mathbf{y}_n] \supseteq \mathbf{y}(t_n;[\mathbf{y}_{n-1}]) = \{ \mathbf{y}(t_n;\mathbf{y}_{n-1}) : \forall \mathbf{x}_{n-1} \in [\mathbf{x}_{n-1}] \}$$ and $$[\mathbf{x}_n] \supseteq \mathbf{x}(t_n;[\mathbf{x}_{n-1}]) = \{ \mathbf{x}(t_n;\mathbf{x}_{n-1}) : \forall \mathbf{x}_{n-1} \in [\mathbf{x}_{n-1}] \}$$

## Contributors

Julien Alexandre dit Sandretto, Alexandre Chapoutot, Olivier Mullier

This research benefited from the support of the "Chair Complex Systems Engineering - Ecole Polytechnique, THALES, DGA, FX, DASSAULT AVIATION, DCNS Research, ENSTA ParisTech, Télécom ParisTech, Fondation ParisTech and FDO ENSTA"

## Features and Limitations

• We can solve IVP for ODE and DAE in Hessenberg index 1 form
• Implement different Runge-Kutta methods: Heun's method, Midpoint method, Radau-IIA (order 3), RK4, Lobatto-IIIA (order 4), Lobatto-IIIC (order 4)
• With affine arithmetic, uncertain ODE can be considered without putting parameters in the state-vector
• With affine arithmetic, QR-method is not used
• The current algorithm to compute LTE has exponential complexity so the dimension of ODE to solve is limited (< 5)

## Ongoing and future work (To be aware of the latest versions of DynIbex, subscribe RSS feed)

• Reduce the complexity of the computation of LTE (ongoing, testing/debuging phase)
• Extension of DynIbex to deal with hybrid systems as described in the report (future)
• Embed DynIbex in contractor programming framework (future)

## History on Interval Runge-Kutta

This work on validated Runge-Kutta methods takes its roots in two previous works:
1. GRKlib by Olivier Bouissou and Matthieu Martel which defined a validated version of the classical RK4 methods
2. The extension to deal with any explicit Runge-Kutta methods as defined in the article published at NASA Formal Methods 2013 written by Olivier Bouissou, Alexandre Chapoutot and Adel Djoudi.
Note that other tentative to make Runge-Kutta methods guaranteed have been done, mainly:
• Andrzej Marciniak et alii work on this topic since 1999
We go a step forward with our approach by bringing an algorithmic answer to the remark found in the survey article saying that the LTE cannot be written in general form for an arbitrary order.
• Michael Hartmann and Knut Petras at ICIAM 1999 presented some work on this topic but we cannot find no more than 5 lines in an abstract (more information is welcomed).

## List of publications on DynIBEX

DynIbex is made of two main ingredients: an extension of the affine arithmetic and the validated Runge-Kutta methods.

We assume that the following sequence of commands are executed in Ibex directory.

Put the archive into the Ibex directory, then write the following commands:

      Ibex foo$tar xvfz tmp_install.tar.gz -C ibex-2.1.x Ibex foo$ cd ibex-2.1.x
Ibex/ibex-2.1.x foo$python ./install_plugin.py Ibex/ibex-2.1.x foo$ ./waf configure [options]
Ibex/ibex-2.1.x foo\$ ./waf install


# Affine arithmetic

## Description

Affine arithmetic is an extension of interval arithmetic which can track linear correlation between variables. It is suitable to reduce the pessimism due to the dependency problem and the wrapping effect. These properties are very important to implement validated numerical integration methods.

## Example

A simple example showing the ability to reduce dependency problem using affine arithmetic. Note a more complete example on the use of affine arithmetic is given in the distribution of DynIbex in the examples repository, file test_affine_full.cpp.

#include <iostream>;
#include "ibex.h"

using namespace std;
using namespace ibex;

int main(){

Interval x(0.0, 10.3);
Interval y(-1.0, 1.2);
Affine2 ax(x);
Affine2 ay(y);
Interval res = (x + y) - x - y;
Affine2 resa = (ax + ay) - ax - ay;

cout << "Interval result = " << res << endl;
cout << "Affine result = " << resa << endl;

return 0;
}



# Validated Runge-Kutta methods for ODE (Research report and Preprint)

## Description

A Runge-Kutta method, starting from an initial value $$\mathbf{y}_n$$ at time $$t_n$$ and a finite time horizon $$h$$, the step-size, produces an approximation $$\mathbf{y}_{n+1}$$ at time $$t_{n+1}$$, with $$t_{n+1}-t_n = h$$, of the solution $$\mathbf{y}(t_{n+1}; \mathbf{y}_n)$$. Furthermore, to compute $$\mathbf{y}_{n+1}$$, a Runge-Kutta method computes $$s$$ evaluations of $$f$$ at predetermined time instants. The number $$s$$ is known as the number of stages of a Runge-Kutta method. More precisely, a Runge-Kutta method is defined by $\mathbf{y}_{n+1} = \mathbf{y}_n + h \sum_{i=1}^{s}b_i \mathbf{k}_i\enspace,$ with $$\mathbf{k}_i$$ defined by $\mathbf{k}_i = f\left(t_0 + c_i h, \mathbf{y}_n + h \sum_{j=1}^{s} a_{ij}\mathbf{k}_j \right)\enspace.$ The coefficient $$c_i$$, $$a_{ij}$$ and $$b_i$$, for $$i,j=1,2,\cdots, s$$, fully characterize the Runge-Kutta methods and their are usually synthesized in a Butcher tableau of the form as follows: $\begin{array}{c|cccc} c_1 & a_{11} & a_{12} & \dots & a_{1s}\\ c_2 & a_{21} & a_{22} & \dots & a_{2s}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ c_s & a_{s1} & a_{s2} & \dots & a_{ss} \\ \hline & b_1 & b_2 & \dots & b_s\\ \end{array}$ In function of the form of the matrix $$A$$, made of the coefficients $$a_{ij}$$, a Runge-Kutta method can be
• explicit the computation of an intermediate $$k_i$$ only depends on the previous steps $$k_j$$ for $$j < i$$;
• diagonally implicit the computation of an intermediate step $$k_i$$ involves the value $$k_i$$ and so non-linear systems in $$k_i$$ must be solved;
• fully implicit the computation of intermediate steps involves the solution of a non-linear system of equations in all the values $$k_i$$ for $$i=1,2,\cdots,s$$.
To make a Runge-Kutta validated, the local truncation error (LTE) need to be bound, that is $\parallel \mathbf{y}\left(t_n; \mathbf{y}_{n-1}\right) - \mathbf{y}_n \parallel \leqslant c h^{p+1}$ $$c$$ is a constant depending on the coefficients of the Runge-Kutta methods and $$p$$ is the order of the Runge-Kutta method. An analytic form can be defined using the theorems stating the order condition of a Runge-Kutta method see our report for more details.

## Example: basic integration

This example shows the simple way to solve IVP with DynIbex.

#include "ibex.h"
using namespace ibex;

#define __PREC__ 1e-11
#define __METH__ RK4
#define __DURATION__ 10.0

int main(){
// Dimension of the IVP
const int n= 2;

// State variables
Variable y(n);

// Initial conditions
IntervalVector yinit(n);
yinit[0] = Interval(1.5,2.0);
yinit[1] = Interval(-2.0,-1.5);

// Derivative function
Function ydot = Function (y,Return (y[0]-2*y[1],
3*y[0]-4*y[1]));

// Ivp contruction (initial time is 0.0)
ivp_ode problem = ivp_ode (ydot, 0.0, yinit);

// Simulation construction and run
simulation simu = simulation (&problem, __DURATION__, __METH__, __PREC__);
simu.run_simulation();

return 0;
}



# Validated Runge-Kutta methods for DAE (Preprint)

## Description

The validated numerical integration of DAE in Hessenberg index 1 form can be also based on Runge-Kutta scheme. Nevertheless, it implies the definition of a new operator for the computation of the a priori enclosure. For more details, see the article.

## Example: a simple DAE


#include "ibex.h"

using namespace ibex;

int main(){
// State variables
Variable y(3);

// Initialisation of y
IntervalVector yinit(3);
yinit[0]= Interval(1.0);
yinit[1]= 1.0;
yinit[2]= 3.0;

// Algebraic variables
Variable x(2);

// Initialisation of x
IntervalVector xinit(2);
xinit[0]= Interval(0.5);
xinit[1]= 1.0;

// Derivative function
Function ydot = Function(y, x, Return(y[1]+x[0],
y[0]-y[1]*x[0],
y[0]*y[2]-x[1]));

// Explicit constraints
Function g = Function(y, x, Return(y[0]-x[1], y[1]-2*x[0]));

// IVP construction (initial time is 0.0)
ivp_dae_h1 problem = ivp_dae_h1(ydot, g, 0.0, yinit, xinit);

// Simulation construction and run (final time is 0.5)
simulation simu = simulation(&problem, 0.5, RADAU3_DAE, 1e-14);
simu.run_simulation();

return 0;
}



# Towards a contractor programming with differential equations

## Description

This example has been written with the help of Gilles Chabert. It demonstrates the ability to embed validated Runge-Kutta methods for ODE into the contractor programming framework. Nevertheless, more work has to be done to have a complete and robust framework of constraint programming techniques with differential equations.

## Example: a simple forward-backward operator for ODEs.


#include "ibex.h"

using namespace ibex;

// Technical stuff to automatically compute the negation of an Ibex function
Function operator-(Function& f) {
Array <const ExprSymbol> args(f.nb_arg());
varcopy (f.args(), args);
const ExprConstant* cst = dynamic_cast<const ExprConstant*> (&f.expr());
if (cst) {
if (cst -> dim.is_vector ()) return Function (args, -cst->get_vector_value());
else return Function(args, -cst->get_matrix_value());
} else {
const ExprVector* vec = dynamic_cast<const ExprVector*>(&f.expr());
assert(vec); // TODO: manage matrix-valued function
Array<const ExprNode> minus_vec(vec -> nb_args);
for (int i = 0; i < vec -> nb_args; i++) {
minus_vec.set_ref (i, -ExprCopy().copy(f.args(), args, vec->arg(i)));
}

return Function (args, ExprVector::new_(minus_vec,vec->row_vector()));
}
}

// IVP Forward contractor
class CtcIVPFwd : public Ctc {
public:
CtcIVPFwd(const Function& f, double t0, double t1) : Ctc(f.nb_var()), f(f), t0(t0), x0(f.nb_var()), t1(t1) {
}

void set_x0(const IntervalVector& x0) {
this->x0 = x0;
}

void contract(IntervalVector& x1) {
int n = f.nb_var();
ivp_ode ode = ivp_ode(f, t0, x0);
simulation simu = simulation(&ode, t1, RK4, 1e-9);
simu.run_simulation();
x1 &= simu.get_last();
}

Function f; // TODO: avoid copy?
IntervalVector x0;
double t0;
double t1;

};

// IVP Forward-Backward contractor
// Backward contractor is an IVP Forward contractor with dynamics -f
class CtcIVP : public Ctc {
public:
CtcIVP(Function& f, double t0, double t1) : Ctc(f.nb_var()*2), fwd(f,t0,t1), bwd(-f,t0,t1) {

}

void contract(IntervalVector& x0x1) {
int n = x0x1.size() / 2;
IntervalVector x0 = x0x1.subvector(0,n-1);
IntervalVector x1(n);
IntervalVector x1_save = x0x1.subvector(n,2*n-1);

fwd.set_x0(x0);
fwd.contract(x1);

if (!x1_save.is_superset(x1)) {
x1 &= x1_save;
if (!x1.is_empty()) {
bwd.set_x0(x1);
bwd.contract(x0);
}
}

if (x0.is_empty() || x1.is_empty()) {
x0x1.set_empty();
throw EmptyBoxException();
} else {
x0x1.put(0,x0);
x0x1.put(n,x1);
}
}

CtcIVPFwd fwd;
CtcIVPFwd bwd;
};

int main(){
// IVP dimension
const int n = 2;

// Some uncertainty
Interval eps = 1000.0 * Interval(-1,1);

// Initial condition at time 0.0
IntervalVector yinit(n);
yinit[0] = 10.0 + eps;
yinit[1] = 0.0 + eps;

// Final condition a time 2.0
IntervalVector yfinal(n);
yfinal[0] = -9.62 + eps;
yfinal[1] = -19.62;

// Numerical contractor
NumConstraint dist ("x0[2]","x1[2]","x0(1)+x1(1)=0");
CtcFwdBwd c2 (dist);

// IVP contractor
Variable y(n);
Function ydot = Function(y,Return(y[1], -Interval(9.81)));
CtcIVP c1 (ydot,0.0,2.0);
IntervalVector box(4);
box.put(0,yinit);
box.put(n,yfinal);

// Classical composition of contractors
CtcCompo compo(c1, c2);
CtcFixPoint fix(compo);

// Bisection strategy
LargestFirst lf(0.0001);
CellStack _stack;

// Default solver
Solver solver(fix, lf, _stack);

vector<IntervalVector> sols = solver.solve(box);

if (sols.empty()) {
cout << "no solution" << endl;
}
for (int i = 0; i < sols.size(); i++) {
cout << "sols n°" << i << "= " << sols[i] << endl;
}

return 0;
}



# Simulation of a car in open loop, continuous-time closed loop and sampled closed loop

## Description

We consider a simple nonlinear dynamics of a car given by equations \begin{aligned} \dot{v}(t) & = \frac{ -50.0 v - 0.4 v^2}{m} \\ & \text{with}\quad m \in [990,1010] \end{aligned}

## Example: open loop system


#include "ibex.h"

using namespace ibex;

int main(){

const int n = 1;
Variable y(n);

IntervalVector state(n);
state[0] = 10.0;

// Dynamique d'une voiture avec incertitude sur sa masse
Function ydot(y, ( -50.0 * y[0] - 0.4 * y[0] * y[0])
/ Interval (990, 1010));
ivp_ode vdp = ivp_ode(ydot, 0.0, state);

// Integration numerique ensembliste
simulation simu = simulation(&vdp, 100, RK4, 1e-6);
simu.run_simulation();

//For an export in order to plot
simu.export1d_yn("export-open-loop.txt", 0);

return 0;
}



## Example: continuous closed loop system with PI controller

• We increase the dimension of the dynamical system to embed the integral part of the PI controller $\begin{pmatrix} \dot{v} \\ \dot{w} \end{pmatrix} = \begin{pmatrix} \frac{k_p (10.0 - v) + k_i w -50.0v - 0.4v^2}{m} \\ 10.0 - v \end{pmatrix} \quad\text{with}\quad m\in[990, 1010], k_p=1440, k_i=35$ The control input is $$10$$

#include "ibex.h"

using namespace ibex;

int main(){

const int n = 2;
Variable y(n);

IntervalVector state(n);
state[0] = 0.0;
state[1] = 0.0;

// Dynamique d'une voiture avec incertitude sur sa masse + PI
Function ydot(y, Return ((1440.0 * (10.0 - y[0]) + 35.0 * y[1]  - y[0] * (50.0  + 0.4 * y[0]))
/ Interval (990, 1010),
10.0 - y[0]));
ivp_ode vdp = ivp_ode(ydot, 0.0, state);

// Integration numerique ensembliste
simulation simu = simulation(&vdp, 10.0, RK4, 1e-6);
simu.run_simulation();

simu.export1d_yn("export-closed-loop.txt", 0);

return 0;
}



## Example: continuous closed loop system with PI controller with safety property


#include "ibex.h"

using namespace ibex;
using namespace std;

int main(){

const int n = 2;
Variable y(n);

IntervalVector state(n);
state[0] = 0.0;
state[1] = 0.0;

// Dynamique d'une voiture avec incertitude sur sa masse + PI
Function ydot(y, Return ((1440.0 * (10.0 - y[0]) + 35.0 * y[1]  - y[0] * (50.0  + 0.4 * y[0]))
/ Interval (990, 1010),
10.0 - y[0]));
ivp_ode vdp = ivp_ode(ydot, 0.0, state);

simulation simu = simulation(&vdp, 10.0, RK4, 1e-6);
simu.run_simulation();

// verification de surete
IntervalVector safe(n);
safe[0] = Interval(0.0, 11.0);
bool flag = simu.stayed_in (safe);
if (!flag) {
std::cerr << "\t\t\033[31mERROR SAFETY VIOLATION\033[0m" << std::endl;
}

return 0;
}



## Example: closed loop system with discretized PI controller


#include "ibex.h"

using namespace ibex;
using namespace std;

int main(){

const int n = 2;
Variable y(n);

// State vector is a zonotope
Affine2Vector state(n);
state[0] = 0.0;
state[1] = 0.0;

double t = 0;
const double sampling = 0.005;
Affine2 integral(0.0);

while (t < 3.0) {
Affine2 goal(10.0);
Affine2 error = goal - state[0];

// Controleur PI discret
integral = integral + sampling * error;
Affine2 u = 1400.0 * error + 35.0 * integral;
state[1] = u;

// Dynamique d'une voiture avec incertitude sur sa masse
Function ydot(y, Return((y[1] - 50.0 * y[0] - 0.4 * y[0] * y[0])
/ Interval (990, 1010), Interval(0.0)));
ivp_ode vdp = ivp_ode(ydot, 0.0, state);

// Integration numerique ensembliste
simulation simu = simulation(&vdp, sampling, RK4, 1e-6);
simu.run_simulation();

// Mise a jour du temps et des etats
state = simu.get_last();
t += sampling;

cerr << state.itv() << endl;
}

return 0;
}