Examples of application of

Examples presented in this documentation

Example 2

The following `Maple` program enable to get the equations of this
example.

with(linalg): #location of the linear actuator on the ground xa1:=0:ya1:=0:xa2:=10:ya2:=0:xa3:=3:ya3:=10: #location of the linear actuator on the end-effector xb1:=0;yb1:=0:xb2:=4:yb2:=0:xb3:=2:yb3:=2: for i from 1 to 3 do OA.i:=array([xa.i,ya.i]): CBr.i:=array([xb.i,yb.i]) od: rot:=array([[cos(teta),-sin(teta)],[sin(teta),cos(teta)]]): for i from 1 to 3 do CB.i:=multiply(rot,CBr.i): od: OC:=array([x,y]): for i from 1 to 3 do AB.i:=evalm(OC-OA.i): AB.i:=evalm(AB.i+CB.i): ro.i:=dotprod(AB.i,AB.i,'orthogonal'): ro.i:=simplify(ro.i): od: eq1:=ro1-50:eq2:=ro2-26:eq3:=ro3-25:The system admit the two solutions:

Here we may use the general solving algorithm (section 2.3), the general solving algorithm with the gradient (section 2.4), the general solving algorithm with the gradient and Hessian (section 2.5).

These algorithms are implemented in the test program
`Test_Solve_General1`, `Test_Solve_Gradient_General1`,
`Test_Solve_JH_General1`

Example 3

The following `Maple` program enable to get the equations of this
example.

with(linalg): xa1:=0:ya1:=0:xa2:=10:ya2:=0:xa3:=3:ya3:=10: xb1:=0;yb1:=0:xb2:=4:yb2:=0:xb3:=2:yb3:=2: for i from 1 to 3 do OA.i:=array([xa.i,ya.i]): CBr.i:=array([xb.i,yb.i]) od: rot:=array([[cos(teta),-sin(teta)],[sin(teta),cos(teta)]]): for i from 1 to 3 do CB.i:=multiply(rot,CBr.i): od: OC:=array([x,y]): for i from 1 to 3 do AB.i:=evalm(OC-OA.i): AB.i:=evalm(AB.i+CB.i): ro.i:=dotprod(AB.i,AB.i,'orthogonal'): ro.i:=simplify(ro.i): od: eq1:=ro1-50:eq2:=ro2-26:eq3:=ro3-25: en1:=eq2-eq1:en2:=eq3-eq1: sw:=solve({en1,en2},{x,y}): assign(sw): eq1:=numer(eq1): eq1:=simplify(eq1,trig): eq1:=convert(eq1,horner,[sin(teta),cos(teta)]):This equation has 2 solutions 0, -0.8067834380 (or equivalently 5.47640187). Here we may use the general solving algorithm (section 2.3), the general solving algorithm with the gradient (section 2.4), the general solving algorithm with the gradient and Hessian (section 2.5) or the procedure devoted to trigonometric equation (section 2.13).

These algorithms are implemented in the test program
`Test_Solve_General2`,
`Test_Solve_Gradient_General2`,
`Test_Solve_JH_General2`, while for the the specific procedure
devoted to solving trigonometric equations we may use the program
`Test_Solve_Trigo`. In the later case we may use or not the
procedure which compute bounds on the roots of the equation.

We may then investigate the computation time when we have a pretty
good idea of the solution. For example we assume that we have a root
in the range [5.3, 5.7]. We may now use Newton method
(see section 2.9)
implemented in
the test program `Test_Newton2`, Brent's method
(see section 2.8 and Ridder's method (see
section 2.7)
implemented in the test program `Test_Ridder2`.An interesting point is that although Newton is faster than the others
methods we may encounter convergence problem in this example. For
example if the initial guess is 5.85 (which is closer to the solution
5.476 than the second solution 0) Newton will converge toward the 0
solution.

Example 4

