Open access peer-reviewed chapter

Path Tracking of a Wheeled Mobile Manipulator through Improved Localization and Calibration

By Tao Song, Fengfeng (Jeff) Xi and Shuai Guo

Submitted: February 4th 2018Reviewed: June 18th 2018Published: November 5th 2018

DOI: 10.5772/intechopen.79598

Downloaded: 153


This chapter focuses on path tracking of a wheeled mobile manipulator designed for manufacturing processes such as drilling, riveting, or line drawing, which demand high accuracy. This problem can be solved by combining two approaches: improved localization and improved calibration. In the first approach, a full-scale kinematic equation is derived for calibration of each individual wheel’s geometrical parameters, as opposed to traditionally treating them identical for all wheels. To avoid the singularity problem in computation, a predefined square path is used to quantify the errors used for calibration considering the movement in different directions. Both statistical method and interval analysis method are adopted and compared for estimation of the calibration parameters. In the second approach, a vision-based deviation rectification solution is presented to localize the system in the global frame through a number of artificial reflectors that are identified by an onboard laser scanner. An improved tracking and localization algorithm is developed to meet the high positional accuracy requirement, improve the system’s repeatability in the traditional trilateral algorithm, and solve the problem of pose loss in path following. The developed methods have been verified and implemented on the mobile manipulators developed by Shanghai University.


  • mobile manipulator
  • localization
  • tracking
  • path following

1. Introduction

Recently, mobile manipulators have been used in various industries including aerospace or indoor decoration engineering, which requires a large workspace [1, 2, 3, 4]. The said mobile manipulator consists of an industrial manipulator mounted on a mobile platform to perform various manufacturing tasks such as drilling/riveting in aerospace industry or baseline drawing in decoration engineering. Wheeled mobile platforms with Mecanum wheels that can easily navigate through crowded spaces due to their omnidirectionality with a zero turning radius are commonly used. Path tracking is one of the important issues for mobile manipulators, in particular for performing manufacturing tasks. This chapter addresses this issue from the aspect of localization and calibration.

Localization is a key functionality of mobile manipulator in order to track and determine its position around the environment [5]. Many methods are proposed to address this issue [6]. They can be divided into two categories: absolute localization and relative localization.

Absolute localization relies on detecting and recognizing different features in the environment to obtain the position and posture. The features can be normally divided into two types: artificial landmarks and natural landmarks. Compared with the natural landmarks, artificial landmarks have advantages of high recognition, which will lead to high accuracy. There is no cumulative error problem when a localization method based on artificial landmarks is used. The key challenge is to identify and extract the needed information from the raw data of landmarks. For the relative localization, dead reckoning and inertial navigation are commonly carried out to obtain the systems’ position. It does not have to perceive the external environment, but the drift error accumulates over time.

Researchers have proposed several solutions, such as fuzzy reflector-based localization [7] and color reflector-based self-localization [8]. The main drawbacks of these two methods are that the anti-interference ability is poor and the computation is huge. To solve these problems pertaining to the localization method based on artificial landmarks, Madsen and Andersen [9] proposed a method using three reflectors and the triangulation principle. Betke and Gurvits [10] proposed a multichannel localization method with the three-dimensional localization principle and the least squares method. Because of the unavoidable errors in position and angle measurement of reflectors, the use of only a trilateral or triangular method will not achieve high accuracy [11]. Nevertheless, there are still many challenges for mobile manipulator working in industrial environments such as aerospace manufacturing or decoration engineering, which requires high maneuverability and high accuracy at the same time. The stationary industrial manipulator has high repeated localization accuracy, which the mobile manipulator cannot achieve. This chapter focuses on the improvement of path tracking of a mobile manipulator through enhanced localization combined with calibration. Calibration is required for the system kinematic model established based on nominal geometry parameter to improve motion accuracy. Muir and Neuman [12] proposed a kinematic error model for Mecanum wheeled platform and applied actuated inverse and sensed forward solutions to the kinematic control. Wang and Chang [13] carried out error analysis in terms of distribution among Mecanum wheels. Shimada et al. [14] presented a position-corrective feedback control method with vision system on Mecanum wheeled mobile platform. Qian et al. [15] conducted a more detailed analysis on the installation angle of rollers. An improved calibration method is presented in this chapter to improve the tracking accuracy of a mobile manipulator.

