FEEDFORWARD CONTROL OF MULTISTAGE ASSEMBLY

PROCESSES USING PROGRAMMABLE TOOLING

L. Eduardo Izquierdo1, Jianjun Shi2, S. Jack Hu1, and Charles W. Wampler3

1
Department of Mechanical Engineering
The University of Michigan
Ann Arbor, Michigan
2
Department of Industrial and Operations Engineering
The University of Michigan
Ann Arbor, Michigan
3
Manufacturing Systems Lab
General Motors Corporation
Warren, Michigan

KEYWORDS variation in automobile body may lead to wind

noise and water leakage, thus variation should
Programmable Tooling, Multistage Assembly, be minimized whenever possible. The autobody
Process Control, Quality Improvement assembly process involves up to 150 parts
assembled in up to 100 stations, where the
sources of variation may come from any parts or
ABSTRACT assembly operations. Therefore, in such a
complex process, determining the deviations
The combination of feedforward control and and the appropriate correction for variation
programmable tooling has emerged as a reduction are always difficult and time
promising method to reduce product variation in consuming tasks.
multistage manufacturing systems. Feedforward
control allows compensation of deviations on a There are three approaches commonly used in
part-by-part basis using programmable tooling. variation reduction in manufacturing: robust
This paper addresses the problem of designing design, Statistical Process Control (SPC) and
an optimal feedforward control law that improves active error compensation. Robust design helps
quality. The controller design involves deviations to develop products and processes that are less
estimation (from in-line measurements), sensitive to part/process variation and
variation propagation modeling and analysis, disturbances. However, robustness does not
and process/parts constraints. Therefore, a mean complete rejection of errors; therefore,
control law is obtained using constrained production variation (parts and process
optimization. A case study is conducted on a variation) may still impact final product quality.
multistage assembly of a vehicle side frame to SPC has been successfully used to detect
illustrate the developed methodology. changes in product/process caused by mean
shift or variation changes, and also to identify
some predetermined root causes of variation.
1. INTRODUCTION However, SPC cannot be used to compensate
deviations on a part-by-part basis. This type of
Variation reduction is an important but compensation can only be achieved using active
challenging task in multistage manufacturing deviation compensation, which allows the
processes. As an example, dimensional

Transactions of NAMRI/SME 295 Volume 35, 2007

control of product/process deviations through to the joining by adjusting the position of the
corrections for each assembly. PTs. Following this approach deviations are
compensated, and quality is improved. Other
The enablers of active compensation are: benefits of using a feedforward control in
• Programmable Tooling (PTs): PTs allows assembly processes include reducing process
the automatic adjustment of locators and ramp-up time (time to market) and improving the
clamps used to hold parts. Because of the disturbance response time. Those advantages
high precision of the PTs and their capability not only improve quality, but also enhance
to perform part-to-part adjustments, they process responsiveness and reduce cost.
provide the capability to compensate
part/tooling deviations. One example of a PT One of the first attempts to use feedforward
is the Fanuc robot F-200iB (Fanuc, 2006), control on assembly was done by Svensson
which was the first robot introduced in (1985). With the help of a vision system, he
assembly to serve as a fixture carrier to modified the trajectory of a robot to achieve
allow the assembly of mixed models in the better fit of doors and windshields in car
same line. assembly. Similar applications were reported by
Sekine et al. (1991), Wu et al. (1994), and
• In-line dimensional measurement sensors: Khorzard et al. (1995), where different
The development of accurate non-contact techniques were used to determine the
sensors that can endure real process appropriate fitting of parts.
conditions has brought the possibility to
obtain in-line quality information on the The aforementioned feedforward control
assembly stages (Perceptron, 2006). strategies were related to the variation in one
• Stream-of-Variation (SoV) modeling tools: particular stage, without considering down
SoV tools allow modeling the variation stream processes. This single stage approach is
propagation process in multistage assembly only effective on reducing variation in multistage
processes (Hu 1997; Jin and Shi 1999; Shi assembly processes if the stage involved is the
2006). Those models can be used to last one, or the Key Product Characteristics
determine the impact that deviations and (KPC) of the product controlled in the stage will
control actions have on the final product not be affected by later processes. However, if
quality. none of those conditions hold, the single stage
control is not appropriate because it does not
Active error compensation in multistage consider the impact of deviations and control
assembly processes can be approached in two actions on the downstream part dimension.
ways: feedback control and feedforward control. Therefore, the control actions obtained for the
Feedback control implies that the control actions single stage scheme may not be optimal to
(corrections) are determined using downstream improve final product quality.
measurements usually obtained at the end of
the process or in certain intermediate stages. On A feedforward control in a multistage process
the other hand, feedforward control uses needs a model to determine the impact that
distributed sensors to determine deviations of control actions at one intermediate stage have
parts/process, and then apply control actions on the final product. Mantripragada and Whitney
before the joining takes place. In this way, (1999) proposed a multistage model and the use
feedforward control proactively compensates of optimal control theory to determine control
current deviations instead of reacting to past actions during the assembly. Using
deviations as feedback control does. measurements of the parts before assembly,
they were able to calculate the control actions
Product and process deviations in assembly that minimize the final product variation. They
can be understood as mean shifts and variance assumed that parts are the only source of
changes. Due to the usual absence of auto variation in the process. More recently,
correlation of the variation sources in multistage Djurdjanovic and Zhu (2005) proposed the use
assembly systems (Hu and Wu 1990), feedback feedback and feedforward control using the
control can only be used to compensate mean state space model to control deviations in
shifts, but not to reduce variability. Thus, multistage machining applications by modifying
feedforward control scheme is preferable in the position of the fixtures and tool path.
assembly processes to perform corrections prior However, those papers do not address the