with(linalg): xa1:= -9:ya1:=9:xa2:=9:ya2:=9: xa3:=12:ya3:=-3:xa4:=3:ya4:=-13: xa5:=-3:ya5:=-13:xa6:=-12:ya6:=-3: xb1:= -3:yb1:=7:xb2:=3:yb2:=7: xb3:=7:yb3:=-1:xb4:=4:yb4:=-6: xb5:=-4:yb5:=-6:xb6:=-7:yb6:=-1: rot:=array([[cos(p)*cos(h)-sin(p)*cos(t)*sin(h) ,-cos(p)*sin(h)-sin(p)*cos(t)*cos(h), sin(p)*sin(t)],[sin(p)*cos(h)+cos(p)*cos(t)*sin(h), -sin(p)*sin(h)+cos(p)*cos(t)*cos(h),-cos(p)*sin(t)], [sin(t)*sin(h),sin(t)*cos(h),cos(t)]]): for i from 1 to 6 do OA.i:=array([xa.i,ya.i,0]): CBr.i:=array([xb.i,yb.i,0]): CB.i:=multiply(rot,CBr.i): od: cen:=array([x,y,z]): for i from 1 to 6 do AB.i:=evalm(cen-OA.i+CB.i): ro.i:=dotprod(AB.i,AB.i,`orthogonal`): ro.i:=simplify(ro.i,trig): od: for i from 1 to 6 do eq.i:=ro.i-60^2: od: er1:=eq2-eq1:er2:=eq3-eq1:er3:=eq4-eq1:er4:=eq5-eq1:er5:=eq6-eq1: sw:=solve({er1,er2,er3},{x,y,z}): assign(sw): en1:=numer(eq1):en2:=numer(er4):en3:=numer(er5): for i from 1 to 3 do en.i:=simplify(en.i,trig): en.i:=convert(en.i,horner,[sin(p),cos(p),sin(t),cos(t),sin(h),cos(h)]): od:We are looking for a solution in the domain:

The system has a solution which is approximatively:

Here we may use the general solving algorithm (section 2.3), the general solving algorithm with the gradient (section 2.4), the general solving algorithm with the gradient and Hessian (section 2.5).

These algorithms are implemented in the test program
`Test_Solve_General`, `Test_Solve_Gradient_General`,
`Test_Solve_JH_General` (in that later case Kantorovitch
theorem is applied only at level 1).

`ALIAS` offers tools that may be used to solve many problems and
interval analysis has a broad range of application. In this section we
will present examples of realistic problems that have been solved by
interval analysis. Some of them have been solved using directly the solving
algorithms of `ALIAS`, for others it has been necessary to write
additional C++ code. Other examples may be found in the HEPHAISTOS web
page.

The following example originated from the INRIA Geometrica project. It exemplifies that focusing on the problem at hand while keeping in view that interval analysis will be used allows one to develop a very efficient algorithm that is at the same time faster and safer (in term of quality of the result) than conventional approaches.

Consider 2 ellipses E1, E2 in the plane that are perfectly defined and the
so-called *bi-tangent* i.e. the lines that are tangent to both
ellipsis. There is 4 such lines but we assume that some criteria
allows to choose two bi-tangents L1, L2. L1 intersects E1 at point P,
E2 at and L2 intersects E1 at point R and E2 at point S. We construct
now the vectors **PQ**, **RS** and then the determinant
D=**PQ** **RS**. The problem is to determine the sign of D in a
guaranteed manner (this sign is then used for another algorithm and a
mistake in the sign has a very negative influence).

A classical approach to solve this problem is to consider that the 8 coordinates of the points P, Q, R, S are the unknowns, write 8 equations to determine these unknowns, solve the system and then determine the sign of D. But this approach has drawbacks:

- the system of equations is complex
- it is difficult to solve
*exactly*

