Needs of Industry & Applications

The Needs of Industry Task

The purpose of the Needs of Industry task was to study the relationship between the current state of the art in polynomial system solving and industrial needs in this area. The goal of this study was to collect the industrial requirements for polynomial system solving by means of the identification of key problems and necessary algorithms and implementations required for their solution. A report was produced at the end of the first year to guide the algorithms tasks, with a further supplement added at the end of the second year.

The main purpose of this chapter is, therefore, to describe our perception of the particular needs of European Industry in this field. This information, gathered mainly during the first two years of the FRISCO project, has strongly determined the choice of algorithms which have been developed and implemented during the project.

This chapter is divided into three parts. In the first we describe the procedure used to approach European companies, industries and researchers; in the second we present our perception, shown in an unified way, of how users from industries regard and manage polynomial system solving; and in the third we describe the most relevant cases we encountered.

Strategy for Contacting Companies/Industries/Researchers

The first step was to decide what strategy to follow in order to get the required input data from users dealing with polynomial system solving in companies, industries and scientific laboratories. The following strategy was decided upon:

• To prepare a questionnaire containing several questions related to how polynomial system solving is carried out by the user. The questionnaire, which can be found in Appendix A, is divided into five sections:
  1. the first one is devoted to obtaining information about the industry, company, etc where the researcher works and about the mathematical software used there;
  2. the second one tried to obtain information about how the equations describing the mathematical model of the problems under consideration were generated;
  3. the third one asked about the characteristics of the equations (type, number, unknowns, coefficients, ...);
  4. the fourth one, when dealing with a nonlinear system of equations, asked about the kind of solutions looked for (real, complex, ...) and about the existence of parameters;
  5. the last one was devoted to determining where the nonlinear solving takes place: in real time during the production process, in a simulation stage, etc.
• To contact companies, enterprises and R&D laboratories dealing with polynomial system solving, and to send them the questionnaire which would be completed by a member of the FRISCO project during a meeting with the researcher involved.
• To organise an Open Workshop on the needs of industry for polynomial system solving (Barcelona, October 23-25, 1996) where the different enterprises, companies etc would be invited to present their cases and where FRISCO would show its present and future capabilities for polynomial system solving. A short description of the workshop together with its final program can be found on pages 127-129 of the first edition of ``The SAC Newsletter''.

The information gathered from the answers to the questionnaire and from the contacts in the Open Workshop are described in section 8.1.2. About one hundred companies, enterprises and R&D laboratories have already been contacted. Most of them were from Spain and France since, with the Open Workshop taking place in Barcelona, it was easier to persuade them to spend one or two days discussing their needs and experiences of polynomial system solving in ``real-world problems''. Another reason was that the main Spanish vendor of mathematical software, Addlink, was strongly involved in the preparation and organisation of the Open Workshop by providing the FRISCO project with a list of all the Spanish companies, enterprises and R&D laboratories using mathematical software such as Mathematica, Maple, Matlab, NAG Libraries, Axiom, etc. In the second year the activity was mainly concentrated in France and culminated in an industrial day which took place during the MEGA conference at St. Malo where several researchers from the FRISCO project and various companies exchanged their points of view about the role of polynomial system solving in industrial problems. A summary of this meeting, produced by B. Mourrain and M.-F. Roy, will appear in the ``SMAI Bulletin'' (SMAI: Societe pour les Mathematiques Appliquees et Industrielles).

There were three main difficulties found when contacting the companies, enterprises and R&D laboratories:

1. Time is a very valuable and expensive commodity for researchers working in industry. Thus it was very difficult to find enough time for a meeting, even if the researcher was strongly interested.
2. It is worth noting that the answers obtained started to be sufficiently detailed only after a third or fourth meeting.
3. Finally, as expected, researchers in academia and industry speak very different languages: we use the same words for the same kind of objects and problems but with very different meanings.

Summarising the answers

This section is devoted to presenting in an unified way the answers to the questionnaire obtained from the different companies, enterprises and R&D laboratories. It is organised according to the structure of the questionnaire.

Main Areas involved

The researchers who answered the questionnaire and were strongly interested in using the capabilities of FRISCO came from the following fields:

• Computer Aided Design (CAD), Computer Aided Modelling (CAM) and Computer Aided Engineering (CAE).
• Mechanical Systems Design.
• Signal Processing and Filter Design.
• Civil Engineering.
• Fluids, Materials and Durability.
• Robotics.
• Simulation.

In general, all these researchers arrived at a nonlinear problem since linearisation (which is, in general, the usual first step) was either not possible due to the characteristics of the problem, or did not provide accurate enough results.

The nonlinear problems to be solved appear mainly in two different ways: either as a nonlinear system of equations, or as an optimisation problem where the objective function is nonlinear or the feasible set is defined by some set of nonlinear conditions.

Mathematical Software used

Concerning the software used to deal with nonlinear systems of equations (their generation or their resolution), the following products appear in almost all the answers:

• Symbolic Software:
  • Mathematica and Maple
• Numerical Software:
  • Matlab, NAG Library and IMSL Library
• CAD/CAM:
  • AutoCAD and PRO-ENGINEER
• CAE:
  • ABAQUS, PRO-ENGINEER and Mathcad
• Development Tools:
  • C, C⁺⁺ and Fortran

In general, but not always, Symbolic and Numeric software such as Matlab is either used in a first stage during the preprocessing of the modelling equations, or as a checking procedure which provides a simulation of how the software which comprises the final solution should behave. It is worth mentioning that in several cases (for example, for the tolerance analysis in the design of a variational CAD/CAE system at LABEIN), some subroutines from the NAG or IMSL libraries are incorporated in the final product. In other cases a less sophisticated, but also useful, solution comes from using the `` Numerical Recipes'' book (in the design of a CAD system at CANDEMAT).