2. System modeling

The wheeled mobile manipulator, as shown in Figure 1, is built with a manipulator onto a wheeled mobile platform with four Mecanum wheels. This system aims to carry out fuselage drilling/riveting tasks at assembly stations in an adaptive and flexible way in the aerospace industry.

Figure 1.

The wheeled mobile manipulator for fuselage drilling/riveting.

This system needs to localize in real time during machining process. The global frame {XG,YG,ZG}is attached to the ground in the environment. The platform frame {XP,YP,ZP}is attached to the center of the mobile platform. The tool frame XMYMZMis attached to the tool of the manipulator. Two laser range finders are equipped and their frames XL1YL1ZL1and XL2YL2ZL2are attached to the left-front and right-rear of the mobile platform, respectively. The vision frame XVYVZVis attached to the industrial camera.

The position and posture of the system in the global frame can be defined as


where xGand yGare the positions in two directions, respectively, and θGis the azimuth angle, as shown in Figure 2.

Figure 2.

The position and posture of system.

3. Accuracy analysis of a wheeled mobile manipulator

3.1. Problem formulation

Figure 3 shows the kinematic model of the mobile platform with Mecanum wheels. According to kinematic analysis [16], the motion equation can be obtained as


Figure 3.

The kinematic model of a Mecanum wheeled mobile platform.

where V=VXVYω0T, D0=14R0R0R0R0R0R0R0R0R0l0+LR0l0+LR0l0+LR0l0+L, W=ω1ω2ω3ω4T.

VXVYω0Tis defined as the linear and angular speeds of the platform; L is the half distance from the front axle to the rear axle as shown in Figure 3; l0is the transverse distance from the wheel centers to the platform center line; R0is the radius of the Mecanum wheel; andω1, ω2, ω3, and ω4are angular velocities of the four wheels, respectively.

The Mecanum wheel and its roller are shown in Figure 4. The roller is fitted on the edge of the wheel at a certain angle (generally 45°) with its central axis. Each roller can rotate freely around the central axis. The wheel relies on the friction between the roller and the ground to move. The material of roller’s outer rim is usually rubber, which will deform under pressure between ground and the wheel.

Figure 4.

Mecanum wheel and roller.

Figure 5 shows the distribution of roller deformation. F and T are the force and driving torque on the roller, respectively. The radiuses of wheels reduce differently under different pressures. The deformation zone and the position of the wheel center change with the action of driving torque and the shifting [17]. Furthermore, with such deformation, R0and l0will change.

Figure 5.

Roller deformation.

As shown in Figure 6, based on the kinematic analysis and consideration on deformation of the roller, Eq. (2) can be revised as follows: [18].


where D1=14R1R2R3R4R1R2R3R4R1l1+LR2l2+LR3l3+LR4l4+L; here D1is defined by R1, R2, R3,and R4, which are individual radiuses of the four Mecanum wheels, respectively. In order to improve the system motion accuracy, the relative errors Ribetween Riand R0, and relative errors libetween liand l0(i=1,2,3,4) should be obtained to revise matrix D1(Figure 6).

Figure 6.

Motion model of the mobile platform.

3.2. Error modeling

By multiplying time t, Eq. (3) becomes a displacement equation as


where X=Vtis XYθT. θ=Wtis θ1θ2θ3θ4T. Individual geometric parameters including l1, l2, l3, l4, R1, R2, R3,and R4are variables in Eq. (4). The relative errors can be determined from


leading to


where ∆X=∆X∆Y∆θT, ∆B=R1R2R3R4l1l2l3l4T, and


With the least squares method, the geometric errors can be solved as


Riis generally defined as the tolerance resulting from manufacturing or assembly. With the deformation of the roller and lilimited to hhand jj(h=3mm, j=5mm) respectively, liand Rican be defined as:


where liand Riare in l0jl0+jand R0hR0+h, respectively, and r is a random number between 0 and 1. The angular speeds of all wheels are set to 20rad/s, and the time t is set to 1 s.

The rank of matrix Gin Eq. (7) is generally 3, so Gis a full rank matrix and ∆Bcan be obtained. The following steps can be carried out. First, the displacement error measurement is analyzed. In order to obtain ∆B, it is needed to obtain the displacement error matrix ∆Xfirst.