Now we may note that the coordinates of the points are bounded as the points are included in the bounding box of the ellipses. Being given ranges for the unknowns we may compute the interval evaluation of the determinant. If the lower bound is positive or the upper bound is negative we may directly state the sign of the determinant. If not we bisect one the unknown coordinates: we has thus a list of boxes with 8 ranges in each box, one for each coordinate. For each box we check if the coordinates ranges allows the points P, Q, R, S to belong to the ellipsis and to the bi-tangents, otherwise we reject the box. Then we compute the interval evaluation of D: if the lower bound is positive we store the box in a special array of boxes, called the positive array and discard it from the list. Similarly if the upper bound is negative we store the box in the negative array and discard it from the list. When the list of boxes is exhausted we look at the number of elements in the positive and negative array: if one array has 0 element, then we have determined the sign of D. Otherwise we consider the array that has the lowest number of elements and restart the bisection procedure on this new list, eliminating boxes that do not describe points on the ellipsis or that are not on the bi-tangent. At the end of this process we may have eliminated all the boxes of the list or, in the worst case, we will have computed the coordinates of P, Q, R, S i.e. we will end up with the classical approach. Our test show however that this is seldom the case and the above procedure is much more faster than solving using the classical approach.

Still in some cases the sign of D cannot be safely computed: this may happen for example when some coordinates have very large values while other have very small one. Numerical round-off errors then prohibit the exact determination of the sign of D, a fact that is detected by interval analysis. But we may still in some cases solve the problem: for that purpose instead of using Cartesian coordinates we may use polar coordinates with different ranges for the unknowns in the determinant or we may use an extended arithmetics. In any case the algorithm is safe as it provides either a sign, which is guaranteed, or will state that the sign cannot be determined with the current arithmetics.

A robot kinematics example

We present here a simplified version of the problem (see [10] for a more detailed version).

Consider a plane in the 3D space, called *the base*, on which lie
6 points
whose coordinates in an arbitrary reference frame.

Consider another plane, called the *platform*, on which lie
similarly 6 points
. To this plane we attach a
frame, called the mobile frame. The coordinates of
the 6 points in the mobile frame are supposed to be known. These
6 points constitute a rigid body which has a position/orientation with
respect to the reference frame (i.e the location of any point of the
rigid body may be calculated as soon as we know the coordinates of a
point of this body in the reference frame and the rotation matrix that
relates the mobile frame to the reference frame).

Assume now that the distances between each pair is fixed and is known. The problem is to determine what are the possible position/orientation of the rigid body such that these distance constraints are satisfied.

This problem is very important in practice. The platform is the "hand" of the robot and the distances between the may be changed by using actuators. The only information we have in the robot is obtained through sensors that measure the distances between the . Solving the above problem allows one to determine where is the hand of the robot by using the sensory information (but this problem has also other application such as robot calibration).

This problem is called the *direct kinematics* problem. It has a
counterpart, the *inverse kinematics* problem: being given the
position/orientation of the rigid body computes the distances between
the pairs . This inverse kinematics problem has a
straightforward solution: each squared distance can be calculated as a
function of the parameters that describe the position/orientation of
the rigid body. Hence solving the direct kinematic problem is
equivalent to solve the system of equations of the inverse kinematics
in the position/orientation parameters.

This problem is considered as one of the most difficult robot kinematics problem. Before 1992 the only known method to solve this problem was to use the Newton scheme and there was no theoretical result on this problem. Afterward with the use of sophisticated mathematics it has been possible to show theoretical results:

- there will be at most 40 solutions, complex and real
- examples with 40 real solutions have been exhibited
- it is possible to reduce this problem to the the solving of a 40th order univariate polynomial

The fastest algorithm to solve this problem is the use of the Gröebner basis package of Faugère, Rouillier that allow to compute all the solutions in a computation time that ranges between a few seconds to less than 30 seconds.

Interval analysis allows to solve this problem. Basically we have
noticed that we may choose as unknowns the coordinates of 3 points on
the platform (called the *reference points*). Indeed the
coordinates of the 3 other points may be
obtained as linear combination of the 3 chosen points (i.e. they are
virtual points). Furthermore the only equations we have are only
distance equations: distances between the pairs and
distances between the reference points and the virtual points. Hence
the distance solving algorithm described in section 2.16
is convenient for this problem.