Characteristics of the polynomial systems of equations

In general, the nonlinear systems of equations to be solved have the same number of equations as unknowns, and over-determined systems appear quite rarely in practice. However, it is important to mention that the way the system of equations is generated is the main reason why no over-determined systems appear. One example of an over-determined nonlinear system arises in the filter design procedure of the CCETT Company.

The number of equations (and of unknowns) goes from 16 in the case of the mathematical modelling of the mechanical components of a shock absorber from LABEIN, to the several thousands arising from some nonlinear systems coming from the application of Finite Elements to solve some differential or integral equations (forging process modelling from LABEIN and nuclear plant simulation from TECNATOM) or from the resolution of a variational problem (bridge design from APIA XXI).

The equations appearing in the nonlinear system are usually quite sparse, in many cases being presented in a non-expanded form which produces this sparsity. The degree is not usually very big: either the total degree or the degree in every unknown is bounded by two or three.

With respect to the coefficients, these were real numbers in every example except the EdF case where the coefficients belong to $\mathbb {{Z}}$₃. In most cases they come from experimental data and thus are known only up to a limited accuracy. We also mention that, in the design of a variational CAD/CAE system at LABEIN, the coefficients are exact rational numbers.

Finally, parameters also appear quite often, but they are usually specialised in a first stage and, thus, the solution of the nonlinear system is obtained by interpolating the solutions of the specialised systems.

Algorithms and Solutions

The way the nonlinear system of equations is generally solved is by a method based on Newton Methodology (Newton-Raphson, Newton relaxation, gradient procedures for optimisation, predictor-corrector, etc). One of the reasons for this widespread use of Newton-like methods is that generally a good starting point can be given by experimental observations.

It is important to quote the CCETT case where the design of a new filter bank was accomplished by using Gb and RealSolving, two packages developed by members (J.-C. Faugere and F. Rouillier) of the FRISCO consortium.

There are two possibilities for the time to be spent on the resolution of the nonlinear system under consideration: the solution must be available in real time (in the inverse kinematics problem for a manipulator) or more time is allowed (when fitting the parameters of the theoretical model according to the experimental data available).

The solutions to be computed are always real (with the exception of the EdF example) and the best accuracy of the answers is around 10^-7. This arises from the design of a variational CAD/CAE system at LABEIN where the coefficients were exact rational numbers. It is a surprising fact that, in at least two cases (tolerance analysis at LABEIN and simulation at TECNATOM), the real-world problem implied that the particular nonlinear system of equations had only one real solution.

Explicit needs from the answers

This subsection is devoted to listing some of the concrete requirements stated explicitly in the completed questionnaires. They are very different in type and content due to their strong dependence on the particular area or problem they come from.

One of the most common request from companies, enterprises and R&D laboratories, apart from speed and accuracy, is the need of interactivity for the user and interconnectivity in the manipulation of different mathematical software packages. The provision of good-quality interfaces, together with the possibility of exporting functions to be integrated easily in specialised C⁺⁺ libraries, would greatly simplify the task of researchers in industry.

The availability of up-to-date methods for polynomial system solving would be very useful. For example, a nonlinear system of equations containing parameters is only solved by specialising the parameters several times and solving the corresponding specialised system, while some parametric computations (Gröbner Basis, resultants, etc) could be useful. Another example of this situation appears when convergence problems arise in the application of a Newton-like method: instead of replacing the algorithm by another one which is more powerful or better adapted (symbolic method, homotopy, ...), the mathematical model is simplified or changed.

A particular problem from a fluid structure interaction at LABEIN shows how the FRISCO capabilities can be useful in the near future. The problem is no more than the optimisation of an objective function whose analytical expression is not available and whose value is obtained when a concrete system of polynomial equations is solved (the input values of the function provide the coefficients of the system). This implies that an analytical expression of the function could either be computed exactly, or approximated by performing an elimination procedure.

Finally, and related to the previous point, a good collection of test suites, not only with statements of the problems and results, but also containing calling sequences allowing the user to experiment with every system on a WWW server, would be really useful. A preliminary version of auch a WWW server can be found at the FRISCO test suite site.

Case Studies

In this section we provide several concrete examples collected from the different companies and enterprises contacted by FRISCO. Some of them are presented together with their corresponding solution, or with the possible way FRISCO methodology could help solve them.

The shock absorber from LABEIN
Communicated by S. Casado (LABEIN)

The system under consideration has 27 equations and 27 unknowns with several parameters whose value is only known approximately. Two kind of unknowns are considered:

• the pressures: {P_i : 1 $\leq$ i $\leq$ 14}.
• the flows: {q_j : 1 $\leq$ j $\leq$ 15}.

The total degree of the equations is equal to one or two. They come from two different physical conditions: flow equilibrium in the different components of the electric-hydraulic circuit and pressure value in every node of such a diagram. The flow equilibrium equations are:

$$\matrix{ \hfill q_1&=&q_2\hfill \cr \hfill q_2&=&q_3+q_7\hfill \c... ...fill \cr \hfill q_6+q_{10}&=&q_{11}\hfill \cr \hfill q_{11}&=&q_{12}\hfill \cr}$$

and the pressure equations are:

$$\matrix{ \hfill P_2-P_1&=&q_1^2\displaystyle{{\rho}\over{2}}k_{\e... ..._{8}-S_9}\over{S_8^3}}+k_{\eirm gr}\displaystyle{1\over{S_8^2}}\bigg]\hfill\cr}$$

$$\matrix{\hfill P_{10}-P_9&=&q_{10}^2\displaystyle{{\rho}\over{2}}... ...{\eirm exp}\displaystyle{{(S_{14}-S_{13})^2}\over{S_{14}^2S_{13}^2}}\hfill\cr }$$

