SCALAR  Simultaneous Calibration of 2D Laser And Robot’s Kinematic Parameters Using Three Planar Constraints
Abstract
Industrial robots are increasingly used in various applications where the robot accuracy becomes very important, hence calibrations of the robot’s kinematic parameters and the measurement system’s extrinsic parameters are required. However, the existing calibration approaches are either too cumbersome or require another expensive external measurement system such as laser tracker or measurement spinarm. In this paper, we propose SCALAR, a calibration method to simultaneously improve the kinematic parameters of a 6\acdof robot and the extrinsic parameters of a 2D \aclrf which is attached to the robot. Three flat planes are placed around the robot, and for each plane the robot moves to several poses such that the \aclrf’s ray intersect the respective plane. Geometric planar constraints are then used to optimize the calibration parameters using LevenbergMarquardt nonlinear optimization algorithm. We demonstrate through simulations that SCALAR can reduce the average position and orientation errors of the robot system from 14.6mm and 4.05 to 0.09mm and 0.02.
lrfLRFLaser Range Finder \newacronymsvdSVDSingular Value Decomposition \newacronym[longplural=Degrees of Freedom]dofDoFDegree of Freedom \glsunsetdof
I Introduction
In traditional robotics applications such as pick and place, spraypainting and spotwelding, the robots either do not need very high accuracy or they are programmed by teaching, where the repeatability of the robot is more important than the accuracy. Repeatability refers to the robot’s capability to return precisely to the same location as previously taught, whereas accuracy refers to the robot’s capability to precisely reach a pose computed based on the robot’s kinematic model.
However, there are many applications where the accuracy of the robot becomes very crucial, given that the robot has to adapt to each task automatically with a great precision. Consider, for example, a robot drilling task in [1] where the robot is supposed to drill several holes at preciselydefined locations on a workpiece. The workpiece can be different for each task, and the placement within the workspace may not be precisely known. Programming by teaching in this case requires manual reprogramming for each workpiece which is very inefficient. To program the robot automatically for such task, the robot has to do a few things accurately: the robot has to scan the workpiece, determine the location of the holes, and finally move to that location accurately. The accuracy of such a robotic system depends on at least two things: The accuracy of the robot and the accuracy of the measurement system.
The accuracy of the robot is determined by how closely the kinematic parameters of the robot’s model resemble the actual kinematic parameters of the physical robot. This is affected by the manufacturing process, the assembly process, and the wear and tear during the operation of the robot. Robot kinematic calibration is usually conducted to achieve a better accuracy, either by using an external measurement system (such as motion capture system or coordinate measuring machines) or by constraining the motion of the endeffector.
The accuracy of the measurement system can be divided into two parts: the accuracy of the measurement device itself and the accuracy of the relative pose between the robot frame and the measurement frame. The accuracy of the measurement device depends on the type of device that is used and its specification. For example, a camera system is generallly less precise as compared to a laser system, although a camera can give more information. The second part of the accuracy comes from the fact that the data from a measurement system is always obtained in the measurement device coordinate frame, and it needs to be transformed to the robot coordinate system. Hence, the relative pose between the robot coordinate frame and the measurement device frame needs to be obtained. This relative pose is often called the extrinsic parameters of the measurement device, and the method to estimate the extrinsic parameters is called extrinsic calibration.
Researchers  Robot  Measurement Device  Initial Accuracy[mm]  Final Accuracy[mm] 

Ginani and Mota [2]  ABB IRB 2000  ROMER Measurement Arm  2.20  1.40 
Ye et al. [3]  ABB IRB 2400/L  Faro Laser Tracker  1.764  0.640 
Nubiola and Bonev [4]  ABB IRB 16006/1.45  Faro Laser Tracker ION  2.158  0.696 
Newman et al. [5]  Motoman P8  SMX Laser Tracker  3.595  2.524 
In this paper we present SCALAR, a calibration method to simultaneously improve the accuracy of the robot and the measurement system. SCALAR calibrates simultaneously the kinematic parameters of the robot and the extrinsic parameters of a 2D \acllrf (\aclrf) using only the information provided by the \aclrf attached to the robot endeffector. An \aclrf is chosen because it gives very accurate measurement data both for the calibration and for the subsequent tasks (such as drilling).
The overall calibration procedure is as follows:

The \aclrf is attached to the robot and three perpendicular planes are placed around the robot as shown in Fig. 2. Only the rough estimate of the position and orientation of the planes are necessary to be known, so the setup can be easily done.

For each plane, the robot is moved to several poses such that the \aclrf’s 2D ray is projected onto the plane. The data from the \aclrf as well as the robot’s joint angles information are collected.

An optimization algorithm is used to find the optimal robot’s kinematic parameters together with the \aclrf extrinsic parameters, using the geometric constraint that the projected \aclrf data should be located on the three planes.
The remainder of the paper is as follows. In Section II we discuss the existing approaches to the calibration problem both for the robot’s kinematic parameters and the \aclrf’s extrinsic parameters, and how SCALAR differs from the other approaches. In Section III, SCALAR is explained in detail. A simulation study is presented in Section IV to verify SCALAR, and finally we conclude with a few remarks in Section V.
Ii Related works
Iia Calibration of robot’s kinematic parameters
Robot kinematic calibration has been researched for a long time; some of the earliest works began in 1980s. The calibration procedures can be categorized into unconstrained and constrained calibration. In unconstrained calibration, the robot moves its endeffector to several poses while an external measurement system measures the pose. The measured pose is then compared to the one computed from the kinematic model, and the model is updated to minimize the difference between the predicted pose and the measured pose. In constrained calibration, some constraints are applied to the motion of the endeffector, and the constraints yield several calibration equations by which the robot’s kinematic parameters are calibrated.
Examples of unconstrained calibration works can be seen in Table I. The issues with such calibration method are the difficulty in setting up the calibration setup and the expensive cost of the external measurement system. For example, the cost of a laser tracker is more than USD [4]. Therefore, many researchers try to find calibration methods which only rely on the internal sensors of the robot by constraining the motion of the endeffector such as in constrained calibration.
In [6], Ikits and Hollerbach propose a kinematic calibration method using a planar constraint. A touch probe attached to the robot’s flange is moved to touch random points on a plane. When the touch probe is in contact with the plane, the joint angles are recorded. The kinematic parameters of the robot model are then updated to satisfy the planar constraint. While the approach is promising, they also report that some of the parameters are hardly observable when the measurements are noisy or when the model is not complete. The unobservability of the parameters means that some of the parameters cannot be obtained accurately from the calibration procedure.
In [7], Zhuang et al. investigate robot calibration with planar constraints, in particular the observability conditions of the robot’s kinematic parameters. They prove that a singleplane constraint is insufficient for calibrating a robot, and a minimum of three planar constraints are necessary. Using three planar constraints, the constrained system is proved to be equivalent to an unconstrained pointmeasurement system under three conditions: a) All three planes are mutually nonparallel, b) the identification Jacobian of the unconstrained system is nonsingular, and c) the measured points from each individual plane do not lie on a line on that plane. They verify the theory by doing a simulation on a PUMA560 robot.
In [8], Joubair and Bonev calibrated both the kinematic and nonkinematic (stiffness) parameters of a FANUC LR Mate 200iC industrial robot by using planar constraints, in the form of a high precision 9inches granite cube. The robot is equipped with an MP250 Renishaw touch probe, which is then moved to touch four planes of the granite cube. The granite cube’s face is flat to within 0.002mm. They improved the maximum plane error from 3.740mm to 0.083mm.
IiB Calibration of extrinsic 2D \aclrf parameters
Extrinsic calibration of an \aclrf consists of finding the correct homogeneous transformation from the robot coordinate frame to the laser coordinate frame. Most of the works on extrinsic calibration of an \aclrf involves a camera, since both sensors are often used together in many applications. The works in this field are largely based on Zhang and Pless’ work [9]. They propose a method to calibrate both a camera and an \aclrf using a planar checkerboard pattern. First, the camera is calibrated by a standard handeye calibration [10] using a checkerboard pattern. The calibrated camera is then used to calculate the pose of the pattern. Next, the robot is moved to several poses with the \aclrf pointing to the pattern. By using the geometric constraints that all the data points from the \aclrf should fall on the pattern plane, the extrinsic parameters of the \aclrf can be obtained. Finally, the same constraints are used to optimize both the intrinsic and extrinsic parameters of the camera and the extrinsic parameters of the \aclrf. The nonlinear optimization problem is solved by using LevenbergMarquardt optimization algorithm.
Unnikrishnan and Hebert [11] use the same setup as [9], although they do not optimize the camera parameter simultaneously due to the nonlinearity of the resulting cost function. Li et al. [12] use a specially designed checkerboard to calibrate the extrinsic parameters between a camera and an \aclrf, and claim that the result is better than [9]. Vasconcelos et al. [13] develop a minimal closedform solution for the extrinsic calibration of a camera and an \aclrf, based on the work in [9].
IiC Novelty of the proposed method
SCALAR can be seen as a combination of the algorithm for extrinsic calibration of an \aclrf [9] and the algorithm for calibration of robot’s kinematic parameters using three planar constraints [8]. It has the following advantages as compared to the other calibration approaches:

SCALAR simultaneously calibrates both the \aclrf extrinsic parameters and the robot’s kinematic parameters. Given that calibration process is often cumbersome, this saves a lot of time and effort. Moreover, the errors in the robot’s kinematic parameters affect the extrinsic calibration accuracy, and vice versa. Hence, calibrating both parameters simultaneously results in better accuracy.

SCALAR does not need an additional camera to calibrate the \aclrf, unlike [9].

SCALAR does not need another expensive external measurement system. The measurement is done using the \aclrf that will also be used in the subsequent robot task, hence it does not incur additional cost. Moveover, an \aclrf can achieve very high accuracy at much lower cost (more than ten times cheaper) as compared to the commonly used measurement systems such as Vicon or Faro Laser Tracker.

SCALAR does not need a precisely manufactured calibration object such as the granite cube in [8], which requires the planes’ position and orientation to be known accurately. SCALAR only requires three flat surfaces which are oriented roughly perpendicular to each other and the rough estimate of their locations. This also means that the calibration setup can be done easily.

The calibration poses can be distributed throughout the whole workspace, instead of being confined only in a local region such as in [8].
Iii Method
The calibration setup is depicted in Fig. 2, where three roughly perpendicular planes () are placed around the robot. An \aclrf is attached to the robot flange. For each plane, the robot is moved to poses such that the \aclrf’s ray is directed to the respective plane. One data set from the \aclrf consists of hundreds of data points, so data points are selected from the \aclrf data for each pose, and the robot’s joint angles are recorded.
This section describes the detail on how to calibrate both the robot’s kinematic parameters and the \aclrf’s extrinsic parameters. First, the initial estimate of the LRF’s extrinsic parameters is obtained using the linear leastsquares method with the data from one of the planes. This is based on the algorithm in [9], although presented differently for better clarity. After that, the robot’s kinematic parameters and the LRF’s extrinsic parameters are optimized simultaneously to satisfy the three planar constraints using LevenbergMarquardt nonlinear optimization method. Finally, we explain how \acsvd can be used to analyse which calibration parameters are identifiable, and the steps to handle the unidentifiable parameters are then presented.
Iiia Initial Estimate of the \aclrf Extrinsic Parameters
To obtain an initial estimate of the \aclrf extrinsic parameters, only the data from one plane is necessary. Arbitrarily, the bottom plane is chosen. The extrinsic parameters of the \aclrf , i.e. the homogeneous transformation from the robot flange coordinate frame to the \aclrf coordinate frame, is estimated by the following calculations.
Let the subscript/superscript , , and denote the coordinate frame of the robot base, the robot flange, and the \aclrf, while the subscript , , and refer to the \aclrf data point index, the robot pose index, and the plane index respectively. Let be one of the data points from the \aclrf which lies on the , be the normal unit vector of , and be the perpendicular distance from the origin of the robot coordinate system to . Since is located on the , it has to satisfy the following constraint,
(1) 
depends on the robot pose at pose index and the \aclrf extrinsic parameter , so (1) can be expanded,
(2) 
and are known approximately ( and ), can be computed from the robot’s joint angles at pose index , and is obtained from the laser. Let
(3) 
then
(4) 
The only unknown in (4) is which has 12 elements , where u and v denote the column and the row index of the matrix. Note that the fourth row of only consists of 0 and 1. Without loss of generality, let’s assume that the data points from the \aclrf lie on the XZ planes of the laser frame , so . If we expand (4) and rearrange such that the components of are stacked together as a vector , we have
(5) 
where
(6) 
and
(7) 
For each data point , we obtain such equation as in (5). With data points per pose and a total of robot poses, there are such equations. The equations can be stacked together to form the following matrix equation,
(8) 
where
(9) 
and
(10) 
Equation (8) can be solved by a linear leastsquare procedure to obtain . can then be computed from as follows:

The parameters and are required to be unit vectors, so they have to be normalized. They consistute the first and the third column of the matrix .

The parameters constitute the position component of the matrix (the 4th column).

The parameters can be calculated as the cross product of and .
has 12 parameters, but only 6 parameters are independent. To reduce the redundancy in the subsequent steps, the rotation part of is represented by the axisangle representation , while the position part is represented by .
IiiB Optimizing both the \aclrf Extrinsic Parameters and Robot’s Kinematic Parameters
In the second step, the data from all the three planes are used to optimize the extrinsic parameters of the \aclrf, the robot’s kinematic parameters and the plane parameters. The objective function is described as follows:
(11) 
The parameters consist of the following:

Robot’s kinematic parameters. We use DH parameters [14] to represent the robot’s kinematics, so there are 24 DH parameters for a 6\acdof robot arm.

\ac
lrf’s extrinsic parameters. As mentioned in the previous section, we use the axis angle representation for the rotation part , and for the position part.

Plane parameters. Each plane can be described by a unit vector normal to the plane and its perpendicular distance from the robot base’s coordinate system origin , so there are 12 parameters for 3 planes.
In total, there are 43 parameters to be optimized by minimizing the objective function . To do that, the number of data points have to exceed the number of parameters. The optimization problem is then solved using a LevenbergMarquardt nonlinear optimizer [15]. The objective function uses the geometric constraints that all data points from the \aclrf have to fall on the respective plane. Zhuang et al. [7] prove that a calibration process with such constraints is equivalent to the calibration of a robot using endpoint measurement in unconstrained calibration.
For the unit vector parameters ( and ), the following constraints are added to the optimization solver:
(12) 
(13) 
Further analysis on the observability of the parameters will be presented in the next section.
IiiC Identifiability of the calibration parameters
Depending on the chosen robot calibration poses and the robot’s kinematic model, some of the calibration parameters might not be observable due to the linear dependency among the parameters. This is a critical problem in calibration, as it will result in some of the parameters assuming erratic values which gives us unstable calibration result. To prevent that, we have to first analyse which calibration parameters are identifiable and which are not.
Following the approach in [8] and [16], SVD is applied on the identification Jacobian matrix . can be computed as follows. Let be the geometric constraint equation on the data point at the robot pose and on the plane k,
(14) 
Then can be computed by differentiating (14) for all the data points at the robot poses and for all the planes , then stack them together as a matrix,
(15) 
We can then apply SVD to the matrix ,
(16) 
Note that for this identification step, the parameters are excluded from the parameters vector , since those four parameters are obtained as linear combinations of other parameters (Equation (12) and (13)). That leaves us with 434 = 39 parameters in , which correlates to the 39 singular values in . The number of zero singular values in is then equal to the number of unidentifiable parameters in the calibration procedure. For a given zerovalue singular value , the rth column vector of the matrix is the linear combination of the parameters which cannot be identified independently.
i  

1  0.0  0.0  345.0  
2  90.0  0.0   90.0  0.0 
3  0.0  305.0  + 90.0  0.0 
4  90.0  10.0  300.0  
5  90.0  0.0  0.0  
6  90.0  0.0  70.0 
In this paper, we use a Denso VS060 6\acdof industrial manipulator with its DH parameters presented in Table II. The \aclrf’s frame is defined such that the rotation part , and the position part . Applying the identifiability analysis to the system, we found that there are 7 sets of linearly dependent parameters out of the 39 parameters.

The parameters (the translation along the zaxis of the 6th link frame on the flange) and (the z coordinate of the \aclrf frame) are linearly dependent. Physically this means that if we shift the origin of the 6th link’s frame in its z direction by changing , we can compensate by shifting the origin of the \aclrf frame in the opposite direction by changing .

The parameters and are linearly dependent. These correspond to the rotation of the 6th link’s frame and the rotation of the \aclrf’s frame around the same zaxis.

