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 and Limitations

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

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:

List of publications on DynIBEX


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 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 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


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

Some applications