A Traffic Model for Machine-Type Communications Using Spatial
Point Processes
Henning Thomsen, Carles Navarro Manchón, Bernard Henri Fleury

arXiv:1709.00867v1 [cs.NI] 4 Sep 2017

Department of Electronic Systems, Aalborg University, Denmark
Email: {ht, cnm, fleury}@es.aau.dk

Abstract—A source traffic model for machine-tomachine communications is presented in this paper. We
consider a model in which devices operate in a regular
mode until they are triggered into an alarm mode by
an alarm event. The positions of devices and events are
modeled by means of Poisson point processes, where the
generated traffic by a given device depends on its position
and event positions. We first consider the case where
devices and events are static and devices generate traffic
according to a Bernoulli process, where we derive the
total rate from the devices at the base station. We then
extend the model by defining a two-state Markov chain
for each device, which allows for devices to stay in alarm
mode for a geometrically distributed holding time. The
temporal characteristics of this model are analyzed via the
autocovariance function, where the effect of event density
and mean holding time are shown.

I. I NTRODUCTION
Machine-type communications (MTC), also called
machine-to-machine (M2M) communications, is a type
of transmission which does not directly involve human intervention. Examples include smart metering,
telemetry and sensors [1]. This type of transmission is
projected to be an important part of fifth generation (5G)
communication systems [2].
It is expected that the number of Machine-type Devices (MTDs) will be at least two orders of magnitude
higher than the number of human-type users in a given
cell [3]. Also, MTC traffic from several devices can be
correlated in time and space, since e.g. sensors can react
to external events such as a fire or traffic accident, which
have a spatial characteristic. This sort of behaviour can
result in traffic spikes or bursts which can put stress on
the receiver (e.g. a base station (BS)) since MTC traffic
is usually uplink-dominated. It is therefore desirable to
have tractable models that can be used to assess the
impact of this type of traffic on the BS.
In this paper, we formulate a traffic model based on
tools from stochastic geometry that takes into account
the spatial aspect of devices and events. Motivated by the
observation that traffic from individual MTDs is often
small, while the total traffic at the BS can be large [1],
we derive analytical expressions for the expected total
rate of the devices. It is assumed that each device can
be in two modes, regular and alarm, depending on the
proximity of the device location to the event locations.
978-1-5386-3531-5/17/$31.00 c 2017 IEEE

We also derive a Markov chain model and provide an
analytical expression which is a good approximation for
the expected total rate when the time spent in the alarm
state is low.
II. BACKGROUND AND R ELATED W ORKS
Models for M2M traffic can broadly be categorized
into two groups, source models and aggregated models [4]. In a source model, each source has its own set
of tunable parameters which control its individual traffic
characteristics. In an aggregated traffic model all sources
have a common set of parameters. A source model thus
allows for more fine-grained control; however this comes
at a higher cost in computational complexity because of
a higher number of parameters.
Traffic modeling for MTC has been considered by the
3rd Generation Partnership Project (3GPP) in [5]. They
propose two models, one accounting for uncoordinated
arrivals and one accounting for coordinated arrivals. In
the former, the packet arrivals are specified to follow a
uniform distribution in [0, Tu ], where Tu is 60 sec. The
latter model uses a scaled Beta distribution over [0, Tb ],
where Tb is 10 sec, with parameters α = 3 and β = 4.
In [6], other values of α and β are derived numerically
for different types of alarm events.
Modeling of human-type traffic (such as telephone
traffic) is usually done by means of Poisson processes
(PPs). However, PPs are not well suited for modeling coordinated traffic, which is expected to play a
large role in MTC applications. To address this, many
papers consider Markov modulated Poisson processes
(MMPPs). In such a model, a source can be in a number
of states S = {s1 , . . . , sm } and the state transitions are
governed by an irreducible Markov chain. When the process is in state si ∈ S , it generates Poisson-distributed
traffic with rate λi . A variant of this is the Markov
modulated Bernoulli process (MMBP), where Bernoullidistributed traffic with state-dependent probability pi is
generated [7]. Another model used especially in Internet
traffic modeling is the Markovian arrival process (MAP).
A MAP can be in a number of states si ∈ S , also
termed phases [7]. The process can then either enter
another phase with rate γ, or, with rate λ , enter another
phase and generate a packet. Contrary to the MMBP
and MMPP, where the interarrival times of packets are
geometrically and exponentially distributed respectively,