The initial data is given by the knowledge of the pressure and flow at the entry node:

q₁ = Q,    P₁ = 0

In this case, only real solutions are required.

The MESFET transistor
Communicated by A. Mediavilla (DICOM) and elaborated by L. Gonzalez-Vega

The design of the MESFET transistor requires the solution of the following nonlinear system of equations for n $\in$ {2, 3, 4}:

$$\sum_{j=1}^{n}$$a_ijcos(x_j + $$\beta_{ij}^{}$$) = b_i,        1 $$\leq$$ i $$\leq$$ n

This system arises when computing the coefficients linked to the derivatives in the Taylor development of the nonlinearities of the capacity appearing in the circuit model for the microwave transistor MESFET. These coefficients are very important when determining the transistor behaviour in the intermodulation distortion and, currently, almost no model takes care of this nonlinearity.

This system is converted to an algebraic one by expanding the cosine function and by introducing the new variables:

c_i = cos(x_i)    s_i = sin(x_i),        1 $$\leq$$ i $$\leq$$ n

In this way the algebraic problem is reduced to the resolution of the following polynomial system of 2n equations with 2n unknowns:

$$\matrix{ \displaystyle{\sum_{j=1}^n A_{ij}c_j+\sum_{j=1}^n B_{ij}s_j=b_i},& 1\leq i\leq n\cr\cr c_i^2+s_i^2=1,& 1\leq i\leq n\cr}$$

where the A_ij's and the B_ij's are polynomials with integer coefficients in the cos($\beta_{ij}^{}$), sin($\beta_{ij}^{}$) and a_ij. In this case, only real solutions are required.

The case n = 2 is very easily solved and the solution is given by the following:

$$\matrix{ U_4 s_2^4+U_3 s_2^3+U_2 s_2^2+U_1 s_2+ U_0=0\cr \cr c_2=... ...(s_2)}\over{V_0(s_2)}}\cr \cr s_1=\displaystyle{{V_3(s_2)}\over{V_0(s_2)}}\cr}$$

where the U_k's are polynomials in $\mathbb {{Z}}$[A_ij, A_ij, b₁, b₂] and the V_i(s₂)'s are polynomials in $\mathbb {{Z}}$[A_ij, B_ij, b₁, b₂][s₂]. For example, V₀(s₂) is:

$$\matrix{ 2(-B_{11}A_{21}+A_{11}B_{21}) (A_{11}^2A_{22}B_{22}s_2-A... ...2+B_{11}^2A_{22}B_{22}s_2+A_{12}B_{12}A_{21}^2s_2-A_{12}b_1 A_{21}^2)\hfill\cr}$$

Two interesting problems arise when dealing with these solutions:

1. The simplification of the polynomials U_i or V_j(s₂) would provide a better solution: for example in the previous expression of V₀(s₂) the first factor can easily be represented as
   - B₁₁A₂₁ + A₁₁B₂₁ = a₁₁a₂₁sin($$\beta_{21}^{}$$ - $$\beta_{11}^{}$$)

2. The solution of the cases n = 3 and n = 4 requires the manipulation of very big parametric expressions and one possible way of solving this problem could be the simplification question addressed before.

Finally, we note that the proposer of this problem considered the solution for the case n = 2 to be very important for his practical purposes, since he never imagined that such a solution existed.

The filter bank design case from CCETT
Communicated and elaborated by F. Moreau, J.-C. Faugere and F. Rouillier

The CCETT is a joint research centre between CNET (research organisation for France Telecom and TDF, a broadcasting operator, part of France Telecom Group). Its activities concern audio and video services, television and multimedia services, broadcasting techniques and so on. This includes research on signal processing and signal processing applications.

Description of the problem

Subband coding based on filter banks is a very efficient technique, but the design of the filter bank is still an open issue, since two main properties, orthogonality and linear phase, cannot be simultaneously achieved in the commonly used dyadic separable schemes, except for the Haar wavelet. Therefore designing efficient filter banks yielding simultaneously orthogonality and phase-linearity is a relevant issue in subband coding. It is possible to achieve these properties with the standard separable sampling scheme with non-separable filters. The design of such filters is however a difficult task. Three main approaches have been considered in the design of orthogonal non-separable wavelets:

• Straightforward formulation of the design as an optimisation of the filters' coefficients under the quadratic constraints for orthogonality (implying perfect-reconstruction) [9].
• Optimisation of cascade forms ensuring structurally orthogonality [4,5,6]. It is to be noticed that no complete cascade exists in the multidimensional case, due to the lack of a factorisation theorem. The main difficulty here appears when trying to optimise the parameters of the cascade.
• State-space matrix representations [7] look a very promising approach to the design of orthogonal multidimensional filter banks [8].

The design technique we propose derives from the cascade form approaches. It views the optimisation issue as solving a set of polynomial equations using the computer algebra techniques known as Gröbner bases.

We consider the design of M-band bidimensional FIR filter banks, i.e. M polynomials H₀,...H_{M - 1} exhibiting some desirable properties. A filter bank implements signal decomposition onto a basis, which we want to be orthogonal. It is also desirable in image processing to use linear-phase filters, which means (at least) centro-symmetric or centro-antisymmetric. It is well known that orthogonality and phase-linearity cannot be achieved for 2-band systems, expect in the case of Haar filters. Thus these properties cannot hold simultaneously for separable filter banks with sampling matrix 2I, which are the most common in image processing. This means that the simplest system allowing simultaneously orthogonality and phase-linearity is made of non-separable filters and sampling matrix 2I. Examples under cascade form are proposed in [4]. We consider a particular family of non-separable filter banks for sampling matrix 2I, holding structurally orthogonality and centrosymmetry. This family [4] is defined by polynomial matrix products. These matrices include some angles that can be chosen arbitrarily.

