Do it Yourself!

The examples in this page are presented under the form of “labs” so that you can use them for practicing.

The complete source codes are available under the examples/ folder.

In these labs, we use Vibes to plot boxes but you can easily adapt the code to use your own graphical tool.

Fast instructions for installing and using Vibes are given here.

Set image

The complete code can be found here: examples/lab/lab1.cpp.

Introduction

The goal of this first lab is to calculate the image of a box by a function using interval arithmetic:

If we denote by f the function and [x] the initial box, then the set S to calculate is:

\[S := \{ f(x), x\in[x] \}.\]

Applying directly an interval evaluation of the function f to [x] will give a single box that only represents an enclosure of S.

To fight with the wrapping effect, we will split [x] into smaller boxes and evaluate the function with every little box as argument. This will result in a better description of the set that will eventually converge to S as the size of the boxes tend to zero.

This will consist in three tasks:

  • creating the function f
  • creating the initial box [x] and splitting it into small boxes
  • evaluating the function of every boxes
  • plotting the results

Question 1

Create in the main the function

\[\begin{split}f:(x,y)\in\mathbb{R}^2\mapsto \begin{pmatrix}\sin(x+y)\\\cos(x+0.9\times y)\end{pmatrix}.\end{split}\]
show/hide solution

Question 2

Create the box ([x],[y])=([0,6],[0,6]) and split each dimension into n slices, where n is a constant.

show/hide solution

Question 3

Evaluate the function on each box and plot the result with Vibes.

show/hide solution

Question 4

Compare the result with n=15, n=80 and n=500.

You should obtain the following pictures:

_images/lab1-15.png _images/lab1-80.png _images/lab1-500.png
n=15 n=80 n=500

Set inversion (basic)

The complete code can be found here: examples/lab/lab2.cpp.

Introduction

The goal of this lab is to program Sivia (set inversion with interval analysis) [Jaulin & Walter 1993] [Jaulin 2001], an algorithm that draws a paving representing a set E defined implicitely as the preimage of an interval [z] by a non-linear function \(f:\mathbb{R}^n\to\mathbb{R}\) (here n=2).

\[E:=\{(x,y)\in\mathbb{R}^2, \ f(x,y)\in[z] \}.\]
_images/sivia-basic.png

Sivia (basic variant). Result obtained with f(x,y)=sin(x+y)-0.1xy and [z]=[0,2], by simply alternating an evaluation and bisection phase. For a precision of \(\varepsilon=0.1\), the number of boxes generated by the algorithm is 11891.

The Sivia algorithm performs a recursive exploration of the initial box by applying the following steps:

  • inner test: if the image of ([x],[y]) by f is a subset of [z], the box is painted in green;
  • outer test: if the image does not intersect [z], the box is painted in red;
  • if none of these test succeeds and if ([x],[y]) has a maximal diameter greater than \(\varepsilon\), the box is split and the procedure is recursively called on the two subboxes.

Question 1 (Initialisation)

Create the Function object that represents

\[(x,y)\mapsto \sin(x+y)-0.1\times x\times y.\]

and the initial bounding box ([-10,10],[-10,10]).

show/hide solution

Question 2 (Initialisation)

We shall use a stack for implementing the recursivity. This stack is a container that will be used to store boxes.

Create a C++ stack and set the precision of bisection to 0.1.

Push the initial box in the stack. Define the image interval [z] and initialize it to [0,2].

show/hide solution

Question 3

Create the loop that pop boxes from the stack until it is empty. Define a local variable box to be the current box (the one on top of the stack).

Hint: use the top() and pop() functions of the stack class.

show/hide solution

Question 4

Implement the inner test (see above).

Hint: use is_subset.

show/hide solution

Question 5

Implement the outer test (see above).

Hint: use intersects.

show/hide solution

Question 6

If none of these test succeeds, split the box. We will split the box on the axis of its largest size. Finally, the two subboxes are pushed on the stack.

Hint: use extr_diam_index and bisect.

show/hide solution

Set inversion (with contractors)

The complete code can be found here: examples/lab/lab3.cpp.

Introduction

We will improve the Sivia algorithm by replacing in the loop the inner and outer tests by contractions. This leads to a more compact paving and a smaller number of boxes (see figure below).

The first part of the code is unchanged:

int main() {
  vibes::beginDrawing ();
  vibes::newFigure("lab3");

  // Create the function we want to apply SIVIA on.
  Variable x,y;
  Function f(x,y,sin(x+y)-0.1*x*y);

  // Build the initial box
  IntervalVector box(2);
  box[0]=Interval(-10,10);
  box[1]=Interval(-10,10);

  // Create a stack (for depth-first search)
  stack<IntervalVector> s;

  // Precision (boxes of size less than eps are not processed)
  double eps=0.1;

  // Push the initial box in the stack
  s.push(box);
  ...

The idea is to contract the current box either with respect to the constraint

\[f(x,y)\in[z],\]

in which case the contracted part will be painted in red, or

\[f(x)\not\in[z],\]

in which case the contracted part will be painted in green.

Given a contractor c, the contracted part is also called the trace of the contraction and is defined as \([x]\backslash c([x])\).

_images/sivia-full.png

Sivia (with contractors). Result obtained with f(x,y)=sin(x+y)-0.1xy and [z]=[0,2]. For a precision of \(\varepsilon=0.1\), the number of boxes generated by the algorithm is 5165.

Question 1

Build forward-backward contractors for the four constraints (see the tutorial):

\[\begin{split}f(x)<0, \quad f(x)\ge0, \quad f(x)\le2 \quad \mbox{and} \quad f(x)>2.\end{split}\]
show/hide solution

Question 2

Thanks to the composition, build a contractor w.r.t. \(f(x)\in[0,2]\).

Similarly, thanks to the union, build a contractor w.r.t. \(f(x)\not\in[0,2]\).

show/hide solution

Question 3

Create the function contract_and_draw with the following signature:

void contract_and_draw(Ctc& c, IntervalVector& box, const char* color);

This function must contract the box box in argument with the contractor c and plot the trace of the contraction (see above) with Vibes, with the specified color color.

Hints: use the diff function of IntervalVector to calculate the set difference between two boxes.

show/hide solution

Question 4

Replace in the loop the inner/outer tests by contractions.

show/hide solution

Set Inversion (using “Sets”)

The complete code can be found here: examples/lab/lab4.cpp.

Introduction

The purpose of this exercice is just to get familiar with the structure proposed in Ibex for representing sets (or pavings).

The set inversion is naturally one of the main features proposed in this part of the library. We will solve the same problem as before but this time with the Set class directly. This will take only a few lines of code.

Give first a look at the documentation on sets.

Question 1

Create the function \((x,y)\mapsto \sin(x+y)-0.1\times x\times y\) and the forward-backward separator associated to the constraint

\[0\le f(x,y) \le 2.\]
show/hide solution

Question 2

Build the initial set [-10,10]x[-10,10] and contract it with the separator with a precision of 0.1.

show/hide solution

Question 3

Plot the set with Vibes using a SetVisitor.

Solution: copy-paste the code given here.

Parameter Estimation

The complete code can be found here: examples/lab/lab5.cpp.

Introduction

This exercice is inspired by this video.

The problem is to find the values of two parameters \((p_1,p_2)\) of a physical process that are consistent with some measurements. Measurements are subject to error and we want a garanteed enclosure of all the feasible parameters values.

The physical process is modeled by a function \(f_{p_1,p_2}:t\mapsto y\) and a measurement is a couple of input-ouput \((t_i,y_i)\). We assume the input has no error. The error on the output is represented by an interval.

The model is:

\[f_{p_1,p_2}:t \mapsto 20\exp(-p_1t)-8\exp(-p_2t).\]

We have the following series of measurements:

t y
1 [4.5,7.5]
2 [0.67,4.6]
3 [-1,2.8]
4 [-1.7,1.7]
5 [-1.9,0.93]
6 [-1.8,0.5]
7 [-1.6,0.24]
8 [-1.4,0.09]
9 [-1.2,0.0089]
10 [-1,-0.031]

Question 1

Build the function f as a mapping of 3 variables, p1, p2 and t.

show/hide solution

Question 2

Build two interval vectors t and y of size 10 that contain the measurements data (even if the input has no error, we will enter times as degenerated intervals).

Hint: build interval vectors from array of double.

show/hide solution

Question 3

Build a system using a system factory. The system must contain the 10 constraints that represent each measurements and the additional bound constraints on the parameters:

\[0\le p_1\le 1, \quad 0\le p_2\le 1.\]
show/hide solution

Question 4

Calculate the parameter values using set inversion (see lab n°4). You should obtain the following picture:

_images/param-estim.png
show/hide solution

Parameter Estimation (advanced)

The complete code can be found here: examples/lab/lab6.cpp.

Introduction

This lab is a follow-up of the previous one.

We now introduce uncertainty on the input variable t. However, we will try somehow to make our parameter estimation robust with respect to this uncertainty. This means that the values of our parameters should be consistent with our output whatever is the actual value of the input. Mathematically, we require p1 and p2 to respect the following constraints.

\[\forall i, \quad \forall t\in[t_i], \quad f(p_1,p_2, t) \in [y_i].\]

As a contractor-oriented library, Ibex does not provide quantifiers at the modeling stage. This means that you cannot write directly a constraint like this one. You have to build contractors and apply “quantifiers” on contractors. Read the documentation about contractors and quantifiers.

Question 1

Like in the previous lab, create the function and the vector of measurements.

show/hide solution

Then, define a constant delta_t (the uncertainty on time) and inflate the vector of representing input times by this constant:

show/hide solution

Question 2

We will use the generic constructor of CtcForAll. Create a bitset that will indicate among the arguments of the function f which ones will be treated as variables and which ones will be treated as quantified parameters.

show/hide solution

Question 3

Create the inner and outer contractor for the “robust” parameter estimation problem. Set the splitting precision \(\varepsilon\) of the parameter to tdelta/5.

Hints: follow the same idea as in the contractor variant of set inversion (see lab n°3). Also, we need here to deal with n inner/outer contractors. To perform the composition/union of several contractors, see here.

show/hide solution

Question 4

The pictures below show the results obtained for increasing values of deltat (0,0.1,0.2 and 0.3). As expected, the larger the error on input, the smaller the set of feasible parameters.

_images/param-estim-1.png
_images/param-estim-2.png
_images/param-estim-3.png
_images/param-estim-4.png
deltat=0 deltat=0.1 deltat=0.2 deltat=0.3
show/hide solution

Stability

The complete code can be found here: examples/lab/lab7.cpp.

Introduction

The goal of this lab is to cast a classical problem in control theory into a set inversion problem.

We have a dynamical system y(t) governed by the following linear differential equation:

\[y^{(4)}(t)+a y^{(3)}(t)+b y^{(2)}(t)+(1-b) y'(t) + a y(t)=0.\]

where a and b are two unknown parameters.

Our goal is to find the set of couples (a,b) that makes the origin y=0 stable. It is depicted in the figure:

_images/stability.png

Hint: apply the Routh-Hurwitz criterion to the caracteristic polynomial of the system.

show/hide solution

Unstructured Mapping

The complete code can be found here: examples/lab/lab8.cpp.

Introduction

A robot is moving in a rectangular area [-L,L]x[-L,L] (with L=2) and tries to build a map while avoiding obstacles. The map is precisely the shape of obstacles inside the area.

The only information we have is when its euclidian distance to an obstacle get smaller than 0.9. It receives an alert, a vector of measurements which contain its own position (named x_rob and y_rob in the code) and the position of the detected obstacle point (named x_obs and y_obs in the code). We also know that all the points that are less distant than 0.1 from the detected point belong to the obstacle.

The robot has a series of n=10 measurements. The goal is to build an approximation of the map using set intervals.

You can first copy-paste the data:

// positions of the robot
double x_rob[n]={2, 1.61803, 0.618034, -0.618034, -1.61803, -2, -1.61803, -0.618034, 0.618034, 1.61803};
double y_rob[n]={0, 1.17557, 1.90211, 1.90211, 1.17557, 0, -1.17557, -1.90211, -1.90211, -1.17557};
// positions of the obstacle point
double x_obs[n]={0.9, 0.728115, 0.278115, -0.278115, -0.728115, -0.9, -0.728115, -0.278115, 0.278115, 0.728115};
double y_obs[n]={0, 0.529007, 0.855951, 0.855951, 0.529007, 0, -0.529007, -0.855951, -0.855951, -0.529007};
_images/mapping.png
show/hide solution