The system is moved in the YGdirection and stopped every 2 s to obtain its position and posture XkGYkGθkG. The principle of measurement is shown in Figure 7 and Xcan be calculated as


Figure 7.

The principle of measuring displacement.

According to Eq. (2), XT, YT, and θTare theoretical values which can be obtained as follows: XT=Vxt, YT=Vyt, θT=0. Finally, through the experiment, the displacement errors can be obtained.


Substituting Eq. (11) in to Eq. (7), ∆Bcan be determined.

The Monte Carlo analysis is applied to deal with 50 samples. While Xhas been given in Eq. (11), the corresponding Bshould also satisfy its own tolerance, i.e., 3Ri3and 5Li5.

The results are given in Figure 8. These data are averaged to form a new result: B=0.5411.2370.6051.0550.4580.6830.0480.683T.A process capability index is then used to estimate the value of B. Process capability refers to the ability of a process to meet the required quality. It is a measure of the minimum fluctuation of the internal consistency of the process itself in the most stable state. When the process is stable, the quality characteristic value of the product is in the range μ3σμ+3σ, where μ is the ensemble average of the quality characteristic value of the product, and σ is the standard deviation of the quality characteristic value of the product.

Figure 8.

The values obtained with the Monte Carlo method.

There are two kinds of situations for which the process capability index has to be solved. First, μ=Mas shown in Figure 9: USLis the maximum error of quality and Lslis the minimum error of quality; M=USL+LSL/2; μ is the average value of process machining; and Cpindicates the process capability index, and Cp=USLLSL/6σ.

Figure 9.


Second, μMas shown in Figure 10: Cpkindicates the process capability index, and Cpk=USLLSL/6σMμ/3σ. Only process capability index greater than 1 is valid. As 3∆Ri3and 5∆Li5, M=0and μM, the following results can be obtained: Cpk∆R1=1.58, Cpk∆l1=0.57, Cpk∆R2=1.06, Cpk∆l2=0.58, Cpk∆R3=1.54, Cpk∆l3=0.63, Cpk∆R4=1.17, andCpk∆l4=0.5.

Figure 10.


The values of R1, R2, R3, and R4are closer to the real values than the values of l1, l2, l3, and l4. Now, taking B1a=0.5411.2370.6051.0550000T, B1b=00000.4580.6830.0480.683Tin Eq. (5), it can be seen that X1=0.4119,40.0003T, X1b=000.0003T. According to the change rate relative to Eq. (10), the influence of R1,R2, R3,and R4is much bigger than that of l1, l2, l3,and l4. Although the values of l1, l2, l3, and l4are not accurate enough for the real values, B1is valid.

The interval analysis is also carried out. In Eq. (6), the value ofGis uncertain because of the change of ∆Riand ∆liin 33and 55, respectively. Now, the increment of ∆Riand ∆liis set to 1 mm; there are 7 values of ∆Ri, and 11 values of ∆li, so there are 74×114groups of combinations of matrix G.Take all the combinations of matrix Gto Eq. (6) to obtain ∆Band exclude the cases for which ∆Riand ∆liare not in 33and 55, respectively. Thus, 277 groups of ∆Bcan be obtained. As shown in Figure 11, the average values of R1, R2, R3, R4, l1, l2, l3, and l4, lead to a new B2=0.2331.4670.9320.8240.20.2130.0680.039T, which is close to B1.

Figure 11.

The values obtained in interval analysis.

The same method is used to solve the process capability index of ∆B2. One obtains: Cpk∆R1=3.34, Cpk∆l1=0.68, Cpk∆R2=1.22, Cpk∆l2=0.59, Cpk∆R3=2.5, Cpk∆l3=0.63, CpkR4=1.73, andCpkl4=0.6. The values of R1, R2, R3,and R4in ∆B2are closer to the real ones than the values of l1, l2, l3,and l4in ∆B2. As the influence of R1, R2, R3, andR4is bigger than that of l1, l2, l3,and l4, ∆B2is valid.

The result is verified as well. ∆B1is used to revise the parameters of D1in Eq. (2) as


