Application of interval analysis

Categories: solving a system of inequalities

Applications:: verification of the controllability of a system

## 1  Introduction

An interval for a variable x will be denoted x = [xInf,xSup] and represent the set of possible values for x such that xInf £ x £ xSup. The width of an interval x is defined as xSup-xInf.

A box for a set of n variables x1, ... xn is a set {x1, ..., xn} of intervals for all these variables. The width of a box is the largest width of all its interval.

An interval evaluation of a function f(x1,...,xn) of n variables x1, ... xn over a box for these variables is an interval f=[fInf, fSup] such that for all xi in xi
fInf <= f(x1,...,xn) <= fSup

There are different way to determine an interval evaluation of a function f. The simplest is the natural evaluation which consist to substitute any mathematical operator appearing in f by its equivalent interval arithmetic operator. For example consider the function f defined by f=x2-x and the box for x defined by {x=[2,3]}. We have
f([2,3])=[2,3]2-[2,3]=[4,9]-[2,3]=[1,7]
The natural evaluation allows one to determine efficiently lower and upper bound for the real values of a function whatever are the mathematical operators appearing in the function. These bounds may be overestimated (i.e. there is no value of x such that f(x) is either fInf or fSup but the difference between the real lower and upper bound and the interval evaluation will decrease with the width of the intervals in the box.

There are numerous packages that implement interval arithmetics (for example BIAS/Profil). All of them allows to take round-off errors into account, i.e. the real lower and upper bounds of the function are guaranteed to be included in the computed interval evaluation.

## 2  The example

We consider a system in the variable x11,x12,x21,x22,x23,x24,x41,x42, q14, v4+ with x1=[x11,x12], x2=[x21,x22,x23,x24], x4=[x41,x42] defined by
 Ax1+ Bx2+ C=0 (1) Dx4+ E=0 (2) Fx4+ G >=0 (3) Hx1+ I(x4) >=0 (4)
where

A= æ
ç
ç
ç
ç
è
 -1 -1 0 0 1 0 0 -1 0 0
ö
÷
÷
÷
÷
ø

B= æ
ç
ç
ç
ç
è
 0 1 -1 1 1 0 0 0 0 -1 0 0 0 0 1 0 0 0 0 -1
ö
÷
÷
÷
÷
ø
C= æ
ç
ç
ç
ç
è
 0.01-q14 0.01 0 0.01 q14
ö
÷
÷
÷
÷
ø
D=(-0.81, 0.9, 0)     E=-v4++0.1

F= æ
ç
ç
ç
ç
è
 -0.09 -0.9 1 0.9 -1 0 1 0 0
ö
÷
÷
÷
÷
ø
G= æ
ç
ç
è
 -1 0.9 -0.01
ö
÷
÷
ø

I= æ
ç
ç
ç
ç
ç
è
q14
0.01+ Kx4
0
0
ö
÷
÷
÷
÷
÷
ø

H= æ
ç
ç
ç
ç
ç
ç
ç
ç
ç
ç
è
0
-1
0.01+ Kx4

-1
0.1-x4,1
1
0.1-x4,1

1 0
ö
÷
÷
÷
÷
÷
÷
÷
÷
÷
÷
ø
K= (
 0.09 0.9 0
)

By solving the system of 6 equalities (1,2) we reduce the problem to the solving of a system S of 6 inequalities fi(x11, x41, q14, v4+) <= 0 in the 5 unknowns x11, x41, q14, v4+.

We will consider two different problems:
• problem P1: shows the existence of at least one solution to S
• problem P2: determine an approximation A of the set P of all solutions to S and an index allowing to measure the difference between A and P

## 3  Basic algorithms

The algorithm will use a list of boxes L, each box in the list being a set of intervals {x11, x41, q14, v4+}. Bi will denote the i-th box in the list and n will be the number of boxes in the list. This list will be initialized with a box B0 in which we are looking for a solution of a problem.

To solve P1 we may use the following algorithm:
1. define i = 0
2. if i>n, then EXIT. There is no solution to P1 in the box B0.
3. compute the interval evaluation [fInfi, fSupi] of the fi's over the box Bi
4. if there is at least one interval evaluation such that fInfi>0, then i =i+1 and go to step 2.
5. if for all interval evaluation we have fSupi <= 0, then EXIT. Any value of the unknowns in the box Bi is a solution to P1.
6. if for all interval evaluation we have fInfi <= 0 and fSupi >= 0, then
1. if the width of all the intervals in the box is lower than a threshold e then i =i+1 and go to step 2
2. select the unknown xj in the box whose interval has the largest width.
3. define 2 new intervals xj1, xj2 for xj as xj1=[xjInf, (xjInf+xjSup)/2], xj2=[(xjInf+xjSup)/2,xjSup]
4. create the 2 boxes { x1, ..., xj1, xj+1, ..., xn}, { x1, ..., xj2, xj+1, ..., xn} and add them to L. Then n=n+2, i =i+1 and go to step 2
The principle of the algorithm is quite simple: for a given box we compute the interval evaluation of the inequalities. If the lower bound of one evaluation is positive, then the box cannot contain a solution to S, so we consider the next box in the list. If the upper bound of all the evaluation is negative or 0, then any point within the box is a solution. Otherwise if at least one evaluation has a negative lower bound and strictly positive upper bound, then we split the largest interval in the box and create 2 new boxes having the same interval than the initial box, except for the variable that has been bisected. This bisection is done only if the width of at least one interval in the box is larger than a given threshold e. We proceed until one solution has been found or until all the boxes have been processed. In the later case if the threshold e has been chosen very small, then we may assume that there is no solution to the problem.

To solve P2 we use almost the same algorithm. The differences are the following:
• at step 5 a solution box is stored in a result file and the algorithm does not exit. The volume of the solution box box is computed and this volume is added to a number VIn, initialized to 0.
• at step 6-a when a box has all its intervals of width lower than e the volume of the box is computed and this volume is added to a number VNeg, initialized to 0.
The respective value of VIn, VNeg is an indication of the quality of the approximation. Note that it is easy to modify the algorithm so that if for a given e the ratio VIn/VNeg is not satisfactory we may decrease the value of e and process only the boxes that have been neglected during a previous run.

## 4  Improvements of the basic algorithms

The scheme described in the previous section is the simplest numerical scheme based on interval arithmetics. But these algorithms can easily be improved.

### 4.1  Evaluation improvement

We have already computed the interval evaluation of the function f=x2-x for the range {x=[2,3]}. The result was [1,7] which is clearly overestimated. This is due that when computing the natural evaluation we consider each term and compute its bound independently from the other terms. Another example is that the evaluation of x-x for x in [-1,1] is not 0 but [-2,2]. A possible way to improve the evaluation is to use the gradient of the function. In our example the derivative of f is 2x-1. The interval evaluation of this derivative is:
2[2,3]-1=[4,6]-1=[3,5]
Hence whatever is x in the range [2,3] the derivative will always be positive. Hence we will compute the interval evaluation of f by setting the lower bound of the evaluation to f(2)=2 and the upper bound to f(3)=6. The interval evaluation is now [2,6] which is the correct answer. This process may be used for functions having multiple unknowns but has to be recursive: indeed fixing a variable to a constant value may change the interval for an already computed evaluation of the derivative of the function with respect to another unknown. This procedure has not been used in the example.

### 4.2  The 2B method

Assume that we are considering a system in two unknown x, y in which we have the inequality:
a x + b y +c x2y2 <=0
where a, b are positive reals. This inequality may be written as
x <=
-by-c x2y2
a

Let U=[UInf,USup] be the interval evaluation of the right hand side of this inequality computed for a given box {x,y}. If xInf > USup, then the inequality will never be satisfied and the box can be rejected. Otherwise if xSup > USup, then the inequality will never be satisfied for any x in [USup,xSup]: hence we may change the x interval of the box from [xInf,xSup] to [xInf,USup], thereby reducing the width of this interval. Note that changing the interval x will lead to a change in the interval U: hence the procedure may be repeated until no noticeable change is observed on x. We may then start the same procedure for the unknown y.

This principle may evidently be extended to equalities and to consider inequality such that the left hand side term is not directly a variable but an invertible function of one variable (e.g. x2).

### 4.3  The 3B method

Consider a system of inequalities fi in the unknowns x1, ... xn and a particular box Bi= {x1, ..., xn}. Assume now that you change one interval in this box, e.g. the x1 interval, from [x1Inf,x1Sup] to [x1Inf,x1Inf+a] where a is supposed to be small. Proceed then to the interval evaluation of the fi for this new box. If the lower bound of one of the fi is strictly greater than 0 you have determined that there is no solution of the system for this box: you may then change the interval x1 in Bi from [x1Inf,x1Sup] to [x1Inf+a,x1Sup], thereby reducing the width of this interval. You may then repeat the procedure for the new box Bi.

Clearly you may also change the interval x1 on the right side (i.e. changing it from [x1Inf,x1Sup] to [x1Sup-a,x1Sup] and repeat the procedure for the other variables. Note also that if some intervals have been changed during the procedure it may be interesting to repeat it.

## 5ALIAS

ALIAS is a C++ library using the interval arithmetics of the package BIAS/Profil. An originality of this library is that it is partially interfaced with Maple. A specific Maple library allows to generate automatically the C++ code of the algorithms described in this document (including the 2B and 3B methods), to compile it, to execute a C++ program and to get the result without having to leave Maple. The results presented in the following section have been obtained using this process. An example of Maple use is the following, assuming that INEQ is the list of inequality and VAR the name of the unknowns
```
# automatic generation of the C++ code for the 2B method

HullIConsistency(INEQ,VAR,"SIMP_BROGLIATO");

#to indicate that we want to stop as soon as a solution has been found

`ALIAS/stop_first_sol`:=1:

#to indicate that we are looking for a solution with the width of the
#range of order 1

`ALIAS/epsilon`:=1:

#now solve

SOL:=GeneralSolve(
INEQ, #the set of inequalities
VAR, # the unknowns
[
[0.001,10], #search range for x11
[0.01,10], #search range for x41
[-10,10], #search range for x43
[-10,10], #search range for q14
[-10,0.91] #search range for v4p
],
"SIMP_BROGLIATO" #2B simplification
);
#SOL is the first found solution
```

## 6  Some results

### 6.1  Problem P1

Finding a solution of problem M1 with a width of the solution box of order 1 using the Maple interface takes a total time of 6.51s on a Sun Blade workstation. Here we use only the 2B method as simplification rule. This time include the generation, compilation and execution of the C++ code. The first solution found is:
```
[.63080396924504, 1.2536673076923], [2.5075, 3.131875],
[-5.69437499925, -4.5454999995], [-10, -9.375],
[-10,-9.4255625000676]
```
Note that it is possible to search for a solution within another search space without having to restart the code generation and compilation. The executable created during the first run can then be re-used with another search space. In that case the computation time for finding the first solution is 0.08s.

### 6.2  Problem P2

Using the basic algorithm and neglecting the boxes with a width of 3 ALIAS find in 9.13s solution boxes that have a total volume VIn of 12771.1 while the total volume of the neglected boxes is VNeg=195823 (ratio VIn/VNeg=0.0652)

Using the 2B method we get in 9.92s VIn=13737.5, VNeg=159872 (ratio VIn/VNeg=0.08592). Using both the 2B and 3B method we get VIn=12624.1 and VNeg=132142 (ratio VIn/VNeg=0.0955). With the same algorithm and with e=2 we get in 1m35.53s VIn=37652.9 and VNeg=85035.4 (ratio VIn/VNeg=0.4427).

Finally for e=1 we get VIn=61915.3, VNeg=48135.5 (ratio VIn/VNeg=1.2862) in a computation time of 29m15.81s
This document was translated from LATEX by HEVEA.