Trials have shown that the computation time ranges between 10 to 40
seconds, which is competitive with the algorithm of Faugère,
Rouillier. Furthermore in most practical cases we are not interested
in computing *all* the solutions but only the one that corresponds
to the actual position/orientation of the platform and as the direct
kinematics problem has to be solved almost continuously and the robot
speed is limited we may
determine a small search domain in which the current solution should
be located. In that case the interval analysis based algorithm is very
fast as it is sensitive to the search space (at the opposite of the
Gröebner approach that has to compute all the solutions and
then sort which one is the current solution and hence is working in a
constant time). The algorithm is also safe compared to the Newton
scheme that may converge toward a solution that is not the current
one. Indeed we may guarantee that the current solution is included in
the algorithm result. Furthermore if more than one solution is found
then it is better to stop immediately the robot as the result of the
calculation is used to control it.

The problem presented here is for a robot described in the previous
section. We have seen that the hand of the robot was moving with
respect to its base by using linear actuators that allow to change the
distance between the pairs of points . But in practice these
actuators have a limited stroke and consequently the distance between
the points has to lie within a given range:

Furthermore the actuator is attached at point by joints (typically ball-and-socket or universal joints) that have limited motion. Hence the platform may reach only a limited set of position/orientation that is called its

Assume now that the hand has to follow a time-dependent trajectory: it
is clearly important to verify that all the points of this trajectory
lie within the workspace of the robot.
The trajectory is defined in the following manner: the
position/orientation parameter are written as analytic function of the
time T (which is assumed to lie in the range [0,1] without lack of
generality). Using the solution of the inverse kinematics (see the
previous section) it is then possible to express the distance
as functions of the time. For the trajectory to lie in the workspace
we have to verify the 12 inequalities:

when T lie in the range [0,1]. As the analytical form of the position/orientation parameters may be arbitrary we are looking for a generic algorithm that may deal with such arbitrary trajectory.

This can easily be done with interval analysis
(see [9] for a detailed version). First we define the
trajectory in Maple and compute the analytical form of (and
of any other constraint that may limit the workspace of the robot). We
get a set of inequalities that has to be satisfied for any T in [0,1]
if the trajectory
lie within the workspace. The analytical form of these inequalities
are written in a file: this allow their interval evaluation for any T
range by using the `ALIAS` parser. Then the general solving
procedure of `ALIAS` may be used to determine if there is a T such
that at least one constraint is violated.

An important point is that the algorithm allow to deal with the uncertainties in the problem. A first uncertainty occurs when controlling the robot along its nominal trajectory. Indeed the robot controller is not perfect and there will be a positioning error: for a nominal value of the position/orientation parameters the reached pose will be where can be bounded. A second uncertainty source is due to the differences between the theoretical geometrical model of the robot and its real geometry. Indeed to solve the inverse kinematics we use the coordinates of the in the reference frame and of the in the model frame. In practice however these coordinates are known only up to a given accuracy: hence these coordinates for the real robot may have any value within given ranges. Hence the inequalities of the problem have not fixed value coefficients but interval coefficients. But this no problem for interval analysis as the general solving procedures may deal with such inequalities. Hence if the algorithm find out that all inequalities are verified for any T in [0,1], then this means that whatever is the real robot and the positioning error of the robot controller the trajectory followed by the robot will fully lie within the robot workspace.

Robot singularity analysis

In the robot presented in the previous sections we have seen that the
inverse kinematic relates the
the distances between the to the
position/orientation parameters of the
platform **X** by
. An important matrix is the
jacobian matrix of this relation
. Indeed the forces in the legs are related to the
forces/torques applied on the platform by

where is a matrix that is never singular. Hence for given forces/torques applied on the platform the force in each leg may be computed as a ratio whose denominator if . A configuration where is called a

Usually this workspace is defined in the following manner: a geometrical object that will include all possible location of a specific point of the platform and ranges for the parameters that describe the orientation of the platform.

The determinant of can usually be expressed as a function of the position/orientation parameters. Interval analysis is appropriate to design an algorithm that check if there is a singularity within a given workspace. For that purpose we will take a point inside the workspace and calculate the value of at this point. Let us assume that this value is positive. A singularity will occur within the workspace if we are able to find another point in the workspace where will be negative.

We may assume that we are able to calculate a bounding box of the workspace and that we are able to design a test that allow to determine if a box for the position/orientation parameters is fully included in the workspace, is fully outside or part of it is inside and part of it is outside.