By setting two sets of four wheel speeds W1=10π/2110π/2110π/2110π/21Tand W2=4π/74π/74π/74π/7T, the displacement errors can be computed as S1=D1D0W2tand S2=D2D0W2t. The results of the two correction methods are almost identical in theory.

The experiment with four movements is shown in Figure 12(a), while the actual movement is shown in Figure 12(b). It can be found that S1and S2are close to the measured displacement errors, so both ∆B1and ∆B2are satisfied to revise the matrix D1.

Figure 12.

(a) Schematic diagram of measuring and (b) the actual measuring.

4. Localization

4.1. System configuration

The localization component of mobile manipulator includes two laser range finders and a number of reflectors that are of cylindrical shape placed in the environment. Each reflector is covered by a reflective tape with high reflectivity (Figure 13).

Figure 13.

Localization system for the mobile manipulator.

4.2. Dynamic tracking

To achieve accurate localization, the system should have the ability to perceive external information through extracting the reflector features. In this research, as shown in Figure 13, there are n reflectors, and the feature of each reflector Bii=12nis extracted from the raw data. The feature extraction algorithm consists of three steps: (i) filtering and clustering, (ii) identification, and (iii) feature extraction.

The first step is filtering and clustering. The raw data obtained by each finder are a set of discrete data sequences γλii=12n. γ is the distance from the target point to the finder. ∅ is the polar angle. λiis the intensity value of the ithdata point. In the process of data analysis, the outlier points that are contaminated by noise will be filtered out.

The density of the collected data points is proportional to the distance from the target point to the laser range finder. To improve the efficiency of the feature extraction process, an adaptive clustering method is adopted as given by Eq. (12). Unless the distance between two data points is less than the threshold δ, these data points are clustered for one reflector.


where γi1is the distance value of the i1thdata point, is the angle resolution of the laser range finder, β is an auxiliary constant parameter, and σγis the measurement error. The values of parameters σγand β are given as 0.01 m and 10°,respectively.

The second step is identification. After clustering, the data set can be correlated to each observable reflector. Each reflector data set is then used to calculate the position of the reflector in the laser range finder frame [19]. Let λδbe the reflected intensity threshold and DDδD+Dδbe the diameter range of a reflector, where D is a nominal reflector diameter and Dδis the tolerance. The values of λδand Dδare selected based on the actual situation. Considering that Wcrepresents a set after clustering, i.e., Wc=γλii=mn, for each reflector, the measured diameter Dcis calculated as Dc=γn2+γm22γnγmcosnm, where γmand γnare the beginning and end data of a reflector set, respectively. When the set Wcsatisfies the following two conditions


then the set Wcis established for all reflectors.

The third step is feature extraction. The feature extraction algorithm of a reflector extracts its central position γL,BβL,Bin the laser range finder frame XLYLZL, where γL,Band βL,Bare the distance and azimuth of each reflector, respectively. In the process of extracting the feature of the circular reflector, only a small portion of the entire circle was scanned with noise, so the accuracy of the fitting result would be low if a general least square circle fitting method is used. Since the radius of the reflector is known, the known radius is used as a constraint to improve the accuracy.

First, the value of βL,Bis obtained from nmconsecutive data points of the reflector set


As shown in Figure 14, Mirepresents the ith data and OL2B¯is a line between the reflector center B and the laser range finder center OL2; the projected angle of the line OL2B¯in the laser range finder frame is LβB, while the angle between the line OL2B¯and the line OL2Mi¯is θi. The distance from Mito B is approximately the radius of the reflector, i.e., BMi¯=D/2.


Figure 14.

Extraction of the center of the cylindrical reflector.

Furthermore, as mentioned before, the position of the finder in the platform frame is xPL2yPL2. The position of reflector center B γPBβPBin the platform frame can be expressed as


The optimal triangulation localization algorithm based on angle measurement is carried out. A number of experiments were carried out on the reflector feature extraction algorithm before proposing a localization algorithm for the system. A single reflector was placed in the measuring range of the finder. The finder performed repeated scanning when it was placed at difference places. Two different feature extraction algorithms were tested to calculate the position of the reflector, and the results are shown in Figure 15. Algorithm A is the feature extraction algorithm proposed in this chapter, while algorithm B is the least square circle fitting method without a radius constraint. Apparently, the former one yields a better result.