Transactions of NAMRI/SME 296 Volume 35, 2007

feedforward control of the multistage assembly 2.1 SoV Model
process including parts/process requirements
and specific engineering constraints on the The model used here to describe the variation
derivation of the control actions. propagation in multistage assembly of rigid parts
is the state space model developed by Jin and
This paper presents a methodology to design
Shi (1999).
an optimal feedforward control that improves
product quality by considering process/parts
characteristics, multistage variation propagation,
and constraints in process and control actions
due to actuator characteristics, interference with
other components and other factors. Thus, the
determination of the control actions can be
formulated as a constrained optimization FIGURE 1. MULTISTAGE MANUFACTURING
problem, where the requirements to determine PROCESS.
the optimal actions are:
• Obtaining an expression of the final product
deviations (objective function) as a function Figure 1 presents a schematic of a multistage
of the control actions and the estimated assembly process. As the subassemblies are
parts deviations obtained from distributed moved from one stage to the next stage, they
measurements. sequentially accumulate errors (Shiu et al.
1996). This process can be modeled as,
• Defining the search space for the control
actions considering the PT’s constraints x k = A k −1x k −1 + Bk uk + w k (1)
and parts/processes characteristics.
yk = Ck x k + v k ; k = 1" N , (2)
• Determining the control actions that
minimize the effects that the estimated
deviations have on final product quality where Eq. (1) is the state equation, variable
without violating the constraints, by using a x k ∈ ℜn represents the state of the system in
suitable optimization method. stage k (part deviations from nominal). Variables
uk ∈ ℜ p and w k ∈ ℜn represent the fixture
The remainder of the paper is organized as deviations and the disturbances respectively.
follows: Section 2 presents the multistage
Matrix Α k ∈ ℜn×n stands for the reorientation
process model and the control estimation
problem. Section 3 addresses the part/process matrix, which relates the fixture layout of two
deviations estimation problem and the adjacent stages (k-1 and k). The effects of
development of the optimal feedforward control fixture deviations into the state of the system are
law in details. A case study is presented in determined by matrix Βk ∈ ℜn× p . The observation
Section 4, and the conclusions are given in equation, Eq. (2), is used to determine the
Section 5. deviations of the measurement points yk ∈ ℜm ,
which usually corresponds to the KPCs of the
product. Their deviations are obtained from the
2. FEEDFORWARD CONTROL OF
MULTISTAGE ASSEMBLY state using the observation matrix Ck ∈ ℜm×n
and adding the measurement noise v k ∈ ℜm .
This section formulates the optimal For details on how to derive each matrix please
feedforward control problem for multistage refer to Jin and Shi (1999), Ding et al. (2000),
assembly processes that includes part and and Shi (2006).
process constraints. First, the SoV model is
presented, which is used to determine the The state transition matrix Φk ,i describes the
impact that the control actions have on the final
product quality. Second, it is addressed the deviation transmission between stages i and k
determination of the control action as a (Ding et al. 2000), and it is calculated as
constrained optimization problem using Φk ,i ≡ A k −1A k − 2 " A i +1A i , k>i; Φi ,i ≡ I (I is the
estimated deviations. identity matrix). Then, Eq. (2) can be written as,

