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) |
|
| |
|
| A= |
æ
ç
ç
ç
ç
è |
|
ö
÷
÷
÷
÷
ø |
|
| |
|
| B= |
æ
ç
ç
ç
ç
è |
| 0 |
1 |
-1 |
1 |
| 1 |
0 |
0 |
0 |
| 0 |
-1 |
0 |
0 |
| 0 |
0 |
1 |
0 |
| 0 |
0 |
0 |
-1 |
|
ö
÷
÷
÷
÷
ø |
C= |
æ
ç
ç
ç
ç
è |
|
ö
÷
÷
÷
÷
ø |
|
| |
|
D=(-0.81, 0.9, 0) E=-v4++0.1 |
| |
|
| F= |
æ
ç
ç
ç
ç
è |
| -0.09 |
-0.9 |
1 |
| |
| 0.9 |
-1 |
0 |
| |
| 1 |
0 |
0 |
|
ö
÷
÷
÷
÷
ø |
G= |
æ
ç
ç
è |
|
ö
÷
÷
ø |
|
| |
|
| I= |
æ
ç
ç
ç
ç
ç
è |
|
ö
÷
÷
÷
÷
÷
ø |
|
| |
|
| H= |
æ
ç
ç
ç
ç
ç
ç
ç
ç
ç
ç
è |
|
ö
÷
÷
÷
÷
÷
÷
÷
÷
÷
÷
ø |
K= |
( |
|
) |
|
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:
-
define i = 0
- if i>n, then EXIT. There is no solution to P1 in the box
B0.
- compute the interval evaluation [fInfi,
fSupi] of the fi's over the box Bi
- if there is at least one interval evaluation such that
fInfi>0, then i =i+1 and go to step 2.
- 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.
- if for all interval evaluation we have fInfi <= 0
and fSupi >= 0, then
-
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
- select the unknown xj in the box whose interval has the
largest width.
- define 2 new intervals xj1, xj2 for xj as
xj1=[xjInf, (xjInf+xjSup)/2],
xj2=[(xjInf+xjSup)/2,xjSup]
- 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
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.
5 ALIAS
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.