Received 22 March 2016; accepted 23 May 2016; published 26 May 2016
The fundamental concepts and characteristics of Petri Nets (PNs) made them a significant tool for modeling and analyzing asynchronous systems with concurrent and parallel activities. A Petri net is both a graphical and an analytical modeling tool. Carl Adam Petri from Bonn University, Germany, developed it as a special class of generalized graphs or nets. The chief attraction of a Petri net is the way in which it identifies the basic aspects of various systems, conceptually, through a graphical representation, and mathematically, eventually supported by formal programming languages. As a graphical tool, PNs are a visual-communication aid similar to flow charts, block diagrams, and networks. In addition, tokens (or markers) in these nets can simulate the dynamic activities of systems. As a mathematical tool, it allows setting up algebraic equations for formally analyzing the mathematical models governing the systems’ behavior, similarly to other approaches to formal analysis, e.g. occurrences nets and reachability trees. Using Petri nets to model and analyze asynchronous systems with concurrent and parallel activities is simple and straightforward making them a valuable technique for analyzing such complex real time systems. They are capable of modeling material and information flow phenomena, detecting the existence of deadlock states and inconsistency and are well suited for the study of Discrete Event Systems (DES). As such, they have been applied in manufacturing and logistics, hardware design and business processes, as well as in distributed databases and DES simulation  -  .
There is no time measure in the original PNs. Hence, their modeling capability was limited and models often became excessively large, because all data manipulation had to be represented directly into the net structure. To overcome these drawbacks, over the years PNs have been extended in many directions including, time, data and hierarchy modeling. To adapt Petri nets to the study of systems’ performances, Timed Petri Nets (TPNs), allowing for delays in the system, were introduced. There are several possibilities for introducing time in PNs, thus transforming them into TPNs, among them timed transitions and timed places (see e.g.   for detailed presentations).
The objective of this paper is to show different ways to utilize this modeling frame. The paper reviews several published papers where the framework was used in different circumstances such as, estimating expected utilization of resources at steady state in open queueing networks, verifying computerized simulations and batch planning in textile industries.
The paper is organized as follows. Section 2, entitled “From Petri Nets to Timed Petri Nets”, presents basic PN modeling and describes a technique for modeling TPNs, which associates time with places. Section 3 entitled “Decomposing TPNs” reviews a method for decomposing Timed Petri Nets of Open Queuing Networks. Section 4 describes the decomposition application for verifying computerized simulation. Section 5 entitled “Application of TPNs to textile processing” describes an additional published example. Section 6 concludes the paper.
2. From Petri Nets to Timed Petri Nets
The guiding principles of net theories concern states/conditions and changes/events, two fundamentally intertwined concepts. In a Petri Net the relationships between these two concepts are both graphically and analytically described. A PN consists of two parts: 1) its structure and 2) its marking representing the system overall state on the defined structure. The net structure describes its possible states (places) and events (transitions). When one or more markers (tokens) reside in a place, those tokens mark it and thus determine its local state. The state of the entire system is determined by the net marking, which is the distribution of markers to the net places. Tokens move according to firing rules. Firing of transitions and flows of tokens, model the dynamic behavior of the net. Let us first present PNs as a graphical tool and then as an analytical tool.
2.1. PNs as a Graphical Tool
A Petri net is a graph, where:
T is a set of transitions, portrayed by bars, representing instantaneous or primitive events.
P is a set of places, portrayed by circles, representing conditions.
I(T) denotes the set of all input places of transition T. Similarly, O(T) denotes the set of all output places of transition T.
A is a set of directed arcs (T*P) U (P*T) derived through an incidence matrix (see Section 2.2). A “marking”, M of a Petri net is a distribution of tokens (or markers) to the places of a Petri net.
A transition can “fire” (meaning the respective event occurs) if there is at least one marker in each of its input places (the necessary conditions for its occurrence). Following the firing of a transition, the number of tokens in the corresponding input places is depleted by one.
Figure 1 illustrates a simple example of an ordinary PN, PN1.
Its structure is as follows:
Figure 1. An ordinary Petri Net, PN1.
We see that the original marking of the net contains one token in place P0. The firing of transition T1 will insert one token in P1. As P0 and P1 are the two input places of transition T3, they will enable its firing. Firing of transition T3 will remove the tokens in P0 and P1 and insert a token in P3. Then, transition T4 will fire and will thus remove the token in P3 and insert a token in P0. This process describes a firing sequence. Another firing sequence will take place if following the current initial state of the system, transition T2 fires introducing a token in place P2. The respective tokens in places P0 and P2 will enable the firing of transition T5 that will remove the tokens in P0 and P2 and will add a token to place P4. This will be followed by the firing of transition T6 removing the token in P4 and inserting a token in P0.
Assume now that the initial state of the net contains not one token but three tokens, one in each of the three places P0, P1 and P2. This situation represents a conflict as both transitions T3 and T5 are enabled. Which transition should fire, T3 or T5? A rule for solving the conflict should be introduced, or the system should be redesigned to eliminate the conflict.
2.2. PNs as an Analytical Tool
Formal analysis of the model can be carried out when a PN is described by linear algebraic equations.
The structure of a PN is an integer matrix C, called an incidence or flow matrix. C is an m × n matrix whose m rows correspond to the places of the net and whose n columns correspond to the transitions. In ordinary PNs an element of this matrix is defined as follows:
It means that upon firing transition one token is removed from each input place and one token is added to each output place. The incidence matrix C1 of the PN1 in Figure 1 is given below.
By convention, the flow of tokens through transition is denoted by, j = 1, 2, 3, 4, 5, 6.
2.3. Time modeling of PNs
According to the version, introduced in  and adopted by others e.g.  , Timed Petri Nets (TPNs) associate a positive number s(T) to each transition T. In  a different approach to model time in Petri nets was developed, relying on the classical Petri nets approach to a non-primitive event (with time greater than zero) which prescribes its decomposition into two transitions representing two primitive events. These are respectively a commencement and a termination of the original non-primitive event. The author associates time duration (delays) with the places describing the occurrence of the non-primitive events. Other scientists also adopted this approach. The delay is modelled either as a deterministic or as a stochastic variable. However, if a stochastic model is used, calculations only consider the expected values.
Our approach uses timed places to model processing of multiple part classes within a mix, repairs following machine breakdowns and eventually set-ups ensuing a change in mix. Since most of these variables are stochastic by nature, the actual values are their respective expected values. For instance, the input of parts to the system and the machine breakdowns are respectively the expected time between arrivals and the expected time between failures. It is worth mentioning that the current method is not Stochastic Petri Nets, which deal with non-primi- tive events whose duration is assumed to be exponentially distributed  . A Timed Petri Net (TPN) graph looks exactly like a PN graph, except that tokens in timed places become available only after zi time units (zi > 0). This delay is independent of the token’s arrival instant. A TPN is described by its incidence matrix, which is exactly like that of its equivalent PN, to which a diagonal matrix Z, representing time delays, is added.
Let us now turn our attention back to Figure 1. Examining PN1 we observe that places P3 or P4 are as described above, meaning that a token in P3 or in P4 may well represent concurrent, time elapsing activities. We may now interpret PN1 as a PN describing a resource that provides service to two types of customers. The firing of transitions T1 or T2, respectively represents input of type 1 or type 2 customers to the system. A token in place P1 or P2, respectively, represents a type 1 or type 2 customers waiting to be served by the resource. A token in place P0 shows that the resource is available. Firing of transition T3 or T5, respectively, expresses the beginning of servicing a type 1 or type 2 customers. Firing of transition T4 or T6 expresses the end of service. This occurs after z3 or z4, time units have respectively elapsed. Then, one token returns to P0, meaning the resource is again available.
To convert PN1, illustrated in Figure 1, into TPN1, we use incidence matrix C1 and delay matrix Z1. Z1 is defined as a square diagonal matrix with elements Z1 = zi, i = 0, 1, 2, 3, 4 for i = j and zero otherwise. The timed places are z3 and z4 respectively representing service time of type 1 and type 2.
The graphical representation of TPN1 is also Figure 1 to which z3 and z4 are respectively added (eventually in parentheses) below places 3 and 4.
The flow of tokens through transition, j = 1, 2, 3, 4, 5, 6 is vector I.
3. Decomposing Timed Petri Nets
The author of  provides a detailed illustration and proof for decomposing a TPN model of n-part type open queuing network with M workstations, into M independent TPNs. A part (customer) type is a composite of specific jobs to be performed by workstations (resources) in a predetermined order. In open queueing networks, customers enter and leave the network in an unconstrained manner.
The essence of the approach consists in proving that, under steady state conditions as defined in  and applied to open systems, the arrival flow of any such job to its assigned workstation equals the input flow to the system of its parent part type.
Figure 2 graphically illustrates an example of an open queueing network at steady state, as presented in  . It consists of M = 3 workstations and n = 4 part types. Parts of type j, (exogenous customers comprising the mix served by the queueing network) arrive at the system with given flow rate, j = 1, 2, 3, 4. Processing of a part is composed of a number of operations (jobs) to be performed in a given order by previously assigned stations (deterministic routing). Upon arrival, each part is routed to its initial station. Part types 1 and 2 enter the system at workstation 1, while part types 3 and 4 enter the system at workstation 2. Part types 1 and 4 exit the system at workstation 3, while part types 2 and 3 exit at workstation 2. The flow labeling in Figure 2 indicates the preservation of the exogenous input flows throughout the system. According to the conditions detailed in Section 3.1, this means the queueing network operates at steady state.
3.1. TPN Decomposition Rules
An important strategy for managing complex systems is breaking them down into appropriate subsystems and providing means for integrating them.
A TPN model of an n-part type open queueing network with M workstations can be decomposed into M independent TPNs, provided steady state conditions are fulfilled at each work station m,
As stated in  and repeated in  , given a timed Petri net by its incidence matrix C, the net functions at its steady-state rate for a given flow vector iff satisfies the equations:
is a generator of the set of solutions of Equation (1)
Z is the delay matrix (on places) defined as a square diagonal matrix with elements for i = j and zero otherwise (see detailed matrix in Section 2.3).
Q(0) is a vector representing the initial marking of the net.
(C+) is obtained by replacing the “−1” elements in the incidence matrix C by zeros.
Intuitively Equations (1) and (2) can be respectively associated with conservation of flow and conservation of tokens.
3.2. Application of the Decomposition Rules
Let us now apply the decomposition rules to the queueing network in Figure 2. The two steady state conditions,
Figure 2. An open queueing network.
Equations (1) and (2), should be fulfilled by each of the three stations.
Here, we shall only analyze workstation 1. This station processes two part types, type 1 and type 2, thus matching TPN1 in Figure 1. We may notice that station 1 is the entering station of these two part types to the open queueing network. Hence, their arrival flow at this station is their respective exogenous flow to the system, and. According to Figure 2 when part type 1 leave station 1, they proceed to station 3, while part type 2 proceed to station 2.
Applying Equation (1) to the incidence matrix of Figure 1, (), we obtain, respectively, from rows 2 and 3: and. Rows 3 and 5, respectively, yield and. The above results comply with the equation in row 1. Combining the equations, we obtain,
It is thus seen that under steady state conditions the entering (exogenous) flow of type 1 parts to the system, equals I4, the firing rate of transition T4 representing the exit flow of type 1 parts from station 1 and also their entering flow to station 3. Similarly, at steady state the entering (exogenous) flow of type 2 parts to the system, equals I6, the exit flow of type 2 parts from station 1, which is also their entering flow to station 2.
In this example, one generator of the set of solutions of Equation (1) is J1,
The “1” values in J1 represent all the possible steady states of workstation 1: P0-idle, P3-processing part type 1, P4-processing part type 2.
Substituting Equation (4) in Equation (2) and using Equation (3) leads to:
Q0, Q3, Q4 represent the token content of places P0, P3, P4 for any firing sequence. Hence, the Left Hand Side (LHS) of Equation (5) represents the sum of tokens over the steady states of station 1. Since these three states represent all the mutually exclusive states of workstation 1, LHS equals 1. The three summation components of the Right Hand Side (RHS) represent the respective proportion of time station 1 spends in each state under steady-state conditions or, in probabilistic terms, it represents the steady state distribution over the three states. After substitution, Equation (5) becomes:
and are the respective contributions of part type 1 and part type 2 to the utilization of station 1. As and represent the total entering flow to station 1, intuitively, z0 is the station expected idle time between any two consecutive arrivals of parts. Any feasible solution has to satisfy z0 > 0; leading to:
We may now use the above results to estimate the utilization of a station.
3.3. Estimating a Station Utilization Using TPN
Based on the application of the decomposition rules presented in the previous section, the procedure for computing the utilization of a station is summarized as follows:
The total expected steady-state utilization rm of station m, is the sum of its expected utilization by each class of customers. These comprise the n-part types as determined by a given deterministic routing and eventual machine failures or other disturbances including eventual set-ups.
where and represent the respective contributions of, the input flow of part types j, and that of If, the flow of disturbances, to the m station utilization, ,
is the expected operation duration of part j at station m.
is the expected duration of a disturbance.
Provided the constraint rm < 1 is satisfied for each, and all buffers are of unlimited size, Equation (8) is a valid expression of the TPN based expected steady state utilization rm of any processing station m, in the system.
In the next section, we shall use Equation (8) for verifying simulation models results of computerized open queueing network at steady state.
4. Using Timed Petri Nets as a Simulation Verification Tool
This section is a simplified version of a paper describing the usage of decomposed Timed Petri Nets as an analytical approach for verification of simulation models  .
An important and difficult task in any simulation project is the validation and verification of the simulation model. Validation is commonly defined as the process intended to substantiate that the conceptual simulation model, within its domain of applicability, is an accurate representation of the real world system it describes. The verification process seeks to determine that the computerized simulation model was built right, and operates as intended. There is a rich literature on these topics, see e.g.  -  .
From an output analysis perspective, system simulation is typically classified into terminating or steady state. While terminating simulations aim to estimate system parameters for periods well defined in terms of starting and ending conditions, steady state simulations aim to estimate system parameters as limits obtained for infinite simulation time. Before reaching steady state conditions when the system parameters attain stable values, the system behavior undergoes transition periods, reflecting the initial system operating conditions (see e.g.  ). In the analysis of such systems, an important objective is to avoid bias in parameter estimation, eventually introduced by data collected during transition periods. One among the many approaches to simulation verification is comparing simulation outputs with analytical results  . In many simulation studies dealing with manufacturing systems or computers and communication systems, the networks of queues represents the reality modelled by the system simulation. Hence, for simulation verification of such systems analytical approaches describing queueing networks may be appropriate.
The TPN approach puts into evidence the flow realization along a variety of network paths and may consider different network structures. Such a decomposed TPN model of a resource enables a quick calculation of its expected utilization at steady state thus providing a stick yard against which to verify the long run simulation based utilization estimate of the same resource. In other words, it represents a verification method.
4.1. A Numerical Example
There are three stations (M = 3) and four part types (n = 4).
1) The inter-arrival time for part j, j = 1, 2, 3, 4 is Erlang distributed with mean Aj (time units)
The reciprocal of the inter-arrival mean serves to estimate the average input flow of parts. Hence the input flows (in parts per time unit, pptu) are:
2) The service times (in time units) of parts j, j = 1, 2, 3, 4 at stations m = 1 and m = 3 are normally distributed, while the service time at station m = 2 is Uniformly distributed, each with their respective mean, , m = 1, 2, 3.
4.1.1. The Simulation Verification Procedure
We shall first use Equation (9) to calculate, for j = 1, 2, 3, 4, the contribution of each part type to the station’s utilization, and then Equation (8) to calculate rm, the total expected utilization of station m (m = 1, 2, 3).
4.1.2. Simulation Results and Conclusions
The simulation project used Siman V software and the Arena framework. As the published example investigated the effect of several factors such as the planned utilization of the stations and the queue discipline, forty-eight simulation runs were carried out. Here we do not investigate such effects. We solely compared the TPN calculated estimates with the simulation average results after we deleted the transient period (about 400 minutes).The simulation average estimating the station utilizations were close to those calculated using TPN and their differences were no higher than 2%.
The relatively low differences obtained between the numerical results of the two methods provide evidence that the TPN technique is able to verify simulation models. Observed discrepancies between the TPN based and the results of the computerized simulation model can thus detect eventual flaws in the simulation model.
5. Application of TPNs to Textile Processing
We shall briefly present an additional example dealing with the application of TPNs to textile batch planning  .
The need for advanced methods to model and characterize the performance of planning and scheduling is especially apparent in the textile dyeing and finishing operations. There is a huge number of process combinations that a fiber, yarn or fabric may go through on its way to the final product. Well organized planning and control of the entering materials, flowing through the intricate manufacturing operations to reach their final state is critical for meeting customers due dates and thus provide on time delivery. The TPN decomposition approach has been developed and applied to calculate long run resource utilization for manufacturing systems of discrete items. The textile manufacturing conditions are different. Workstations (machines) are of two types: batch processing (especially for dyeing) and continuous processing (for finishing operations such as squeezing or drying). The discrete items entering the system are orders. Orders vary in size, fabric type, fiber blend, color and the kinds of required operations. For batches, an order “size” is measured in weight units while for continuous processing it is measured in length units. An order is characterized by its specific processing requirements and by its specific attributes such as fiber blend and Machine Direction (MD) weight (kg/m). Incoming orders are stored and wait for processing. There are two main operation categories: equipment oriented operations and work force oriented operations such as lot transporting, loading and unloading operations. As the relative duration of the latter is insignificant, here these operations are disregarded.
To apply the TPN decomposition approach for calculating the expected (long run) utilization at steady state of the equipment used in the textile processes, we shall adapt the basic assumptions in  , detailed in Section 3, to a textile dyeing workstation.
5.1. Manufacturing Assumptions of the Textile Dyeing Operations
This section is an extension of Section 3.3, which estimates a station utilization of discrete systems.
Orders of type j, arrive at the system with given flow of weights, , (Kg/hour). Processing of an order is composed of a number of operations (jobs) to be performed in a given order by previously assigned stations (deterministic routing). Upon arrival, each order is routed to its initial station, where it is stored and waits for processing. Typically, the initial station is a dyeing workstation, the example here.
5.1.1. Dyeing Operations
- The dyeing operations are performed in batches, not exceeding the given manufacturing capacity of the dyeing machines. We assume that all dyeing machines have the same manufacturing capacity, G (kg/batch).
- A batch of type j is the maximal accumulated weight of type j orders in store, since the starting of the last batch, not exceeding the machine capacity, G.
- The average input flow of type j weights, , is transformed into, the average input flow of type j batches, and is calculated as, (batches/hour).
- A dyeing machine, m, has a batch processing duration, , for processing operation k of type j order, , ,.
The duration of this operation is dependent on the required “basic dyeing operation”, such as washing, bleaching or dyeing to a specific shade as well as on the specific order attribute in terms of “fiber blend”.
5.1.2. Evaluating the Utilization of the Dyeing Machines at Steady State Using TPN
Let be the contribution of, the input flow of type j batches, for possible multiple visits to station m.
is the given routing of batches j to machine m for performing operation k. It may assume a “0” or “1” value.
is the number of operations of type j batches at machine m.
is the expected processing duration of operation k of batch j at station m.
The total expected steady-state utilization rm of station m, is the sum of its expected utilization by each class of customers to be served by the machine. Here these are the n different types of batches as determined by the given deterministic routing and eventually failures or other disturbances as an additional class of customers. As the batches replace the discrete part types in Section 3.3, the following equation is identical with Equation (8),
Provided the constraint is satisfied for each, and all buffers are of unlimited size, the above equation is a valid expression of the TPN based expected steady state utilization, of any dyeing station m, in the system.
5.2. An Example
The example is based on a real textile dyeing and finishing mill which performs about 25 - 30 specific processes. We have selected eight among these processes, whose accumulated input flow in kg/hour represents about 77% of the total input flow to the system.
5.2.1. Expected Utilization of the Dyeing Machines
The average overall input order rate is 2 orders per hour and the mean weight per order is 170 kg. Accordingly, the total input flow of weights of the selected processes to the first processing machine, which is a dyeing machine (DM) is 262 kg/hour. There are 10 DM machines assigned to the selected processes and the manufacturing capacity of each is G = 250 Kg.
A type j order arriving at a dyeing machine where its first operation (k = 1) will be performed, is characterized by its “basic dyeing operation” (washing, bleaching, pale shade or dark shade, i.e. four possibilities) and by its specific fiber blend (cotton/viscose, cotton/lycra, cotton/polyester, polyamide or acetate, i.e. five possibilities). The batch processing duration at any dyeing machine is dependent on the above combinations.
A type j batch comprises fabrics requiring the same type of basic operation and the same dyeing duration. Following this rule, eight types of batches were defined, each with a given weight flow (kg/hour) and given batch duration, , , , see Table 1. Assuming that only one batch type is assigned to a machine, the results in Table 1 show that for batches of type j, j = 1, 2, 3, 4, 5, 6 and 7, one machine is needed for each, whereas for j = 8, four machines are needed. The total number of machines needed becomes 11, while the total number of machines available for these processes is only 10. Hence, we have to assign more than one type of batch to a machine.
To switch from one basic dyeing operation to another requires a complete cleaning operation of the machine whose duration is 4 hours. To switch between fiber blends within the same basic dyeing operation requires partial cleaning whose duration is 2 hours. It is to be noted that j = 2, 3 need the same basic dyeing operation and so do j = 4, 5, 6 and j = 7, 8.
Considering these set-ups as additional types of customers, we shall calculate their contribution to the machine utilizations using Equation (8) as well.
To obtain a feasible planning solution to our problem (less than 11 machines) we have to select at least two batch types to be processed by the same machine while keeping the expected utilization of each machine below 1. Batch types 5 and 6, respectively allocated to machines 5 and 6, contribute to a low utilization of these two machines and they need only partial cleaning for switching, 2 hours.
Hence we may select batch types j = 5 and j = 6, to be performed by machine 5. The set-up duration for cleaning, , has to be added to the respective processing duration of both batch type 5 and batch type 6.
Accordingly, the expected utilization of machine 5 is calculated as follows:
In the same manner, we may further reduce the number of operating machines.
In this example, we extended the TPN decomposition approach as applied to manufacturing systems of discrete items to accommodate conditions existing in textile processing systems. Under these conditions, the input to the system is through flows of orders. These were transformed into flows of batches, which now replace the flows of discrete items. The procedure also assists in planning the allocation of different types of batches to machines.
Table 1. Machine utilizations by different batch types.
The reviewed published papers showed that Petri nets are a versatile modeling structure. Section 2 of this paper presented the important extension of Petri nets to Timed Petri Nets (TPNs), which are able to time-base evaluate system performances. The section described in some details a technique for introducing time in Petri nets, by associating it with places. Section 3 showed how timed petri nets of open queuing networks were decomposed. Through decomposition, complex networks systems became easier to manage. Long run expected utilization of workstations at steady state, could be quickly calculated using decomposed TPNs. This enabled to use TPNs as a simulation verification method, as described in Section 4. The reviewed paper compared simulation results with TPN easily calculated estimates of stations utilizations. The low differences between the two techniques showed the advantages of TPNs as a simulation verification method. In the last presented paper, the decomposition approach to TPNs was applied to textile industry. Orders of different types are customary input to textile process- ing. The described procedure transformed the input flows of different types of orders into input flows of different types of batches, enabling to estimate the long run machines utilizations. The procedure also assisted in planning the allocation of different types of batches to operating machines.
 Billington, J., Diaz, M. and Rozenberg, G. (1999) Application of Petri Nets to Communication Networks. Springer-Verlag, Berlin.
 Yakovlev, A., Gomes, L. and Lavagno, L. (2000) Hardware Design and Petri Nets. Kluwer, Boston.
 Narahari, Y. and Viswanadham, N. (1985) A Petri Net Approach to the Modelling and Analysis of Flexible Manufacturing System. Annals of Operations Research, 3, 449-472.
 Molloy, M.K. (1982) Performance Analysis Using Stochastic Petri Nets. IEEE Transactions on Computers, 31, 913-917.
 Barad, M. (1994) Decomposing Timed Petri Nets of Open Queueing Networks. Journal of the Operational Research Society, 45, 1385-1397.
 Barad, M. (2003) An Introduction to Petri Nets. International Journal of General Systems, 32, 565-582.
 Barad, M. (1998) Timed Petri Nets as a Verification Tool. Proceedings of Winter Simulation Conference, Washington DC, 13-16 December 1998, Vol. 1, 547-554.
 Sargent, R.G. (1991) Simulation Model Verification and Validation. Proceedings of the Winter Simulation Conference, Phoenix, 8-11 December 1991, 37-47.
 Balci, O. (1994) Validation, Verification and Testing Techniques throughout the Life Cycle of a Simulation Study. Annals of Operations Research, 53, 121-173.
 Kelton, W.D. (1989) Random Initialization Methods in Simulation. IIE Transactions, 21, 355-367.
 Kleijnen, J.P.C. (1995) Verification and Validation of Simulation Models. European Journal of Operational Research, 82, 145-162.