Transactions of NAMRI/SME 297 Volume 35, 2007

 N N determine the control action vector sk, which will
yN = Γ0 x 0 + ∑ Γk uk + ∑ Ψk w k +v N , (3) be applied by the PTs. Since the PT are not
k =1 k =1
perfect, there is going to be an error ek on the
where x 0 represents the deviation of the applied control actions. Finally the state of the
incoming parts, Γk = Ck ΦN ,k Bk , Γ0 = CN ΦN ,0 and system at stage k can obtained as,
Ψk = Ck ΦN ,k .
x k = A k −1x k −1 + B k (uk + s k + e k ) + w k . (4)

The deviations of the incoming parts, fixtures

Using the state space model is possible to
deviations, disturbances and noise are
determine the effect that the estimated
considered as random variables with mean of deviations and control actions have on the
zero and covariances of Σ x0 , Σu , Σw and Σv estimated final product deviations y N / k , given the
respectively. They are also considered to be information available in stage k (the derivation of
independent among them (e.g. Cov (uk , v k ) = 0 ) y N / k is presented in detail in section 3.2). By
and to be independent between different stages doing so, the control action determination can be
(e.g. Cov ( w i , w j ) = 0, ∀i ≠ j ). formulated as a constrained optimization
problem, where the objective function is the
In assembly, in-line measurements are usually weighted sum of the squares of y N / k as
obtained using OCMM sensors (Optical presented in Eq. (5).
Coordinate Measurement Machine), which
provide information on the displacement of the J = min y TN / k Q k y N / k
measurement points in 1D or 2D. The sk
, (5)
measurement points are usually selected to s.t. g( y TN / k , sk ) ≤ 0
coincide with the KPCs. Therefore, they
correspond to features that are important for the
functionality, cost and safety of the product. where matrix Q k ∈ ℜm×m is the weight matrix,
Qk is a diagonal and positive definite matrix. The
values of the weighting coefficients account for
2.2 Feedforward Control Problem the relative importance of the KPCs. The set of
Formulation constraints g(.,.) include the design and
manufacturing requirements for the location of
The feedforward control formulation is based the KPCs and stage/PT characteristics.
on compensating deviations of parts/
subassemblies and process before the joining
process takes place as presented in Fig. 3. 3. DETERMINATION OF THE CONTROL
ACTIONS

This section presents the procedures to

determine the optimal control actions using the
estimated parts deviations.

3.1 Deviation Estimation

Being in stage k, the system equations before

the control actions are applied can be described
FIGURE 3. PROCEDURE TO CORRECT as:
DEVIATIONS.

After the parts/subassemblies are mounted on x Bk = A k −1x k −1 + Bk uk + w k (6)

the fixtures, measurements are performed to
determine the deviations of the parts. Since the yBk = Ck x k + v Bk ; k = 1" N , (7)
measurements are corrupted with noise ( v Bk ),
then, the true state of the system can only be
estimated. Using the estimation, it is possible to

Transactions of NAMRI/SME 298 Volume 35, 2007

 B
where the super index stands for the condition J = min y TN / k Q k y N / k
sk
before applying the control. Using Eq. (7), it is
possible to estimate the state of the system xˆ k s.t. y N / k ∈ ⎡⎣LSLy N / k ,USLy N / k ⎤⎦
using the Weighted Least Squares (WLS) smin ≤ sk ≤ smax (11)
estimation method as,
⎧ s if sk ≥ Δ s
sk = ⎨ k
xˆ k = C†k yBk , (8) ⎩ 0 otherwise.