With that it is easy to design an algorithm that determine if the workspace includes a singularity (see [12] for more details).

First we use interval arithmetics to determine bounds for the value of . If the lower bound of the interval is positive is positive, then there is no singularity within the workspace. Otherwise we bisect the box. For the new boxes we first test if the box is inside the workspace: if one box is fully outside the workspace we discard it. If a box is not fully inside we bisect the box. For a box that is fully inside we compute the interval evaluation of . If the lower bound is positive we discard the box. If the upper bound is negative, then for any position/orientation defined by parameters within the box ranges the determinant will be negative: hence the workspace includes a singularity and we may stop the algorithm. If all the boxes in the list have been processed, then there is no singularity within the workspace.

Note that it is easy to generalize the algorithm so that uncertainties
in the robot geometry are taken into account. A general version of
this algorithm has been implemented using `ALIAS` and Maple and is
able to deal with different robot architectures. Maple is used to
compute the analytical form of . The `ALIAS` parser is used
to perform the interval evaluation of the determinant.

Robot workspace analysis

In practice all the geometrical parameters are submitted to
constraints. For example the distances must satisfy
, the angle must be such that
where is a unit vector that describes the main orientation
of the ball-and-socket joint at and is the maximal
possible rotation of this joint. All these constraints imply that
the value of the elements of are restricted to lie in a
specific region embedded in the 6-dimensional space, this region being
called the *reachable workspace* of the robot as the points of
this region define all the position/orientation that can be reached by
the platform.

In some cases it is possible to find an analytical description of the workspace (for example if the orientation is fixed and only the distance constraints are taken into account) but as soon as all the constraints are taken into account such analytical form is difficult to obtain. Interval analysis allows to obtain an approximation of the workspace as a set of 6-dimensional boxes. It is only an approximation as the algorithm, based on a classical bisection method, discard boxes whose width is lower than a given threshold. At the same time it is possible to obtain a bound on the error between the volume of the approximation and the volume of the workspace. For more details see [11]. A variant of the algorithm allows to solve efficiently a very practical problem: is a given region included in the reachable workspace of the robot ?

Robot synthesis example

The robot presented in the previous sections has numerous advantages
but its performances are very sensitive to its geometry (i.e. to the
location of the points ). The previous sections have
presented *analysis* problems: analyze for a given robot what are
its performance. But we may also consider an even more complex problem
which is the *synthesis* problem i.e. find what should be the
geometry of the robot such that it satisfies some performance
criteria.
This is very complex issue but one which has clearly a large impact in
practice (for example when designing a robot whose load will be over 2
tons and whose accuracy should be better than a micrometer as
considered by the European Synchrotron Radiation Facility in
Grenoble).

Consider for example the following problem: what should be the geometry of the robot (i.e. the location of the and the limits on the ) such that a given set of position/orientation can be reached by the platform ?

The inverse kinematic analysis allows to obtain the constraints as a set of inequalities between the design parameters and the elements of the set . A classical approach to solve this problem will be to determine the set that minimize the cost function . But this approach has many drawbacks:

- we cannot guarantee to find even one solution
- it will give at most a few solutions to the synthesis problem while in general there will be an infinite number of solutions
- for a given solution the geometry of the real robot will be different from the theoretical solution due to the manufacturing tolerances. Hence we cannot guarantee that the real robot will really be able to reach the set

- the lower bound of the interval for one is positive: we discard the box
- the upper bound of the intervals for all is negative: the box is stored as a solution
- the lower bound of the intervals for all is negative while their upper bound is positive: we bisect the box at the variable having the largest range except if all ranges have a width that is lower or equal to twice the manufacturing errors, in which case we discard the box

Using this approach we get not only one solution but ranges for the solution: indeed each point of the solution boxes satisfy all inequalities which means that the theoretical robot will be able to reach all the elements of . But we may take also the manufacturing errors into account. Indeed for a range A= in a solution box we may consider the range AP= that is included in the previous range and exists as the ranges of a solution box have a width at least . If we choose as design parameter a point in AP, then for the real robot the design parameter value will be included in A. As this is true for all the design parameters, then we may guarantee that the real robot will be able to reach all elements of .

