Menu:

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.

Contributors

Julien Alexandre dit Sandretto, Alexandre Chapoutot, Olivier Mullier

Sponsors

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

How to cite DynIBEX

You should cite our paper published in Reliable Computing Journal with the following bibtex entry

@Article{Sandretto:2016:VEI,
  Title                    = {Validated Explicit and Implicit {R}unge--{K}utta Methods},
  Author                   = {Julien Alexandre {dit Sandretto} and Alexandre Chapoutot},
  Journal                  = {Reliable Computing},
  Year                     = {2016},
  Month                    = {Jul},
  Number                   = {1},
  Pages                    = {79--103},
  Volume                   = {22}
}

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:

DOWNLOAD and INSTALLATION

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/$ tar xvfz dynibex-2.0.tar.gz
	    ~/Ibex/$ cd dynibex-2.0
	    ~/Ibex/dynibex-2.0/$ ./waf configure [options]
	    ~/Ibex/dynibex-2.0/$ sudo ./waf install
	  


Publications on DynIBEX

Theory of validated Runge-Kutta methods

Applications using DynIBEX


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

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 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 (Article)

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


#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;
}



Validated Runge-Kutta methods for dynamical systems with state constraints

Description

The validated numerical integration of dynamical systems with state constraints is available. It combines validated Runge-Kutta scheme with contractor operators to ensure the respect of constraints.

Example: a simple Production-Destruction systems


#include "ibex.h"
#include "vibes.cpp"


using namespace ibex;


void plot_simu(const simulation* sim)
{
  std::list<solution_g>::const_iterator iterator_list;
  for(iterator_list=sim->list_solution_g.begin();iterator_list!=sim->list_solution_g.end();iterator_list++) {
    vibes::drawBox(iterator_list->time_j.lb(), iterator_list->time_j.ub(),
		   iterator_list->box_j1->operator[](0).lb(), iterator_list->box_j1->operator[](0).ub(), "red[red]");

    vibes::drawBox(iterator_list->time_j.lb(), iterator_list->time_j.ub(),
		   iterator_list->box_j1->operator[](1).lb(), iterator_list->box_j1->operator[](1).ub(), "blue[blue]");

    vibes::drawBox(iterator_list->time_j.lb(), iterator_list->time_j.ub(),
		   iterator_list->box_j1->operator[](2).lb(), iterator_list->box_j1->operator[](2).ub(), "black[black]");
  }
}


int main(){

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

  IntervalVector yinit(n);
  yinit[0] = Interval(9.98);
  yinit[1] = Interval(0.01);
  yinit[2] = Interval(0.01);
  Interval a(0.3);

  Function ydot = Function(y,Return(-y[0]*y[1]/(1+y[0]),
				    y[0]*y[1]/(1+y[0]) - a*y[1],
				    a*y[1]));

  NumConstraint csp1(y,y[0]+y[1]+y[2] -10.0 = 0);
  NumConstraint csp2(y,y[0]>=0);
  NumConstraint csp3(y,y[1]>=0);
  NumConstraint csp4(y,y[2]>=0);

  Array<NumConstraint> csp(csp1,csp2,csp3,csp4);

  ivp_ode problem = ivp_ode(ydot,0.0,yinit,csp,AUTODIF);

  simulation simu = simulation(&problem,100.0,KUTTA3,1e-6);

  simu.active_monotony();

  simu.run_simulation();

  vibes::beginDrawing ();
  vibes::newFigure("PDS");
  vibes::setFigureProperty("plot","time",10);
  vibes::setFigureProperty("plot","y",10);
  vibes::setFigureProperty("plot","width",600);
  vibes::setFigureProperty("plot","height",600);

  plot_simu(&simu);

  vibes::endDrawing();

  return 0;

}