the MAP allows for more general interarrival time
distributions.
To take into account possible dependency between
MTDs, [8] considers a model where each device can be
in regular or alarm state. Two state-transition matrices
PU and PC are defined, the former representing uncoordinated behaviour and the other, coordinated behaviour.
The model is a source model where each device has its
own transition matrix which is a convex combination
of PU and PC . The coordination between devices is
modeled via a parameter δi for each device i, as well
as a global background process θ (t), which takes values
uniformly in [0, 1]. It is shown that the model has a lower
computational complexity compared to a full source
model. This model is extended in [4] and used for
evaluating a random access scheme supporting humantype and MTD traffic.
A spatio-temporal traffic model for MTDs using
queueing theory and stochastic geometry is analyzed
in [9]. They compare three different access protocols for
MTC and derive success probabilities, average queueing
delays and waiting times. Stability regions for the three
protocols are given, which can be used for scalability
assessment. Compared to this work, our model allows
for devices to be triggered by external events.
Reference [6] formulates a spatial correlation model
of MTDs using the IEEE802.11ah protocol. As in other
works, each device can be in two states, regular or alarm,
with transition probabilities described by the model
in [8]. In [6], they define the background process θ (t)
to be a function of alarm onset and velocity, as well as
location. They give a numerical expression for the total
rate, albeit not in closed form.
A two-state traffic model for automotive M2M applications based on a coupled MAP is proposed in [10].
Here, any car can trigger an event; then other cars go
into event mode if the distance to the car triggering the
event is less than a specified threshold. They show that
this model better captures burstiness than a PP model
In the present paper, we use the theory of spatial
Poisson point processes (PPPs) to model the position of
devices and event epicenters. We also define a function
which models the influence of events on device traffic. Even though spatial modeling has been considered
in [6], [10], our use of PPPs allows us to obtain tractable
analytical expressions for the total rate, and thus avoids
computationally expensive simulations. Also, we include
the temporal aspect in the model via defining a Markov
chain for each device, derive an approximation of the
expected total rate, and show via numerical simulations
that this approximation is quite close. We study the
temporal correlation of the traffic in the Markov chain
model via the autocovariance function. Even though our
model is a source model, the theory of PPPs can still be
applied and provides accurate results.

Figure 1: Conceptual deployment of devices in a cell,
with a BS at the center. Event epicenters are indicated
by red crosses. The influence of an event is coded in the
shade intensity: The darker the shade, the stronger the
influence.
III. P ROPOSED M ODEL
Let C be the coverage disk of a BS (located at the
origin of the disk) of radius R and let N devices be
deployed randomly and independently in C. Each MTD
can be in one of two states: regular (R) and alarm (A).
Time is slotted. In timeslot k, MTD i generates traffic
with rate Ri (k), depending on its current state. In the
regular state, it generates traffic with rate RR,i , and in
alarm state, the rate is RA,i . In this work, all devices are
of the same type, i.e. RR,i = RR and RA,i = RA , ∀i =
1, . . . , N.
A. Modeling the Device Locations and Event Epicenters
The MTDs are deployed in the coverage disk of the
BS according to a homogeneous PPP ΦM with density
λM . The event epicenters are represented by a homogeneous PPP ΦE with density λE in the Euclidean plane.
The processes ΦM and ΦE are assumed independent.
We choose to use PPPs because typical nodes can
be reasonably assumed to be randomly deployed in
the plane, in particular since we are not targeting a
specific application. Furthermore, this choice allows for
analytical tractability and can serve as a good indicator
of performance.
Fig. 1 shows an example of a deployment of devices
in a cell. Event epicenters are indicated with red crosses.
Note that events from outside the BS coverage area
can influence devices inside it. In the figure, three
devices are in alarm mode (transmission indicated with
solid red arrows) and two devices are in regular mode
(transmission indicated with dashed blue arrows).
B. Modeling Alarm Probabilities
In order to capture the effect of a given event on a
device, we define a function with gives the probability p
that an alarm is triggered in a device at location x by an
event with epicenter at location y in the Euclidean plane
R2 . This function depends on the distance between the
two locations. The definition is given below.

Definition (Alarm Triggering Probability Function
(ATPF)). An ATPF is a function f defined as
f : [0, ∞) → [0, 1],