where, C†k is the pseudoinverse of matrix Ck , This general formulation includes the
−1 existence of constraints on the position of the
and it is calculated as C = (C R k Ck ) C R k .
†
k
T
k
T
k KPCs and the control actions. The first
Here matrix R k ∈ ℜm× m is a weighting coefficient constraint ensures that the final product KPCs
matrix, which accounts for differences in the are going to be within the Upper and Lower
importance and characteristics of the measured Specification Limits (USL and LSL). The second
points (KPC or non-KPC points, sensor noise constraint restricts the control actions to be
level, measurement point importance, etc.), and within the upper and lower PT actuation limits
it is a positive definite diagonal matrix. If matrix ( smin and smax ) that can be applied on each
R k contains on its diagonal the inverses of the part/subassembly. The control action limits
sensors noise variances, then, according to consider PTs workspace limitations and
Franklin et al. (1998), xˆ k is the best linear interferences with other stage components.
Finally, the third constraint is an or-type one,
unbiased estimator of x Bk . where there are two possibilities for sk: it is
either bigger than or equal to a threshold Δ s
( Δ s > 0 ), or it is zero. This type of threshold is
3.2 Control Action Determination
used to avoid obtaining control actions that
At stage k, it is possible to write down the cannot be performed by the PTs. The value of
effects that the different variation sources and the threshold can be obtained according to the
the control actions have on the final product accuracy of the PTs.
deviations y N as presented in Eq. (9).
Varying Δ s in Eq. (11) can be understood as
~ using different types of PTs. Therefore, such
y N = Ψk x Bk + Γ k ( s k + e k ) + v N . (9) study can lead to identify the appropriate PTs to
be used based on an effectiveness analysis.
The PTs error vector ek is assumed to be a
random variable with mean of zero and
covariance Σe , where the value of the
covariance depends on the precision
(repeatability) of the PTs utilized.

The expected deviations of the final product

measurements, given the information available
up to stage k, can be obtained by calculating the
expectation of Eq. (9) as,

y N / k = Ψk xˆ k + ΓCk sk . (10)

As presented in section 2.2 the control actions

are obtained based on the constrained
optimization of Eq. (5). Writing down the
constraints, the control problem can be FIGURE 4. PROCEDURE TO DETERMINE THE
CONTROL ACTION.
formulated as follows,

Transactions of NAMRI/SME 299 Volume 35, 2007

 Figure 4 presents the procedure proposed to the complete assembly is moved to stage four
determine the control actions, which is based on for final inspection. The locators sequence used
determining first the unconstrained optimal are: {(P1, P2), (P3, P4)}StageI , {(P1, P4), (P5,
solution of Eq. (11). If this solution does not P6)}StageII , {(P1, P6), (P7, P8)}StageIII and {(P1,
violate any constraints, then control action can P8)}StageIV . It is assumed that (i) all the required
be directly applied. If one or more constraints measurement points (marked in Figure 5b) are
are violated, then, the constrained optimization available at each station, (ii) PTs are used to
problem has to be solved. hold all the parts in stages I and III, and (iii)
there are not fixture errors and disturbances in
The unconstrained solution ( s∗kUnc ) of problem the measurement stage. The parts are assumed
(11) can be obtained by replacing Eq. (10) into to be rigid and the variation happens only in the
Eq. (11) and solving it as a WLS problem x-z plane
(similar to the one in Section 3.1). Following this
approach, the control action can be written down The control performance using the algorithm
in terms of the measurements before control as, derived in Sec. 3.2 are presented through (a)
calculating the Quality Index (QI) defined in Eq.
(14) as the reduction of the 2-norm of the final
s∗kUnc = −K k yBk , (12)
measurement standard deviation (σ) with and
without using control actions. (b) checking if
where the control gain matrix Kk is obtained as, some KPC’s deviations exceed the 2 mm/6σ
limit, which is the standard in the automobile
−1
K k = ⎡⎣(ΓCk )T Q k ΓCk ⎤⎦ (ΓCk )T Q k Ψk C†k . (13) industry.

To solve the constrained optimization problem, σ ywN/ ok Control − σ Control

yN / k
QI = × 100 . (14)
/ 2 2
we used the sequential quadratic programming σ w / o Control
yN / k
algorithm coded in Matlab (fmincon), due to its 2

flexibility to handle different types of constraints.

The first constraint in Eq. (11) may cause the