It is well known that such filter banks may generate wavelet bases [1], [3], and a necessary condition for that is that the filters vanish at some aliasing frequencies. In addition, the resulting wavelets can be N times continuously differentiable only if the polynomials H₀,...H_{M - 1} vanish as well as their derivatives up to order N at their aliasing frequencies. It is known for 1-D dyadic systems that in practice imposing these vanishing moments is an efficient way to obtain some regularity [2], so that we expect that designing maximally flat non-separable filter banks will provide regular wavelet bases.

The polynomial system of equations

The coefficients of polynomials H₀,...H_{M - 1} from the cascade in [4] are polynomials in terms of the cosines and sines of the angles. The flatness equations for the polynomials H₀,...H_{M - 1} are linear combinations of the coefficients of the polynomials H₀,...H_{M - 1}. This shows the principle leading to the use of techniques for solving polynomial systems in this signal processing context.

Actually the straightforward application of existing algorithms to the polynomial system obtained by writing the flatness equations in terms of the cascade parameters cannot design filters with support larger than 6x6 and 2 degrees of flatness. So to design filter banks with higher regularity we propose a change of variables, splitting the system into two smaller subsystems, which make it possible to design filters with support 16x16 and 5 degrees of flatness.

The cascade form we use in the sequel, borrowed from [4], and the formulation of the issue, maximising the flatness of filters of given size, turns it into a polynomial system.

Computer algebra tools used to solve the problem

For dealing with polynomial systems, we used recently-developed computer algebra tools which provide facilities for the computation of Gröbner bases (Gb, RealSolving), and related functionality (Hilbert function, lexicographic Gröbner basis computations by change of ordering, Triangular sets computations, Real Root counting, Rational Univariate representation of zero-dimensional systems). In addition, because of the complexity of the computations (size of the input/output), we designed a special change of variables and an adapted strategy to solve the problem.

Results and Conclusions

The most flat and regular non-separable wavelets ever designed have been computed. We give for each flatness order the minimal size for achieving it, the number of remaining free degrees, a parametrisation of the family of the minimal size filters for this flatness order, the optimisation of the remaining free degrees and the analysis of the resulting filter banks in terms of regularity, frequency characteristics and performance in image compression.

The results of our study concerns four scientific communities:

• From a computer algebra point of view, it is one of the major applications of computer algebra in digital signal processing. This area looks like yielding promising applications, especially in the field of filter banks. In addition, the problem has been solved using the most recent techniques for the Gröbner bases computation.
• From the point of view of the wavelet field, it is the first example of bidimensional orthogonal linear-phase wavelets with such regularity.
• From the point of view of filter bank design, the design examples are convincing: the resulting filters have good frequency characteristics, and it has been shown that the remaining free degrees can be chosen so as to optimise various criteria.
• From the image compression point of view, it is interesting to have new filter banks exhibiting simultaneously orthogonality, centrosymmetry and regularity, and in addition with good frequency characteristics. They give as good results as the best filters known up till know in terms of compression efficiency, without any global optimisation of the compression scheme.

Further work might include the study of other cascade forms for filter banks, or a proof that arbitrarily high flatness and regularity can be achieved in this filter bank family.

Bibliography

1

A. Cohen and I. Daubechies. ``On the instability of arbitrary biorthogonal wavelet packets''. SIAM J. Math. Anal., 24(5):1340-1354, 1993.

2

I. Daubechies. Ten Lectures on Wavelets. Number 61 in CBMS-NSF Series in Applied Mathematics. SIAM, Philadelphia, 1992.

3

J. Kovacevic and M. Vetterli. Non-separable multidimensional perfect reconstruction filter banks and wavelet bases. IEEE Transactions on Information Theory, 38(2):533-555, March 1992.

4

J. Kovacevic and M. Vetterli. Non-separable two- and three-dimensional wavelets. IEEE Transactions on Signal Processing, 43(5):1269-1272, May 1995.

5

D. Stanhill and Y. Y. Zeevi. Two-dimensional linear-phase orthogonal filter-banks and wavelets. pre-print, 1995.

6

D. Stanhill and Y. Y. Zeevi. Two-dimensional orthogonal wavelets with vanishing moments. pre-print, 1995.

7

S. Venkataranam and B. C. Levy. State space representations of @-d FIR lossless matrices. IEEE Transactions on Circuits and Systems-II: Analog and DIgital Signal Processing, 41(2):117-131, February 1994.

8

S. Venkataranam and B. C. Levy. A comparison of design methods for 2-d FIR orthogonal perfect reconstruction filter banks. IEEE Transactions on Circuits and Systems-II: Analog and DIgital Signal Processing, 42(8):525-536, August 1995.

9

E. Viscito and J. P. Allebach. The analysis and design of multidimensional FIR perfect reconstruction filter banks for arbitrary sampling lattices. IEEE Transactions on Circuits and Systems, 38(1):29-41, January 1991.

The ONERA-CERT case
Communicated by A.Burgueño & elaborated by L.Gonzalez-Vega and N.Gonzalez-Campos

ONERA-CERT ( Centre d'Études et de Recherche de l'École Normale Supérieure de l'Aéronautique de l'Espace) is a research laboratory at Toulouse under the auspices of the ``Office National d'Études et de Recherches Aerospatiales''.

One of the critical steps when they implement prototypes modelling slope-parametric hybrid automata (see F. Boniol, A. Burgueño, O. Roux and V. Rosu: Analysis of slope-parametric hybrid automata. Proc. of HART-97, Lecture Notes in Computer Science 1201, 75-80, 1997) is the resolution (real solutions) of a polynomial system of equalities and inequalities. This must be done once for every automata and it can be considered as a preprocessing step and, thus, the computing time it is not important. A first system sent is the one given by the system of non-linear inequalities