d 7→ f (d) = p,

ΦM

ΦE

x

y

(1)

having a finite first moment.
In applications, f (d) is typically non-increasing to
represent a decaying influence of events on devices as
the distance d increases.

|ΦE |

C. Modeling Discrete Rates
Sx

Let Exy be the event that a device at location x is
triggered into alarm mode by an event with epicenter
at y. Let Exy be the complement of Exy and let pxy =
Pr{Exy } be the probability of Exy . The probability of
this device being in alarm mode is
px = Pr{Device at x is triggered by at least one event}
= 1 − Pr{No event triggers device at x}
)
(
\

Exy = 1 − ∏ Pr Exy
= 1 − Pr
x∈ΦE

= 1−

∏ (1 − pxy ),

(2)
(3)

y∈ΦE

(4)

Rx
|ΦM |
Figure 2: Dependence of the rate of the devices on the
two PPPs, illustrated as a directed graphical model.
with px given in (4). Then, each realization of ΦE and
ΦM results in state sequences {Sx (k)}∞
k=1 , which are
i.i.d. Bernoulli processes.

y∈ΦE

where in (3), we used the assumption that alarms
resulting from different events are independent. Note
that 0 ≤ 1 − f (dxy ) ≤ 1, where dxy = kx − yk and
pxy = f (dxy ).
The state of a device at location x, Sx is a random
variable which depends on ΦM and ΦE . The statistical
dependence between Sx , the corresponding rate Rx and
ΦM and ΦE is represented as a directed graphical model,
as shown in Fig. 2. In the figure, the blue box (called a
plate, see e.g. [11, Ch. 8]) around x, Sx and Rx indicates
|ΦM | instances of each of them, where |ΦM | is the
cardinality of ΦM . Similarly, the right plate around y
indicates an infinite but countable number of instances
of y. Note that Sx depends on an infinite number of
parameters, since ΦE contains an infinite number of
points almost surely. However, Sx is still well-defined, as
long as the ATPF is chosen according to the definition
and R < ∞.
Since each device at location x has its own probability
px of going into alarm mode, the present model can be
considered a source model. Also, the model allows for
correlated behaviour, since for two devices at locations x
and x0 respectively, where x 6= x0 , both Sx and Sx0 depend
on all points of ΦE .
We consider two ways of modeling the states Sx (k),
with slot index k. Below, we model them as i.i.d.
Bernoulli random variables, later in Sec. V we model
them as a Markov chain. For the former case,
(
Regular with prob. 1 − px
Sx (k) =
,
(5)
Alarm
with prob. px

IV. D ERIVATION OF THE E XPECTED T OTAL R ATE
A. Bernoulli Traffic Generation
In this subsection, we assume that each device generates traffic according to a Bernoulli process.
The expected rate of the device at location x is
Rx = px RA + (1 − px )RR .

(6)

By substituting the expression of px in (4), we get
Rx = px RA + (1 − px )RR

(7)
!

=

1−

∏ (1 − pxy )

RA +

y∈ΦE

= RA + (RR − RA )

∏ (1 − pxy )RR

(8)

y∈ΦE

∏ (1 − pxy ).

(9)

y∈ΦE

The total rate RT = ∑x∈ΦM Rx of all devices is then
!

∑

x∈ΦM

RA + (RR − RA )

∏ (1 − pxy )

.

(10)

y∈ΦE

The goal is now to obtain an expression for the total rate
averaged over ΦM and ΦE . This expression is given in
the following theorem:
Theorem. Let ΦM and ΦE be independent homogeneous
PPPs with densities λM and λE respectively, and f be
a valid ATPF. Let C be the circular coverage disk of

the BS, and R its radius. Then the expected total rate
averaged over ΦM and ΦE is
RT , E [RT ] = λM πR2 ·


Z
Z
RA + (RR − RA ) exp −2πλE

1 − px



∞

f (r)rdr

(11)
Proof. Taking the expectation over ΦM and ΦE in (10),
we get
"
!#
RA + (RR − RA )

∑

=E

#

= λM

∑

C

(12)
#

"

RA + E

∑

(RR − RA )

∏ (1 − pxy )

y∈ΦE

x∈ΦM

x∈ΦM

Z

∏ (1 − pxy )

y∈ΦE

x∈ΦM

"

RA + (RR − RA )E

∑ ∏

(1 − pxy )

(13)

q

= λM πR2 RA
"

"

∑

##

∏ (1 − pxy )

E

px
Figure 3: State transition diagram of the Markov chain
model describing the temporal behaviour of the ith
device.
It is easy to see that (20) satisfies the conditions in the
Definition of an ATPF. In this case, pxy = exp(−dxy ).
The expected total rate is
(21)

which follows since 0∞ exp(−r)rdr = 1. Using (21),
we can see that λE → 0 (i.e. no event) implies RT =
λM πR2 RR (only regular traffic), while λE → ∞ implies
RT = λM πR2 RA (only alarm traffic), as expected.
R

x∈ΦM y∈ΦE

, (14)

y∈ΦE

x∈ΦM

where in (13) Campbell’s Theorem [12, Thm. 4.1] was
used for a PPP with density λM over a circular disk C
of radius R.
To compute the inner expectation in the second term
of (14), we use the probability generating functional
(PGFL) of a PPP with density λE [12, p.125]. After
changing to polar coordinates, we get for each summand
in (14)
"
#


Z ∞
E ∏ (1 − f (dxy )) = exp −2πλE
f (r)rdr .
0

y∈ΦE

(15)

Then, the outer expectation in (14) equals
"

#
Z ∞
E ∑ exp −2πλE
f (r)rdr .
x∈ΦM

Alarm

RT = λM πR2 (RA + (RR − RA ) exp(−2πλE )),

#

"

+ (RR − RA )E

Reg.

.

0

C

E

1−q

(16)

0

V. M ARKOV CHAIN A NALYSIS
In this section, we extend the model, enabling each
device to stay in alarm mode for a duration geometrically distributed with parameter q.
We characterize the state of a device at location x ∈
ΦM at a given time k as a Markov process, i.e.
Pr{Sx (k) | Sx (k − 1), . . . , Sx (0)} = Pr{Sx (k) | Sx (k − 1)}
(22)
where Sx (k) can take two different values, Regular and
Alarm. The state transition diagram is shown in Fig. 3,
and the state transition matrix is


1 − px px
Px =
.
(23)
1−q q

Using Campbell’s Theorem, the expectation in (16) can
be expressed as


Z
Z ∞
f (r)rdr
(17)
λM exp −2πλE
0
C


Z ∞
= λM πR2 exp −2πλE
f (r)rdr .
(18)

This Markov chain is ergodic; it has a unique steadystate probability vector πx = [πx,A , πx,R ], where πx,A
(πx,R ) is the probability of alarm (regular) state. The
vector πx can be found by solving the system of linear
equations

Combining the above results and rearranging yields



Z ∞
2
λM πR RA + (RR − RA ) exp −2πλE
f (r)rdr
.

together with the normalizing condition πx,A + πx,R = 1:

πx Px = πx ,

0

0

(19)

Example (Exponential ATPF). Consider an exponentially decaying function of the distance d:
f (d) = exp(−d).

(20)

πx,A =

px
,
1 + px − q

πx,R =

(24)

1−q
.
1 + px − q

(25)

We assume that the initial state Sx (0) is distributed
according to the above steady-state distribution.
Note that by setting q = px , we recover the model
defined in (5). Also, by replacing each px in (25)
with p , E [px ], we can approximate the total rate

We comment on this approximation in Sec. VI. Note that
by setting q = 0 and multiplying px by a factor θ (t) in
the Markov chain, we recover the models of [4], [8].

102

101

10-4

10-3

10-2

Event density

We study the temporal characteristics of the Markov
chain model through the autocovariance function (ACF)
of the total rate.
The reason for introducting the ACF is that we want
to see the effect of the parameters λE , λM and q
on the long-term temporal dependence of the model.
We suspect that by letting q → 1, the autocovariance
function will decay slowly, since devices triggered by
an alarm will stay in the alarm mode for a longer period
of time, implying temporally correlated traffic. The ACF
of R(k) is defined as

1 
E (R(k) − µR )(R(k0 ) − µR ) ,
2
σR

Analysis
Simulation

100
10-5

A. Temporal Characteristics of the Model

CRR (k, k0 ) =

103

Expected total rate [bps]

R(k) = ∑x∈ΦM Rx (k) of the Markov chain model in a
low complexity manner:
"

#
1
−
q
p
fT = E ∑
RA +
RR
(26)
R
1+ p−q
x∈ΦM 1 + p − q


p
1−q
2
= λM πR
RA +
RR .
(27)
1+ p−q
1+ p−q

(28)

where k, k0 = 1, 2, . . . and µR and σR2 are the mean and
variance of R(k) respectively. The expectation is taken
over the distribution of R(k).
Note that we use the autocovariance function instead
of the autocorrelation, because we want to compare
the temporal dependence for various settings of the
parameters λE , λM and q on the same scale.
We estimate the ACF of R(k) via Monte Carlo (MC)
simulations as follows. For any trial, the unbiased sample autocovariance
NS −k−1
1
∑ (R(i) − µR )(R(i + k) − µR ),
2
(NS − k)σR i=1
(29)
is computed. Here NS is the number of time samples.
Note that in our case, (29) is conditioned on realizations
of ΦE and ΦM . Also, the mean and variance are replaced
by their respective sample estimates. Averaging CbRR (k)
over the PPP realizations via MC simulations, we get
CRR (k). In the next section, the role of this quantity on
the traffic will be discussed. The derivation of the ACF
is left for future work.

CbRR (k) =

VI. N UMERICAL R ESULTS
In this section, we show some quantitative results
illustrating the behavioural features of the Bernoulli
and Markov chain models, using MC simulations done
in Matlab. In the simulations, devices and events are
deployed independently and traffic is generated either

10-1

100

101

E

Figure 4: Expected total rate for Bernoulli traffic.
using the Bernoulli or Markov chain models. The resulting total rate is then averaged over event and device
realizations. We use the exponential ATPF (20) in all
simulations. The simulation parameters are listed in
Tab. I. They are used unless otherwise specified.
A. Bernoulli Process Model
Fig. 4 shows a comparison between simulation and
analytical results for the Bernoulli model. We see that
the theoretical and simulation curves match, which validates our model. Also, we see that for λE → 0 only
regular traffic occurs and λE → ∞ implies that only alarm
traffic is received at the BS. Three regions are identified,
two saturation regions (when λE → 0 and λE → ∞) and
a transition region.
B. Markov Chain Model
The Markov chain model in Sec. V is numerically
evaluated in this subsection. We compare the total rate
computed from the simulations against the approximafT in (27) and the analytical expression for RT
tion R
in (11). We chose to include RT since we want to
compare the two expressions and identify the numerical
values of q for which they are close.
In Fig. 5 the expected total rate of the Markov chain
model is shown for two values of λE . It can be seen
that for the Bernoulli model, the rate is constant, since it
does not depend on q. As for the Markov chain model,
the rate obtained from the simulation is always upper
bounded by the steady-state approximation, where the
latter is obtained from (27). Further, these two curves
have the same shape. The approximation can be used as
an upper bound on the total rate at the BS; however the
bound is not tight.
Table I: Simulation parameters.
Parameter
s
λM
λE
RA
RR
R

Description
Observation window size
Device density
Event density
Rate in alarm mode
Rate in regular mode
Cell radius

Value
100 m
10−1 nodes/m2
10, . . . , 10−5 events/m2
1 packet/slot
0.01 packet/slot
20 m

55
50

Simulation, Markov chain
Analysis, Bernoulli
Steady State approx., Markov chain

Expected Total Rate [bps]

45
40
35

E

= 10-1

30
25
20
15
E

10
5
0.1

0.2

0.3

0.4

= 10-3

0.5

0.6

0.7

0.8

0.9

Probability q

Figure 5: Expected total rates for Bernoulli traffic and
for the Markov chain model; λM = 10−1 .
1
q=0.1
q=0.5
q=0.9

0.9

VIII. ACKNOWLEDGEMENTS

0.8
0.7

Sample ACF

model enables tuning of the model to suit various MTC
applications and alarm reporting strategies. The models
can be used to assess network performance in terms of
the total rate and to see the impact of the event density.
Furthermore, the Markov chain model, along with the
tunable parameter q, can be used to study the impact of
events on the temporal correlation of the total traffic at
the BS. Also, even though our model is a source model,
the expected total rate is still fast to compute, thanks
to the machinery of PPPs. As future work, we plan to
include time-dependency in the event modeling, as well
as to consider several types of MTDs and events. Also,
we plan to model other types of sources, such as Poisson
and periodic sources, and study the variance of the total
traffic, along with its impact on network performance.
This work has been supported by the cooperative
project VIRTUOSO, partially funded by Innovation
Fund Denmark.

0.6
0.5
0.4

R EFERENCES

0.3
0.2
0.1
0
5

10

15

20

25

Lag k

Figure 6: Sample ACF of the total rate for the Markov
chain model; λE = 10−2 , λM = 10−2 .
The sample autocovariance CRR (k) is shown in Fig. 6.
Focusing first on the case q = 0.1, we see that the ACF
almost behaves like the one from an i.i.d. Bernoulli
process, in that it has short memory. This is because
the px s are low (through the relatively low densities of
the PPPs) and because of the relatively small value of
q. Increasing the value of q, see the dashed red (for
q = 0.5) and solid black (for q = 0.9) curves, increases
the memory, i.e. the total rate at a given time k is
correlated with many past values. This is because the
Markov chain now stays in the alarm state longer when
it enters that state.
Hence, the introduction of the additional parameter
q in the Markov chain model allows for tuning the
temporal correlation of the individual rate processes of
the MTDs and, as a result, that of the total rate process.
VII. C ONCLUSION AND F UTURE W ORK
In this paper, a MTC traffic model based on spatial point processes was given. We also provided a
verification of the model via simulations and derived
an expression for the total rate of the devices. This
expression is very fast to compute, compared to a
full simulation, especially for high device and event
densities. A Markov chain model was also defined. It
was seen that approximating the probability px by p for
all devices, gives a good approximation of the total rate
for low values of q. The parameter q of the Markov chain

[1] M. Z. Shafiq, L. Ji, A. X. Liu, J. Pang, and J. Wang, “Largescale measurement and characterization of cellular machineto-machine traffic,” IEEE/ACM Transactions on Networking,
vol. 21, no. 6, pp. 1960–1973, 2013.
[2] F. Boccardi, R. W. Heath, A. Lozano, T. L. Marzetta, and
P. Popovski, “Five disruptive technology directions for 5G,”
IEEE Communications Magazine, vol. 52, no. 2, pp. 74–80,
2014.
[3] B. Yang, G. Zhu, W. Wu, and Y. Gao, “M2M access performance
in LTE-A system,” Transactions on Emerging Telecommunications Technologies, vol. 25, no. 1, pp. 3–10, 2014.
[4] M. Centenaro and L. Vangelista, “A study on M2M traffic and
its impact on cellular networks,” in Internet of Things (WF-IoT),
2015 IEEE 2nd World Forum on, pp. 154–159, IEEE, 2015.
[5] 3GPP TR 37.868 V11.0.0, Study on RAN Improvements for
Machine-Type Communications, 2011.
[6] G. C. Madueño, Č. Stefanović, and P. Popovski, “Reliable and
efficient access for alarm-initiated and regular M2M traffic in
IEEE 802.11ah systems,” IEEE Internet of Things Journal,
vol. 3, no. 5, pp. 673–682, 2016.
[7] O. Ibe, Markov processes for stochastic modeling. Newnes, 2013.
[8] M. Laner, P. Svoboda, N. Nikaein, and M. Rupp, “Traffic models
for machine type communications,” in Wireless Communication
Systems (ISWCS 2013), Proceedings of the Tenth International
Symposium on, pp. 1–5, VDE, 2013.
[9] M. Gharbieh, H. ElSawy, A. Bader, and M.-S. Alouini, “Spatiotemporal Stochastic Modeling of IoT Enabled Cellular Networks: Scalability and Stability Analysis,” IEEE Transactions on
Communications, vol. 65, no. 8, pp. 3585–3600, 2017.
[10] E. Grigoreva, M. Laurer, M. Vilgelm, T. Gehrsitz, and
W. Kellerer, “Coupled Markovian Arrival Process for Automotive Machine Type Communication Traffic Modeling,” in
Communications (ICC), 2017 IEEE International Conference on,
IEEE, 2017.
[11] C. M. Bishop, Pattern Recognition and Machine Learning.
Springer, 2006.
[12] S. N. Chiu, D. Stoyan, W. S. Kendall, and J. Mecke, Stochastic
geometry and its applications. John Wiley & Sons, 2013.