nonexistence of an optimal solution. This
happens when the incoming parts and
subassemblies at stage k are so severely
deviated from their nominal that it is impossible
to satisfy the first constraint. Therefore, there is FIGURE 5. SCHEMATICS OF AN SUV SIDE
no control action capable to adjust the KPCs to FRAME AND ITS SIMPLIFICATION (Ding et al.
make them be within their specification limits. If 2002).
that is the case, the control action should be set
to zero.
TABLE 1. EFFECT OF THE CONTROL AND THE
4. CASE STUDY CONSTRAINTS ON THE MEASUREMENT POINTS
QUALITY (Units: mm).
The case study simulates the assembly of a
Control w/o Control with
Sport Utility Vehicle (SUV) side frame (Fig. 5), W/o control
const. const.
Meas.
proposed by Ding et al. (2002). The side frame point Std in Std in Std in Std in Std in Std
is formed by four parts, which are assumed to x z x z x in z
be rigid and free to move in the x-z plane only (3 M1 0.16 0.06 0.06 0.04 0.07 0.04
dof per part).
M2 0.16 0.32 0.06 0.10 0.07 0.12
M3 0.34 0.24 0.11 0.11 0.13 0.13
The assembly is performed in three stages with
M4 0.34 0.27 0.17 0.16 0.13 0.17
final measurement taken in a final inspection
M5 0.24 0.22 0.23 0.15 0.23 0.16
stage. The assembly sequence is summarized
as follows: in the first stage the fender is M6 0.28 0.24 0.28 0.15 0.28 0.15

attached to the A-pillar, then the B-pillar is M7 0.22 0.25 0.08 0.10 0.12 0.12
added in the second stage, and in the third M8 0.22 0.07 0.08 0.04 0.12 0.05
stage the rear quarter is attached. Afterwards,

Transactions of NAMRI/SME 300 Volume 35, 2007

 The parameters used in the simulations are corrected by relocating the subassembly formed
Σ x0 = 0.04 · I, Σu = 0.0017 · I, Σw = 0.0001 · I, by parts 1, 2 and 3 in stage III (the subassembly
can not be stretched or compressed to correct
Σe = 0.0017 · I, and Σv = 0.0009 · I, where the
2
the errors). However, a significant portion of the
units are mm , and I stands for the identity deviations in the z direction (~30%) can be
matrix with appropriate dimensions. The USL corrected by relocating the subassembly in
and the LSL were set to 0.8 mm and -0.8 mm stage III.
respectively for all the KPCs. The values of smin
and smax were set to -5 mm and 5 mm In the case with control and with constraints,
respectively to all the PTs, and the value of the different scenarios are analyzed to study the
threshold Δ s was set to 0.1 mm for all the PTs. impact of the threshold Δ s has on the quality
The weighting coefficients matrices Qk and Rk improvement. This analysis may help to select
(k=1, 2 and 3) were set equal to the identity the appropriate PT for a given process. Table 3
matrix. The results reported are based on the and Fig. 7 present the results using different
simulation of 1500 assemblies. values of the threshold. The effect on the
improvement has a sigmoid shape, where for
Table 1 and Fig. 6 present the standard small thresholds, equivalent to using an
deviation of the KPCs in the x and z directions accurate PT, there is a small drop in the QI.
for the cases with and without control. In the However as the threshold increases (greater
with-control case, both with and without than 0.1 mm) the QI tends to decay
constraints scenarios are included. asymptotically to zero.

TABLE 3. EFFECT OF THE THRESHOLD Δs ON

THE QUALITY IMPROVEMENT.

Δs mm QI % Exceed 2 mm (6σ) ?
0 44.5 No
0.05 44.0 No
0.075 42.6 No
0.7 40.2 No
0.15 33.7 No
0.2 23.8 No
0.25 14.8 No
0.3 9.3 No
0.35 4.8 No
FIGURE 6. STANDARD DEVIATIONS OF THE 0.4 1.2 Yes
MEASUREMENTS POINTS (KPCs).

The effect of using control significantly

improves quality. The values of QI are 46.2%
and 41.2% for the unconstrained and
constrained control respectively. As can be
expected, the constrained control improvement
is less than the unconstrained one. However,
the drop is not too big, and it reflects a more
realistic process. By analyzing the figure, it is FIGURE 7. EFFECT OF THE THRESHOLD Δs ON
possible to observe that only the uncontrolled THE QUALITY IMPROVEMENT.
case exceeded the 2mm 6σ (σ=0.33 mm)
threshold. Due to the lack of actuators in stage
II, the effects of both controls do not significantly 5. CONCLUSIONS
improve the position of KPCs M5 and M6 in the x
direction. However, both controllers help in the z This paper proposes a new approach to
direction. The reason is that if part 3 is joined in improving product quality in multistage assembly
stage II in the wrong position, the deviations that processes by deviation compensation using
this part has in the x direction can not be feedforward control. The proposed method

