Received 22 March 2016; accepted 23 May 2016; published 26 May 2016

ABSTRACT

Petri Nets (PNs) are an effective structure for modeling and analyzing asynchronous systems with concurrent and parallel activities. A Petri net models the static properties of a discrete event system concentrating on two basic concepts: events and conditions. Most of the theoretical work on Petri nets is a formal definition of Petri nets structures, which consist of a set of places, repre- senting conditions, a set of transitions, representing events, an input function and an output function. For practical purposes, a graphical representation is more useful. Two types of nodes portray places and transitions. A circle is a place and a bar is a transition. There is no inherent measure of time in a classical Petri net. To approach time-based evaluation of system performances, Timed Petri Nets (TPNs) were introduced. Modeling the notion of time is not straightforward. There are several possibilities for introducing time in PNs, among them timed transitions and timed places. This paper reviews several published examples where Petri Nets were used in different circumstances such as estimating expected utilization of processing resources at steady state in open queueing networks, verifying computerized simulations and batch planning in textile industry.

Keywords:

Petri Nets, Timed Petri Nets, Open Queueing Networks, Simulation Verification, Textile Industry

1. Introduction

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 [1] - [3] .

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. [4] [5] 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:

,

We see that the original marking of the net contains one token in place P₀. The firing of transition T₁ will insert one token in P₁. As P₀ and P₁ are the two input places of transition T₃, they will enable its firing. Firing of transition T₃ will remove the tokens in P₀ and P₁ and insert a token in P₃. Then, transition T₄ will fire and will thus remove the token in P₃ and insert a token in P₀. This process describes a firing sequence. Another firing sequence will take place if following the current initial state of the system, transition T₂ fires introducing a token in place P₂. The respective tokens in places P₀ and P₂ will enable the firing of transition T₅ that will remove the tokens in P₀ and P₂ and will add a token to place P₄. This will be followed by the firing of transition T₆ removing the token in P₄ and inserting a token in P₀.

Assume now that the initial state of the net contains not one token but three tokens, one in each of the three places P₀, P₁ and P₂. This situation represents a conflict as both transitions T₃ and T₅ are enabled. Which transition should fire, T₃ or T₅? 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.

Incidence Matrix

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 [4] and adopted by others e.g. [6] , Timed Petri Nets (TPNs) associate a positive number s(T) to each transition T. In [5] 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 [7] . A Timed Petri Net (TPN) graph looks exactly like a PN graph, except that tokens in timed places become available only after z_i time units (z_i > 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 P₃ or P₄ are as described above, meaning that a token in P₃ or in P₄ 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 T₁ or T₂, respectively represents input of type 1 or type 2 customers to the system. A token in place P₁ or P₂, respectively, represents a type 1 or type 2 customers waiting to be served by the resource. A token in place P₀ shows that the resource is available. Firing of transition T₃ or T₅, respectively, expresses the beginning of servicing a type 1 or type 2 customers. Firing of transition T₄ or T₆ expresses the end of service. This occurs after z₃ or z₄, time units have respectively elapsed. Then, one token returns to P₀, 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 = z_i, i = 0, 1, 2, 3, 4 for i = j and zero otherwise. The timed places are z₃ and z₄ respectively representing service time of type 1 and type 2.

The graphical representation of TPN1 is also Figure 1 to which z₃ and z₄ 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 [8] 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 [5] 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 [10] . 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 [9] and repeated in [5] , 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:

(1)

(2)

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,

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,

(3)

It is thus seen that under steady state conditions the entering (exogenous) flow of type 1 parts to the system, equals I₄, the firing rate of transition T₄ 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 I₆, 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 J₁,

(4)

The “1” values in J₁ represent all the possible steady states of workstation 1: P₀-idle, P₃-processing part type 1, P₄-processing part type 2.

Substituting Equation (4) in Equation (2) and using Equation (3) leads to:

(5)

Q₀, Q₃, Q₄ represent the token content of places P₀, P₃, P₄ 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:

(6)

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, z₀ is the station expected idle time between any two consecutive arrivals of parts. Any feasible solution has to satisfy z₀ > 0; leading to:

(7)

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 r_m 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.

(8)

where and represent the respective contributions of, the input flow of part types j, and that of I_f, the flow of disturbances, to the m station utilization, ,

(9)

is the expected operation duration of part j at station m.

is the expected duration of a disturbance.

Provided the constraint r_m < 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 r_m 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 [11] .

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. [12] - [14] .

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. [15] ). 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 [16] . 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 A_j (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 r_m, 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 [17] .

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 [8] , 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.

(10)

for,.

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 r_m 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.

5.2.2. Conclusion

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.

6. Discussion

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.

Cite this paper

Miryam Barad, (2016) Petri Nets—A Versatile Modeling Structure. Applied Mathematics,07,829-839. doi: 10.4236/am.2016.79074

References