Figure 15.

Effect diagram using different feature extraction algorithms.

After extracting the reflector center, the positions of all the reflectors G=γLBβBLii=12ncan be obtained and used to calculate the distance measurement range and the angler measurement range Rβof the reflector, as given below:


where γLBmaxG, γLBminG, βBLmaxG, βBLminG.

It can be seen from Figure 16 that the angle measurement accuracy is better than the distance measurement accuracy. Therefore, based on these results, this chapter proposes an optimal trilateral localization algorithm based on angle measurement. The idea is to use the azimuth angle βbof the reflector in the platform frame to find the optima γbby the cosine theorem. The global pose of the mobile manipulator is then calculated based on the triangular method. The algorithm details are given as follows.

Figure 16.

The position of single reflector in the laser range finder frame.

First, it is assumed that the finder can cover at least three reflectors at each position. After feature extraction, the positions of three reflectors B1, B2,and B3are calculated as γPB1βB1P; γPR,B2βPR,B2; and γPR,B3βPR,B3, respectively. In the global frame, these positions can also be obtained as B1Gx1y1; B2x2y2G; and B3x3y3G. Since they are measured from the same finder, three circles must intersect at the same point OMPG, and the value of OMPGxMPGyMPG,represents the position of the mobile platform, as shown in Figure 17.

Figure 17.

Schematic diagram of the localization algorithm.

According to the cosine theorem, the relations between the above variables can be expressed by the following equations:


In Eq. (21), the known parameters βPB1, βPB2,and βPB3are used to solve γPB1a, γPB2a,and γPB3a. Since the above equations are nonlinear, the steepest descent method is adopted. The three circle equations can be expressed as


The position OMPGxMPGyMPGof the system can be obtained in the global frame by solving the above equations. Since the actual position of the reflector in the global environment deviates from the theoretical position, these circles may not intersect at the same point. In order to minimize the difference between the calculated position of the system and the actual position, the least squares estimation principle is applied. Assuming that the coordinate of the reflector is BiGxiyii=12n, n is the number of reflectors detected by laser range finder. The position value OMPGxMPGyMPGof the system is calculated as




The posture of the system also includes an azimuth angle θGin the global frame. First, θiGis obtained from the ithreflector as


The azimuth angle θGof the system is the averaged value from all the reflectors, as in


The dynamic tracking algorithm is then carried out. Localization of the system is based on landmarks. The system needs to meet the following two conditions to complete real-time localization:

  1. The system requires real-time localization of the reflector in the environment.

  2. The localization system correctly matches the real-time observed reflectors.

During the matching process, the reflectors observed in real time are made to correspond to all reflectors in the last moment in the environment one by one to extract the effective reflectors [20]. The localization can be achieved.

During localization, owing to the fact that some of the reflectors are obscured by obstacles or confused with a highly reflective metal surface object, the loss of position of the system is observed. To address the above problems, a dynamic tracking algorithm is proposed as shown in Figure 18.

Figure 18.

Schematic diagram of the dynamic tracking algorithm.

After placing n reflectors in the global environment, the system actually observes q reflectors at the tth moment, and the coordinate value of the ith reflector in the platform frame is Mt,iγt,Biβt,Bi0iq. The sampling time of the laser range finder is 0.1 s. Therefore, the theoretical position value Mt1,j1γt1,Bjβt1,Bjof the jthreflector in the platform frame can be deduced by the position Ot1,Rxt1yt1of the system at the t1th moment and the position OGjxjyjof the jthreflector in the global environment.

If the difference between the theoretical value Mt1,jγt1,Bjβt1,Bjand the observed value Mt,iγt,Biβt,Biis less than η, then the ithobserved reflector is an effective reflector and is considered matched with the reference reflector Bjas follows:


Therefore, a set A of effective reflectors at the tth moment can be obtained, and A=Mt,oγt,Boβt,BoMt,iγt,Boβt,Bii=01q. Taking the mobile manipulator’s moving speed into account, the value of η depends on the actual situation. Through the use of the optimal triangulation localization algorithm based on angle measurement, the pose Ot,Rxtytof the mobile manipulator can be calculated at the tth moment.

4.3. Experimental verification