Transactions of NAMRI/SME 301 Volume 35, 2007

makes uses of distributed sensing and Franklin, G., J.D. Powel, and M. Workman
programmable fixturing technologies in (1998). Digital Control of Dynamic Systems, 3rd
determining and correcting deviations on a part- ed. Addison Wesley.
by-part basis. The problem of determining the
optimal control action is formulated as a Hu, S. and S.M. Wu (1990). “In-Process 100%
constrained optimization by considering design Measurement and Statistical Analysis for an
specifications and actuator/process Automotive Body Assembly Process.”
characteristics. A method is proposed to obtain Transactions of NAMRI/SME, Vol. 28, pp. 317-
the solutions to the constrained problem by first 321.
solving the unconstrained one and then
searching the constraint boundaries to find a Hu, S.J. (1997). “Stream-of-Variation Theory in
global optimal solution. A case study that Automotive Body Assembly.” Annals of the
considers the assembly of a SUV side frame in CIRP, Vol. 46/1, pp. 1-6.
three stages is presented considering the
existence of PTs in only two stages. The results Jin, J. and J. Shi (1999). “State Space Modeling
proved that feedforward control reduces the of Sheet Metal Assembly for Dimensional
variation of the final product KPCs by more than Control.” ASME J. of Manuf. Science and
40%. The effect of PTs accuracy on the resulting Engineering, Vol. 121(4), pp. 756-762.
quality improvement is also analyzed. From this
analysis it can be concluded that for high PT Khorzard, D., J. Shi, S.J. Hu, J. Ni, E. Zussman,
accuracy the effect is almost constant. However, and G. Seliger (1995). “Optimization of Multiple
as the PT accuracy diminishes, there is a Panel Fitting in Automobile Assembly.”
significant decrease in the amount of variation Transactions of NAMRI/SME, Vol. 23, pp. 241-
that can be reduced. 246.

Mantripragada, R. and D.E. Whitney (1999).

ACKNOWLEDGMENTS “Modeling and Controlling Variation Propagation
in Mechanical Assemblies Using State
This research was supported in part by the Transition Models.” IEEE Trans. on Robotics
General Motors Collaborative Research Lab in and Automation, Vol. 115, pp. 124-140.
Advance Vehicle Manufacturing at The
University of Michigan. Perceptron Inc. (2006). www.perceptron.com

Sekine, Y., S. Koyama, and H. Imazu (1991).

REFERENCES “Nissan’s New Production System: Intelligent
Body Assembly System.” SAE Conference,
Ding, Y., D. Ceglarek, and J. Shi (2000). Number 910816, Detroit, MI.
"Modeling and Diagnosis of Multi-station
Manufacturing Processes: State Space Model.” Shiu, B., D. Ceglarek, and J. Shi (1996). “Multi-
Proc. of the 2000 Japan/USA Symposium on Station Sheet Metal Assembly Modeling and
Flexible Automation, Ann Arbor MI. Diagnostic.” Transactions of NAMRI/SME, Vol.
24, pp. 199-204.
Ding, Y., D. Ceglarek, and J. Shi (2002). “Fault
Diagnosis of Multistage Manufacturing Shi, J. (2006). Stream of Variation Modeling and
Processes by Using State Space Approach.” Analysis for Multistage Manufacturing
ASME J. of Manuf. Science and Engineering, Processes, CRC Press.
Vol. 124, pp 313-322.
Svensson, R. (1985). ”Car Body Assembly with
Djurdjanovic, D. and J. Zhu (2005). “Stream of ASAE 3D-vision.” Proceedings 15th Int.
Variation Based Error Compensation Strategy in Symposium on Industrial Robots, Tokyo, Japan.
Multistation Manufacturing Processes.” ASME
Int. Mech. Eng. Congress, Orlando, FL. Wu, S., S.J. Hu, and S.M. Wu (1994). “Optimal
Door Fitting with Systematic Fixture
FANUC Robotics Inc. (2006). Adjustment.” Int. Journal of Flexible
www.fanucrobotics.com. Manufacturing Systems, Vol. 6, pp. 99-121.

Transactions of NAMRI/SME 302 Volume 35, 2007