$$\matrix{ 0&\leq& a_1 \cr 0&\leq& a_2 \cr 0&\leq& t \cr 40&\leq& t... ...(K-1) \cr 40-40\,K&\leq& t\,(1-K)-a_2 \cr 30\,K-30&\leq& a_2\,K+t\,K(K-1) \cr }$$ (8.1)

The method used at ONERA-CERT to deal with this problem (a generalisation of the Fourier-Motzkin elimination method) provides the description of a set containing the solutions being sought but also contains other points which are not interesting. In this case the solution is presented in the following way:

$$\matrix{ a_1 &=& 0\cr a_2 &\geq& 0\cr t &\geq& 40\cr K &\leq& 1\cr}$$ (8.2)

But through the use of specific methods for Quantifier Elimination it is possible to describe explicitly and exactly the solution set as:

$$\left\{\vphantom{\matrix{ a_1=0\hfill\cr 0\leq K\leq3/4\hfill\cr t\geq 40\hfill\cr a_2=(40-t)(1-K)\hfill\cr}}\right. \matrix{ a_1=0\hfill\cr 0\leq K\leq3/4\hfill\cr t\geq 40\hfill\cr a_2=(40-t)(1-K)\hfill\cr} \bigcup \left\{\vphantom{\matrix{ a_1=0\hfill\cr K=1\hfill\cr t\geq 40\hfill\cr a_2=0\hfill\cr}}\right. \matrix{ a_1=0\hfill\cr K=1\hfill\cr t\geq 40\hfill\cr a_2=0\hfill\cr}$$ (8.3)

This solution provides much more useful information than the previous solutions when trying to use these automata for the development of communication protocols (see P.-H. Ho and H. Wong-Toi: Automated analysis of an audio control protocol. Proc. of CAV-95, Lecture Notes in Computer Science 939, 381-394, 1995). The solution set looks like:

[scale=.50]aufinal.ps

In the near future it is planned to integrate FRISCO Software and algorithms inside the toolboxes designed at ONERA-CERT to perform model-checking for slope-parametric hybrid automata and testing for temporised automata. This will be used to study an important industrial case: the Philips Communication Protocol.

The CANDEMAT case

Apart from the work accomplished directly for the specific three FRISCO applications described in the next section, a big effort was concentrated on continuing the successful contact with the CANDEMAT company initiated at the beginning of the FRISCO project. It is important to state, as main points, the following facts:

• One concrete person, Jesus Espinola, was hired jointly by the FRISCO project and the company CANDEMAT from July 1st 1998 to March 31st 1999 in order to work on the integration of the FRISCO algorithms and software into the CAD/CAM software CSIS developed (and distributed) by this company.
• In December 1998, a joint project entitled `` Integrating new algebraic-numerical techniques in Computer Aided Geometric Design (CAGD): Developing Problem Solving Environments and Implementation into an industrial CAD/CAM framework'' between the FRISCO team at the University of Cantabria and the company CANDEMAT was funded by the Spanish R&D agency (250,000 Euros). The project started on March 1 1999 and will run until December 31st 2001.
• The concrete problems considered during this year were the generic implicititation applied to sectioning and offsetting, the algebraically guided tracing of algebraic curves and surfaces implicitly defined, and the problem of conversion between the VDA and IGES formats.

Next we describe briefly some of the results obtained during the third year of the FRISCO project in cooperation with the company CANDEMAT.

Some details on the funded project between FRISCO team at the University of Cantabria and the company CANDEMAT

The abstract of the project entitled `` Integrating new algebraic-numerical techniques in Computer Aided Geometric Design (CAGD): Developing Problem Solving Environments and Implementation into an industrial CAD/CAM framework'' is the following:

This research project looks for, as a global objective, the development of several high level Problem Solving Environments specially devoted for CAGD (Computer Aided Geometric Design). This will allow to integrate, already evaluated into the Problem Solving Environments before mentioned, new algebraic and numerical techniques for CAGD into an industrial use and, in particular, in the CAD/CAM framework CSIS at the CANDEMAT company.

According to this approach, three different classes of objectives, well differentiated, can be isolated: first, the study and improvement of the graphical capabilities provided by the most frequently used high level software packages for Scientific Computing (Axiom/Aldor, Mathematica, Matlab, Maple) together with the development, into these packages or through a set of WWW servers, of several Problem Solving Environments specifically designed for Geometric Modelling and Visualisation; second, the development and integration of the new techniques in Polynomial System Solving produced inside the projects PoSSo (ESPRIT/BRA 6846: European Union) and FRISCO (ESPRIT/LTR 21024: European Union) to the resolution of different problems arising when dealing, in an very efficient way, with parametric curves and surfaces; and third, the resolution of a list of ``real-world'' problems to cover new requirements into the automotive industry and to improve already existing solutions into the CAD/CAM/CAE framework CSIS used by CANDEMAT for the design and production process of car pieces in the automotive industry.

If the goals previously described are achieved then a real transfer of basic research results, obtained under public funding (European Union, DGES), to the Spanish industry will be actually obtained: the impact of this technological transfer would be measured in terms of improvements into the production scheme at CANDEMAT through the optimisation and the enlargement of the capabilities of CSIS.

Generic Implicititation applied to sectioning and offsetting

The CAD/CAM CSIS software was prepared in order to support a database containing the different implicit equations of the most usual type of parametric surfaces appearing in the usual activity of the company. A concrete algorithm for sectioning was also implemented in CSIS but now using the implicit equation and this was compared with the previously used technique. The obtained result showed a considerably speed up when using the implicit equation.

In order to see how this technique works, we include a concrete example. The database, for this generic parametric surface,

$$\matrix{ x = x_{00}\displaystyle{t_2-t \over t_2-t_1}+x_{01}\disp... ...aystyle{s-s_1 \over s_2-s_1}\biggr]\displaystyle{t-t_1 \over t_2-t_1}\hfill\cr}$$

contains the implicit equation:

$$\matrix{ (z_{00}-2z_{10}+z_{11})xy +(y_{00}z_{10}+y_{11}z_{10}-y... ...y_{11}z_{00}+x_{00}y_{00}z_{11}-x_{00}y_{11}z_{10}-x_{01}y_{00}z_{10}\hfill\cr}$$

The surface to consider for sectioning is defined by a set of 9 patches given by the parameterisations:

x = p₁(s, t),    y = p₂(s, t),    z = p₃(s, t)

defined on the domain:

s₀ = s₁ < s₂ < s₃ < s₄ = s₅

t₀ = t₁ < t₂ < t₃ < t₄ = t₅

where

$$\matrix{ p_1=&\!\!\!\!\!(x_{00}f_0\!+\!x_{10}f_1\!+\!x_{20}f_2\!+... ...3)g_2\!+\! (z_{03}f_0\!+\!z_{13}f_1\!+\!z_{23}f_2\!+\!z_{33}f_3)g_3\hfill\cr }$$

and

$$\left.\vphantom{\matrix{ \left[\matrix{ f_0&=&{s_2-s \over s_2-s_... ...3 \over t_4-t_3} \cr }\right\}\hfill&\!\!\hbox{if $t_3 \leq t<t_4$}\cr}}\right.$$$$\matrix{ \left[\matrix{ f_0&=&{s_2-s \over s_2-s_0} \cr f_1&=&{s-... ...&=&{t-t_3 \over t_4-t_3} \cr }\right\}\hfill&\!\!\hbox{if $t_3 \leq t<t_4$}\cr}$$

The concrete values of the different parameters are:

$$\matrix{ x_{00}=-402.396210422345\hfill& \!y_{00}=125.60493907414... ...33\hfill& \!y_{33}=-29.8428535428032\hfill& \!z_{33}=51.4463427840486\hfill\cr}$$

and

$$\matrix{ s_0=s_1=0.0\hfill& s_2=0.333333333333333\hfill\cr s_3=0.... ...333333\hfill\cr t_3=0.666666666666666\hfill& t_4=t_5=2.14816198642715\hfill\cr}$$

In this case the generic implicit equation is:

c₁₀²c₂₂² + ((c₁₁² - 2c₁₀c₁₂)c₂₀ - c₁₀c₁₁c₂₁)c₂₂ + c₁₀c₁₂c₂₁² - c₁₁c₁₂c₂₀c₂₁ + c₁₂²c₂₀²

where:

c10=((s2-s0)*t0*t0*x01+(-s2+s0)*t0*t2*x00+((s2-s0)*t0*t2+(-s2+s0)*t0*t0)*x)*z11

+((-s2+s0)*t0*t2*x01+(s2-s0)*t2*t2*x00+((-s2+s0)*t2*t2+(s2-s0)*t0*t2)*x)*z10

+((-s2+s0)*t0*t0*x11+(s2-s0)*t0*t2*x10+((-s2+s0)*t0*t2+(s2-s0)*t0*t0)*x)*z01

+((s2-s0)*t0*t2*x11+(-s2+s0)*t2*t2*x10+((s2-s0)*t2*t2+(-s2+s0)*t0*t2)*x)*z00

+(((-s2+s0)*t0*t2+(s2-s0)*t0*t0)*x11+((s2-s0)*t2*t2+(-s2+s0)*t0*t2)*x10

+((s2-s0)*t0*t2+(-s2+s0)*t0*t0)*x01+((-s2+s0)*t2*t2+(s2-s0)*t0*t2)*x00)*z

c11=(((-2*s2)+2*s0)*t0*x01+((s2-s0)*t2+(s2-s0)*t0)*x00+((-s2+s0)*t2+(s2

-s0)*t0)*x)*z11+(((s2-s0)*t2+(s2-s0)*t0)*x01+((-2*s2)+2*s0)*t2*x00+((s2-s0)*t2

+(-s2+s0)*t0)*x)*z10+((2*s2+(-2*s0))*t0*x11+((-s2+s0)*t2+(-s2+s0)*t0)*x10

+((s2-s0)*t2+(s0-s2)*t0)*x)*z01+(((s0-s2)*t2+(s0-s2)*t0)*x11+2*(s2-s0)*t2*x10

+((-s2+s0)*t2+(s2-s0)*t0)*x)*z00+(((s2-s0)*t2+(-s2+s0)*t0)*x11+((-s2+s0)*t2

+(s2-s0)*t0)*x10+((-s2+s0)*t2+(s2-s0)*t0)*x01+((s2-s0)*t2+(-s2+s0)*t0)*x00)*z

c12=((s2-s0)*x01+(-s2+s0)*x00)*z11+((-s2+s0)*x01+(s2-s0)*x00)*z10

+((-s2+s0)*x11+(s2-s0)*x10)*z01+((s2-s0)*x11+(-s2+s0)*x10)*z00

c20=((s2-s0)*t0*t0*y01+(-s2+s0)*t0*t2*y00+((s2-s0)*t0*t2+(s0-s2)*t0*t0)*y)*z11

+((-s2+s0)*t0*t2*y01+(s2-s0)*t2*t2*y00+((-s2+s0)*t2*t2+(s2-s0)*t0*t2)*y)*z10

+((-s2+s0)*t0*t0*y11+(s2-s0)*t0*t2*y10+((-s2+s0)*t0*t2+(s2-s0)*t0*t0)*y)*z01

+((s2-s0)*t0*t2*y11+(-s2+s0)*t2*t2*y10+((s2-s0)*t2*t2+(-s2+s0)*t0*t2)*y)*z00

+(((-s2+s0)*t0*t2+(s2-s0)*t0*t0)*y11+((s2-s0)*t2*t2+(-s2+s0)*t0*t2)*y10

+((s2-s0)*t0*t2+(-s2+s0)*t0*t0)*y01+((-s2+s0)*t2*t2+(s2-s0)*t0*t2)*y00)*z

c21=(((-2*s2)+2*s0)*t0*y01+((s2-s0)*t2+(s2-s0)*t0)*y00+((-s2+s0)*t2

+(s2-s0)*t0)*y)*z11+(((s2-s0)*t2+(s2-s0)*t0)*y01+((-2*s2)+2*s0)*t2*y00

+((s2-s0)*t2+(-s2+s0)*t0)*y)*z10+((2*s2+(-2*s0))*t0*y11+((-s2+s0)*t2

+(s0-s2)*t0)*y10+((s2-s0)*t2+(s0-s2)*t0)*y)*z01+(((s0-s2)*t2+(s0-s2)*t0)*y11

+(2*s2+(-2*s0))*t2*y10+((-s2+s0)*t2+(s2-s0)*t0)*y)*z00+(((s2-s0)*t2

+(-s2+s0)*t0)*y11+((-s2+s0)*t2+(s2-s0)*t0)*y10+((-s2+s0)*t2+(s2-s0)*t0)*y01

+((s2-s0)*t2+(-s2+s0)*t0)*y00)*z

c22=((s2-s0)*y01+(-s2+s0)*y00)*z11+((-s2+s0)*y01+(s2-s0)*y00)*z10+((-s2+s0)*y11

+(s2-s0)*y10)*z01+((s2-s0)*y11+(-s2+s0)*y10)*z00

And the concrete implicit equation, after substitution, is:

$$\matrix{ H(x,y,z)=\quad-\Big[12.86405481303769\Big]z^3\hfill\cr\c... ...2+13792.839722405979x\hfill\cr\cr \quad\quad +10070450.541378483\Big]\hfill\cr}$$

The quality of the generic implicitations is first checked by replacing the parameterisation into the obtained implicit equation:

$$\matrix{H(p_1(s,t),p_2(s,t),p_3(s,t))= (-1.7970904234100502\cdot1... ...ot10^{\bf -9}s\hfill\cr\cr \qquad+3.7252902984619141\cdot10^{\bf -9}\hfill\cr}$$

but, in practice, it is enough to check a significant set of generated points into the surface to verify the implicit equation:

$${\displaystyle{\sum_{i=1}^N H(x_i,y_i,z_i)}\over N}$$ = 5.5134296417236327 ^. 10^-9

The sectioning by the plane X = k of a B-spline proceeds by performing the following steps:

• First, we consider all the implicit equations H_i(X, Y, Z) of the different patches integrating the considered B-spline surface.
• Second, the starting points in the parametric domain are determined and these are used to determine the interesting part of every curve H_i(k, Y, Z) = 0 through a convenient discretisation.
• The obtained points are interpolated by a spline curve which represents the desired section.

If the sectioning by using the implicit procedure has been a first successful application, this will be much more important when dealing with the offsetting by taking into account that the offset of a B-spline is no longer a B-spline or a parametric surface (but an algebraic surface) and thus the use of implicit equations is a direct method not currently used in any CAD/CAM software.

Format conversion: from rational to polynomial parameterisations.

VGA and IGES are two of the formats most commonly used in CAD/CAM systems to represent and deal with the data required to describe and communicate the essential engineering characteristics of physical objects such as manufactured products. The VDA (Verband der Automobilundustrie) format was created in 1982 by the VDA Committee founded mainly by several German automobile and automotive supply companies. The IGES (Initial Graphics Exchange Specification) was created by the IGES/PDES Organisation (USA).

The main difference between these two formats lies in the fact that VDA only accepts surfaces defined by polynomial parameterisations while IGES accepts both rational and polynomial. This is a very important problem in many companies in the automobile industry since they get the information in the IGES format but the specific CAD/CAM software they use only allows the use of the VDA format. The only way of solving this problem up to now has been by means of a method based upon an uniform subdivision of the parameter domain plus a polynomial interpolation procedure whose efficiency depends strongly on the degree of the polynomials to appear in required polynomial parameterisation (see Patrikalakis N., Approximate Conversion of Rational B-spline Patches, CAGD 6 (1989) 189-204). The Computer Algebra behind this problem is still not completely well understood and new algorithms for solving this problem have been developed and initially checked by using new Hermite-Birkhoff interpolation schemes for multivariate polynomials. These techniques also permit the reduction of the degree of the considered surface which produces a simplification of the implicitation procedure since the degree of the parametric equations are smaller.

Algebraically guided tracing of implicitly defined algebraic curves and surfaces

Many important problems in Computer Aided Geometric Design are reduced to the computation of the graph of a planar algebraic curve presented implicitly. For example if we want to section the surface

x = $${{X(s,t)}\over{W(s,t)}}$$,    y = $${{Y(s,t)}\over{W(s,t)}}$$,    z = $${{Z(s,t)}\over{W(s,t)}}$$;        s, t $$\in$$ [0, 1]

with respect to the plane X = x₀ then we have two possibilities: either we draw ``into the square unit'' [0, 1] x [0, 1] the planar algebraic curve defined by

x₀ = $${\frac{X(s,t)}{W(s,t)}}$$

and then this picture is lifted to the considered surface or, if the implicit equation H(x, y, z) of the considered surface is available, the lifting procedure can be avoided by merely computing the graph of the planar curve H(x₀, y, z) = 0 as it was shown in section 8.2.1.

The problem of computing the graph (even topologically) of a planar algebraic curve defined implicitly has received special attention from Computer Algebra since it has been responsible for many advances regarding sub-resultants, real root counting, infinitesimal computations, etc. These algorithms have been successfully implemented in AXIOM/ALDOR and now we are involved in the process of integrating the final algorithm coming from ALDOR into the CSIS software (but some new problems have arisen since CSIS has just moved to Visual Basic under Windows 95).

The main reason for including this algorithm in CSIS apart from its use in the sectioning procedure is the manipulation of trimmed surfaces: parametric surfaces with holes defined as parametric curves in the parametric domain of the surface. For example the sectioning of this kind of surfaces is not currently available in the CSIS software.

The Applications of the FRISCO project

FRISCO and Constraint Programming: The ILOG Solver (Deliverable 5.2.1)

A first version of a C⁺⁺ library allowing the use of the PoSSo library from the Ilog Solver has been developed and provided as deliverable 5.2.1. This task is combined with the experiment of an application based on the enhanced Ilog Solver that was demonstrated at the First Open Workshop. The sub-contractor Ilog has participated in the development of the library and brought some expertise on the use of the Ilog Solver.

The aim of this task was to provide a way to access the PoSSo library (PL) from the Ilog Solver (IS). The IS and the PL are both written in C⁺⁺ and originally it was foreseen that the PL would be simply linked to the IS. However there is no version of the IS for the g++ compiler and, when we began this work, the only PL version available was a g++ version. As a first experiment we tried to link the PL to the IS via the C interface but this method failed due to the difference of format of the object files generated by CC and g++. After investigating various options it was decided to use an abstract class for the communication part which can be derived into various communication subclasses (pipe, socket, direct linking). The current version of the library is using a client/server style interface which could be replaced by a direct library interface if a CC version of the PL became available. The current interface exchanges data via UNIX pipes between the IL and the PL which run as two independent processes.

The architecture of the library can be summarised by the following picture:

Within the IS, the main class needed for the communication with the PL is Groebner class. This architecture has been drawn from previous extension of the IS. The class Groebner contains three objects:

• GroebSys, the polynomial system
• PoSSoParam, a set of parameters sent to the PoSSoLib
• PoSSoSTUB, the communication class

GroebSys is made of a list of IlcGroebExp which are tree-like expressions to represent the constraints. The class PoSSoParam allows the user to set a list of parameters to be passed to the PL algorithms. The PoSSoSTUB class is an abstract class which can be derived into various communication subclasses (pipe, socket,...), like the PipeSTUB which allows a communication through UNIX pipes. The PoSSoSTUB class keeps track of the Groebner instances it is linked to. This could allow the computation to be distributed, or for there to be more than one PL server.

The library is based on the following software:

• for the PoSSo library: the PoSSoLib version 2.0 and g++ compiler (version 2.7.2). The PoSSo library has been compiled on SUNOS 5.5 (Solaris 2) and on Linux 2.0.X.
• for the Ilog Solver: the Ilog Solver version 3.2 and SUN CC version 4.X. The Ilog Solver has been compiled on SUNOS 5.5 (Solaris 2).

A documentation and reference manual for the library was included in the deliverable D5.2.1

FRISCO for Filters and Filter Banks (Deliverable 5.3.1)

This application uses the standard browser architecture to allow a user to use a number of different FRISCO components to solve problems connected to the design of filter banks for digital signal processing. While the browser architecture turned out not to be ideal for a number of reasons, it does provide an environment where a user who is not an expert in the FRISCO components can use them in an effective way.

The application considers six problems in the literature of Digital Signal Processing (DSP) and shows how the construction of filters and filter banks in these problems can be done by solving algebraic systems of equations. For each problem or model, there is a built-in sequence of computational steps where the user can choose from possibly several tools to carry out a particular operation.

FRISCO in Mechanism Theory and Signal Processing (Deliverable 5.4.1)

This application comprises a number of related parts, outlined below.

Forward kinematics of parallel robots

The forward kinematics problem of planar parallel robots using interval based algorithms has been implemented. It has been shown that they are real time and do no present convergence problem.

Accuracy estimation of parallel robots

For 6 degree of freedom parallel robots an algorithm for computing the extremal values of the positioning errors of the end-effector over a given workspace for given errors on the particular sensors has been implemented. Up to now there was no known method to solve this problem.

Workspace estimation of parallel robots

We have addressed the problem of determining the maximal workspace of a mechanism. Two algorithms have been designed and implemented for this problem. The calculation is still computer intensive but a by-product is an algorithm which enable to verify if a given 6D workspace lie within the reachable workspace of the robot in a few seconds.

Calibration of parallel robots

The purpose of calibration is to obtain a geometrical model of a robot close to the reality. Unfortunately calibration is very sensitive to noise measurements. We have tested methods based on the computation of eigenvalues of approximative matrices. Preliminary tests have shown that these methods are indeed less sensitive to noise.

Trajectory calculation of car suspension mechanisms

We are currently updating our methods to compute the trajectory of various suspension mechanism using three different approaches: interval analysis, certified Newton-Raphson and structured matrices. Preliminary tests have shown that these approaches are faster and more general.

Blind channel identification

In a recent paper we presented a method to identify a channel in presence of an unknown MSK modulated input by resorting only to output second order moments. This approach leads to a system of quartics which is solved with the help of an algorithm based on resultant and linear algebra techniques.

Conclusions

We consider the results arising from the Needs of Industry Task as one of the most important results achieved inside the FRISCO project. Apart from leading the research and software development inside the project it is important to remark that, at least, one concrete project with a real company (CANDEMAT SL) is already in progress due to the initial contacts and further work accomplished inside the FRISCO project.

At this moment another two projects, involving FRISCO partners, are under preparation: the first one dealing with the analysis and design of parallel manipulators (arising from the third application, Deliverable 5.4.1) and the second one devoted to the analysis of hybrid automata by using Quantifier Elimination algorithms (see 8.1.3).