The experimental data are obtained by the LMS 100 laser range finder with a scanning range of 270° and an angular resolution of 0.25°. The experimental platform is shown in Figure 19. The outside of the reflector is wrapped by reflective tape. In the experimental environment, five reflectors are placed around the system. Their global coordinate values are (0,995); (0, 0); (0, 1774); (2905, −2449); and (3956, −2032), and the unit is mm, as shown in Figure 20.

Figure 19.

Experimental platform.

Figure 20.

Experimental environment.

The optimal triangulation method based on angle measurement is used for validation by the repeatability of the system. In the stationary state of the system, the environment is scanned by finders. Each result is indicated by a red dot in Figure 21. The repeatability obtained by the trilateral method is nearly 18 mm, while the repeatability of the optimal method is only 9 mm. It can be shown that the optimal method is better than the traditional method.

Figure 21.

Repeatability localization of the system at the same location.

The mobile manipulator moves in the direction of the arrow in Figure 19, and each time the system moves a certain distance, the localization system will perform an experiment, i.e., it will use the left rear finder to calculate the current position. An average of 30 samples is taken for each experiment.

Figure 22 shows the results of static localization accuracy. The maximum distance error is 18 mm and the maximum angle error is 2°, which satisfies the localization requirement of the system.

Figure 22.

System localization error.

The mobile manipulator moves in the designated route, and it needs to constantly calculate and record its own position in the moving process. As shown in Figure 23, the trajectory of the moving system based on the localization method is smoother.

Figure 23.

Tracking results.

This chapter demonstrates the feasibility of a tracking and localization algorithm for mobile manipulators. The following conclusions can be drawn from this study: (i) In the detection of a reflector in the laser range finder frame, the angle repeatability of the reflector is better than that of the distance repeatability based on the feature extraction algorithm; (ii) The repeatability localization accuracy using the optimal triangulation method based on the angle measurement is nearly 9 mm, which is better than that of the trilateral method; (iii) The localization error of the system is 18 mm, which satisfies the localization requirement of system. Improvements in the location method based on reflectors, such as optimizing the layout of reflectors and the map of reflectors selection strategy for localization, are still needed.

5. Summary

In this chapter, through analyzing the roller deformation of the Mecanum wheel, the changed parameters of the motion equation of mobile system are found. The relative variation of the parameters in the motion equation of the Mecanum motion platform is solved by Monte Carlo analysis and interval analysis. Using the relative variation of the parameters to revise the motion equation, the displacement errors in different spaces in theory are solved for and compared with the measured displacement errors. From the comparison, both the methods are found to satisfy the system’s requirement. Then, the feasibility of a tracking and locating algorithm for mobile manipulator is demonstrated. The following conclusions can be drawn from this study: (i) In the detection of a reflector in the laser range finder frame, the angle repeatability of the reflector is better than that of the distance repeatability based on the feature extraction algorithm; (ii) The repeatability localization accuracy using the optimal triangulation method based on the angle measurement is nearly 9 mm, which is better than that of the trilateral method; (iii) The localization error of the system is 18 mm, which satisfies the localization requirement of system. Improvements in the localization method based on reflectors, such as optimizing the layout of reflectors and the map of reflectors selection strategy for localization, are still needed.

The method in this chapter is also used in the research of MoMaCo (Mobile manipulator in Construction), which can draw baseline for architectural decoration engineering as shown in Figure 24. The application result also verified the effectiveness of the method.

Figure 24.

MoMaCo (mobile manipulator in construction).

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Tao Song, Fengfeng (Jeff) Xi and Shuai Guo (November 5th 2018). Path Tracking of a Wheeled Mobile Manipulator through Improved Localization and Calibration, Applications of Mobile Robots, Efren Gorrostieta Hurtado, IntechOpen, DOI: 10.5772/intechopen.79598. Available from:

chapter statistics

153total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

4WD Robot Posture Estimation by Radial Multi-View Visual Odometry

By Edgar Alonso Martínez-García and Luz Abril Torres-Méndez

Related Book

First chapter

Kinematic Performance Measures and Optimization of Parallel Kinematics Manipulators: A Brief Review

By Abdur Rosyid, Bashar El-Khasawneh and Anas Alazzam

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More about us