The parameters and are linearly dependent. These correspond to the shift in the zaxis direction of the 2nd and 3rd link’s frames respectively, which are along the same direction.

Lastly, we have four sets of linearly dependent parameters due to the linear combinations of the first link’s DH parameters and the three calibration planes’ parameters. Physically, this relates to the fact that we can shift the robot’s base frame freely by changing the value of , and the plane parameters will adjust according to the new base location. In other words, the base coordinate is not constrained (floating base).
For each set of the linearly dependent parameters, we can assign a fix value to one of the parameters. In this case, we fix the value of the parameters [] to their initial model’s values.
These results apply to most existing 6\acdof industrial robots whose structures are the same as that of our Denso robot.
Iv Simulation
We verify SCALAR through simulation of the calibration procedure. The simulation is conducted by using Robot Operating System and Gazebo where the robot model, the \aclrf, and the three planes can be simulated. As shown in Fig. 2, three perpendicular planes are located around the robot, and the robot is moved such that the \aclrf ray intersects each plane. Simulated data from the \aclrf can be obtained and Gaussian noise with zero mean and standar deviation can be added to the data. The data is then used as input to the calibration procedure.
After the calibration procedure, the robot is moved to 10,000 random poses to evaluate the accuracy of the calibrated parameters. Let and be the true and calibrated pose of the \aclrf’s frame w.r.t. the robot base frame at the robot pose index , respectively, then the error of the calibrated model can be computed as follows.
(17) 
Let be the element of with the subscript and refer to the row and column index, then the position error at the robot pose index , , can be computed by
(18) 
Let be the rotation part of the homogeneous transformation matrix . can be represented by using an axisangle notation, . We use as the orientation error at the robot pose index . can be seen as the amount of rotation necessary to rotate the calibrated pose to the true pose. The errors and are then averaged over the 10,000 random poses.