Using this approach we are able to find all the geometries such that
the real robot satisfy the workspace constraints. Now we may consider
another design criteria and use the same approach, that will lead to
another set of design parameters. Taking the intersection of this set
with the set obtained for the workspace we will be able to compute all
robot geometries that satisfy the new design criteria **and** the
workspace constraint.

Here again interval analysis is an elegant approach that allows to solve a very practical problem. This approach has been used by the HEPHAISTOS project in various industrial contracts.

A quantum mechanics example

The Kochen-Specker theorem is one of the more important in quantum mechanics. Its a application involves the solving of the following problem (that has been simplified for the sake of clarity).

Consider unit vectors named `1,2,..9,A,B,..Z` in a four
dimensional space. A *group* of vectors is a set of 4 vectors that
have to be orthogonal to each other and not identical or opposite. A
`m-sequence` is a set of `m` groups such that:

- each vector in the sequence has the color black or white
- in each group of the sequence there must be exactly 3 black vectors and one white vector

1234,4567,789A,ABCD,DEFG,GHI1,35CE,29BI,68FHwhich involves 18 vectors.

Our objective is to to make an exhaustive search of m-sequences. For
that purpose we use a software designed by B.D. McKay, called the *generator*, that provides
for any set of `a` vectors and any sequence with `b` groups
all the `b`-sequences that satisfy the color constraints and are
not isomorphic (for example changing the name of the vectors will lead
to the same sequence).
Then we consider each potential `b`-sequence and we must check if
there are at least one set of coordinates for the vectors such that
the orthogonality
constraints for each group are satisfied together with the
non-collinearity between the vectors.

Such verification is not so easy. Indeed the orthogonality constraints
for a group involves 6 equations: consequently for a `b`-sequence
with `a` vectors
we have 6`b` orthogonality equations, `a` equations
that indicate that the vectors are unitary. If we consider for example
the Cabello's system the number of equations is therefore already 74
for 72 unknowns.
Although the problem is algebraic the solving method in algebraic
geometry have difficulties to handle such number of equations (for a
sequence that is one of the smallest).

A second problem is that an exhaustive enumeration of all potential
`b`-sequences leads to a very large number of sequences. For
example for `a`=10 and `b`=12 the number of sequences is 197
885 058

Nevertheless the problem was convenient for interval analysis. Indeed the unknown coordinates are component of unitary vectors and hence will lie in the range [-1,1]. Furthermore if a vector is solution then its opposite is also solution: hence we may restrict one of the coordinates of the vector to lie in the range [0,1]. In the same manner any rotation applied simultaneously to all solution vectors will lead to another set of solution vectors. Hence we may choose arbitrary any group so that its vectors constitutes an orthogonal basis of the four dimensional space. The coordinates of the vectors in this group will then be (1,0,0,0), (0,1,0,0), (0,0,1,0) and (0,0,0,1). Furthermore if one of these vectors appears in another group, then the orthogonality constraint will imply that one coordinate of all vectors of this group (except the basis vector) will be 0.

To solve this problem we have then developed a first very fast
interval analysis algorithm that has been incorporated in the generator
and eliminates sequence that trivially will have no solution. For
example a sequence starting with `1234,1235` will be eliminated as
`5` should be `4` or -`4`.
This algorithm drastically reduces the number of generated sequences.
Then for the remaining sequences we have used a general solving
procedure of `ALIAS`.
Even with the above simplification the size of the systems are large:
for Cabello's system we have reduced the problem to 32 unknowns and 82
equations that can be solve in about 60 seconds on a laptop. But
eliminating sequence without solution is usually faster.
For
example the system

1234,1256,1378,129A,24BC,34DE,5DFG,7HIJ,8BFK,9ADE,CGHI,6EJKwith 59 equations and 26 unknowns is eliminated in 10 seconds.

Although this work has not been completed yet we have already shown that Cabello's system is the smallest sequence (all sequences with a smaller number of groups have no solution).

2018-07-25