The simulated robot’s model is considered as having the true kinematics and extrinsic parameters, and an initial model is generated by introducing random Gaussian errors to the true parameters within the range of 2mm and 1 for the linear and angular parameters, respectively. Note that the initial model’s errors are intentionally set to be large to illustrate the robustness of the calibration method. The average position and orientation errors of the initial model as compared to the true model are 14.6mm and 4.05, while the maximum errors are 103.9mm and 5.95. We run the calibration procedure to improve the initial model with = 120, = 100 and = 0.1mm, and the resulting calibrated model has the average position and orientation errors reduced to 0.09mm and 0.02, while the maximum errors are reduced to 0.19mm and 0.035.
Next, we evaluate the effect of the measurement noise, the number of poses () and points (), and the plane parameters’ initial estimate on the calibration errors.
Iva Effect of the measurement noise
The accuracy of the calibration procedure greatly depends on the accuracy of the measurement system, which is affected by the noise on the data. In this section, Gaussian noises with zero means and varying standard deviations are added to the measurement data in the simulation, and its effect on the calibration errors is shown in Fig. (a)a. As decreases, the calibration errors decrease. At = 0.1mm, the average position and orientation errors are around 0.1mm and 0.02, respectively. For the subsequent sections, is set at 0.1mm.
IvB Effect of the number of calibration poses and the number of points
For each plane, the robot is moved to poses, so in total there are calibration poses. At each pose, we select data points from the \aclrf data. In this section we evaluate the effect of and to the calibration errors. In Fig. (b)b, it can be seen that as the number of poses increases, the error decreases until around =60, beyond which it does not change significantly. It can be concluded that 60 robot poses are sufficient to calibrate the robot model accurately. We also conduct similar analysis on (the data is not presented in this paper) and found that = 20 is sufficient for the calibration.
IvC Effect of the plane parameters’ initial estimate
One of the advantages of SCALAR is that the plane parameters do not need to be precisely known. Here we vary the plane parameters’ estimate to demonstrate the robustness of our method. The initial estimates of the positions and the normals of the planes are disturbed by up to 100mm and 30, as shown in Fig. (c)c and Fig. (d)d. From the figures, it can be seen that the calibration position and orientation errors are not affected by the errors in the plane parameters’ initial estimate. In fact, after calibration, the plane parameters in the calibrated model approach the true parameters within 0.1mm and 0.01.
V Conclusions
In this paper, we have proposed SCALAR, a method to simultaneously calibrate a 6\acdof robot’s kinematic parameters and a 2D \aclrf’s extrinsic parameters using only three flat planes, arranged perpendicularly towards each other around the robot. SCALAR is easier to implement than the previous methods in the literature as it does not require other expensive measurement systems or tedious setup. Through simulations, we have also verified that the method can reduce the average errors in position and orientation from (14.6mm, 4.05) to (0.09mm, 0.02), respectively. This is very useful for many industrial robotics applications that require great accuracy.
The next step after this is to implement SCALAR on the real system. Some challenges that may appear on the real system are the backlash of the robot’s transmission system, the elasticity of the joints, the roughness of the calibration planes and the noises of the \aclrf data which will reduce the calibration accuracy. Moreover, evaluating the calibration errors in the real system is not as easy as in simulation. We will present the calibration result on the real system in our future work.
Acknowledgment
This work was supported in part by NTUitive Gap Fund NGF201601028 and SMART Innovation Grant NG000074ENG.
References
 [1] F. SuárezRuiz, T. Santoso Lembono, and Q.C. Pham, “RoboTSP â A Fast Solution to the Robotic Task Sequencing Problem,” in International Conference on Robotics and Automation, 2018.
 [2] L. S. Ginani and J. M. S. T. Motta, “Theoretical and practical aspects of robot calibration with experimental verification,” Journal of the Brazilian Society of Mechanical Sciences and Engineering, vol. 33, no. 1, pp. 15–21, 2011.
 [3] S. H. Ye, Y. Wang, Y. J. Ren, and D. K. Li, “Robot calibration using iteration and differential kinematics,” Journal of Physics: Conference Series, vol. 48, no. 1, pp. 1–6, 2006.
 [4] A. Nubiola and I. A. Bonev, “Absolute calibration of an ABB IRB 1600 robot using a laser tracker,” Robotics and ComputerIntegrated Manufacturing, vol. 29, no. 1, pp. 236–245, 2013.
 [5] W. S. Newman, C. E. Birkhimer, R. J. Horning, and A. T. Wilkey, “Calibration of a Motoman P8 robot based on laser tracking,” in ProceedingsIEEE International Conference on Robotics and Automation, vol. 4, 2000, pp. 3597–3602.
 [6] M. Ikits and J. Hollerbach, “Kinematic calibration using a plane constraint,” Proceedings of the 1997 International Conference on Robotics and Automation, vol. 2, no. 4, pp. 3191–3196, 1997.
 [7] H. Zhuang, S. Motaghedi, and Z. Roth, “Robot calibration with planar constraints,” Proceedings 1999 IEEE International Conference on Robotics and Automation, vol. 1, pp. 805–810, 1999.
 [8] A. Joubair and I. A. Bonev, “Nonkinematic calibration of a sixaxis serial robot using planar constraints,” Precision Engineering, vol. 40, pp. 325–333, 2015.
 [9] Qilong Zhang and R. Pless, “Extrinsic calibration of a camera and laser range finder (improves camera calibration),” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), vol. 3, 2004, pp. 2301–2306.
 [10] Bouguet J. Y., “Camera calibration toolbox for matlab,” Tech. Rep., 2003.
 [11] R. Unnikrishnan and M. Hebert, “Fast Extrinsic Calibration of a Laser Rangefinder to a Camera,” Robotics, 2005.
 [12] G. Li, Y. Liu, L. Dong, X. Cai, and D. Zhou, “An algorithm for extrinsic parameters calibration of a camera and a laser range finder using line features,” IEEE International Conference on Intelligent Robots and Systems, pp. 3854–3859, 2007.
 [13] F. Vasconcelos, J. P. Barreto, and U. Nunes, “A minimal solution for the extrinsic calibration of a camera and a laserrangefinder,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 34, no. 11, pp. 2097–2107, 2012.
 [14] J. J. Craig, Introduction to robotics: mechanics and control, 3rd, Ed. AddisonWesley, 1986.
 [15] M. Newville, T. Stensitzki, D. B. Allen, and A. Ingargiola, “LMFIT: NonLinear LeastSquare Minimization and CurveFitting for Python,” 2014.
 [16] J. M. Hollerbach and C. W. Wampler, “The calibration index and taxonomy for robot kinematic calibration methods,” International Journal of Robotics Research, vol. 15, no. 6, pp. 573–591, 1996.