An End-to-End Performance Analysis for Service Chaining in a Virtualized Network

TL;DR

This paper presents an analytical framework for end-to-end performance analysis of service chaining in virtualized networks, incorporating MEC and core network servers, to aid in traffic management and system optimization.

Contribution

It introduces a queueing and stochastic modeling approach for analyzing performance in service chaining scenarios with virtual network functions.

Findings

Analytical and simulation results closely match.

Performance metrics vary with different traffic and system configurations.

Insights support traffic flow control decisions.

Abstract

Future mobile networks supporting Internet of Things are expected to provide both high throughput and low latency to user-specific services. One way to overcome this challenge is to adopt Network Function Virtualization (NFV) and Multi-access Edge Computing (MEC). Besides latency constraints, these services may have strict function chaining requirements. The distribution of network functions over different hosts and more flexible routing caused by service function chaining raise new challenges for end-to-end performance analysis. In this paper, as a first step, we analyze an end-to-end communications system that consists of both MEC servers and a server at the core network hosting different types of virtual network functions. We develop a queueing model for the performance analysis of the system consisting of both processing and transmission flows. We propose a method in order to derive analytical expressions of the performance metrics of interest. Then, we show how to apply the similar method to an extended larger system and derive a stochastic model for such systems. We observe that the simulation and analytical results coincide. By evaluating the system under different scenarios, we provide insights for the decision making on traffic flow control and its impact on critical performance metrics.

Equations

$$P_{i,j:u,k}=\mathrm{Pr}\left\{X_{n+1}=i\text{, }{Y}_{n+1}=j|X_n=u,Y_n=k\right\}\text{,}$$

$$\displaystyle\mathbf{P}_{1}=\left[\begin{array}[]{ccccc}\mathbf{B}&\mathbf{C}&&&\\ \mathbf{E}&\mathbf{A}_{1}&\mathbf{A}_{0}&&\\ &\mathbf{A}_{2}&\mathbf{A}_{1}&\mathbf{A}_{0}&\\ &&\ddots&\ddots&\\ &&&\mathbf{A}_{2}&\mathbf{A}_{0}+\mathbf{A}_{1}\end{array}\right]\text{.}$$

$$\displaystyle\mathbf{B}=\left[\begin{array}[]{ccc}\overline{\lambda}_{1}&0&0\\ b_{2}&b_{1}&0\\ 0&b_{2}&b_{1}\end{array}\right]\text{, }\mathbf{C}=\left[\begin{array}[]{ccc}\lambda_{1}&0&0\\ c_{2}&c_{1}&0\\ 0&c_{2}&c_{1}\end{array}\right]\text{, }\mathbf{E}=\left[\begin{array}[]{ccc}0&\overline{\lambda}_{1}\mu_{1}&0\\ 0&\overline{\lambda}_{1}\mu_{1}\mu_{2}&\overline{\lambda}_{1}\mu_{1}\bar{\mu}_{2}\\ 0&0&\overline{\lambda}_{1}\mu_{1}\end{array}\right]\text{,}$$

$$\displaystyle\mathbf{A}_{0}=\left[\begin{array}[]{ccc}\lambda_{1}\bar{\mu}_{1}&0&0\\ a_{2}^{(0)}&a_{1}^{(0)}&0\\ 0&a_{2}^{(0)}&a_{1}^{(0)}\end{array}\right]\text{, }\mathbf{A}_{1}=\left[\begin{array}[]{ccc}\overline{\lambda}_{1}&\lambda_{1}\mu_{1}&0\\ a_{2}^{(1)}&a_{1}^{(1)}&a_{0}^{(1)}\\ 0&a_{2}^{(1)}&a_{1}^{(1)}\end{array}\right]\text{, }\mathbf{A}_{2}=\left[\begin{array}[]{ccc}0&\bar{\lambda}_{1}\mu_{1}&0\\ 0&a_{2}^{(2)}&a_{1}^{(2)}\\ 0&0&a_{2}^{(2)}\end{array}\right]\text{.}$$

$$\displaystyle\mathbf{P}_{1}^{(1)}=\left[\begin{array}[]{cccc}1&0&&\\ \mu_{2}&\bar{\mu}_{2}&&\\ &\ddots&\ddots&\\ &&\mu_{2}&\bar{\mu}_{2}\end{array}\right]\text{, }\mathbf{P}_{1}^{(2)}=\left[\begin{array}[]{ccccccc}0&1&0&&\\ 0&\mu_{2}&\bar{\mu}_{2}&&\\ 0&0&\mu_{2}&\bar{\mu}_{2}&\\ &&&\ddots&\ddots&\\ &&&&0&1\end{array}\right]\text{.}$$

$$\mathbf{B}={\bar{\lambda }}_1{\mathbf{P}}_1^{(1)}\text{, }\mathbf{C}={\lambda }_1{\mathbf{P}}_1^{(1)}\text{, }\mathbf{E}={\bar{\lambda }}_1{\mu }_1{\mathbf{P}}_1^{(2)}\text{,}$$

$${\mathbf{A}}_0={\lambda }_1{\bar{\mu }}_1{\mathbf{P}}_1^{(1)}\text{, }{\mathbf{A}}_1={\bar{\lambda }}_1{\bar{\mu }}_1{\mathbf{P}}_1^{(1)}+{\lambda }_1{\mu }_1{\mathbf{P}}_1^{(2)}\text{, }{\mathbf{A}}_2={\bar{\lambda }}_1{\mu }_1{\mathbf{P}}_1^{(2)}\text{.}$$

$${\lambda }_3=p(1-\alpha )\text{ and }{\lambda }_4=\text{Pr}\left\{Q_3>0\right\}{\mu }_3=\left(\sum _{j=0}^{M_4}\sum _{i=1}^{M_3}{\pi }_{i,j}^{(2)}\right){\mu }_3\text{, respectively.}$$

$$\displaystyle\mathbf{P}_{3}=\left[\begin{array}[]{cccccc}\bar{\lambda}_{5}&\lambda_{5}&&\\ \bar{\lambda}_{5}\mu_{5}&\lambda_{5}\mu_{5}+\bar{\lambda}_{5}\bar{\mu}_{5}&\lambda_{5}\bar{\mu}_{5}&\\ &\bar{\lambda}_{5}\mu_{5}&\lambda_{5}\mu_{5}+\bar{\lambda}_{5}\bar{\mu}_{5}&\lambda_{5}\bar{\mu}_{5}\\ &\ddots&\ddots&\ddots&\\ &&\bar{\lambda}_{5}\mu_{5}&\bar{\lambda}_{5}\bar{\mu}_{5}+\lambda_{5}\end{array}\right]\text{.}$$

$${\pi }_i^{(3)}=\frac{{\lambda }_5^i{\bar{\mu }}_5^{(i-1)}}{{\bar{\lambda }}_5^i{\mu }_5^i}{\pi }_0^{(3)}\text{, }\text{for }i=1\text{, }\dots \text{, }{M}_5\text{ and }{\pi }_0^{(3)}={\left[1+\sum _{i=1}^{M_5}\frac{{\lambda }_5^i{\bar{\mu }}_5^{i-1}}{{\bar{\lambda }}_5^i{\mu }_5^i}\right]}^{-1}\text{.}$$

$$p_{00}$$

$$b_1$$

$$b_2$$

$$b_3$$

$$b$$

$$\displaystyle\mathbf{P}_{4}=\left[\begin{array}[]{ccccccc}p_{00}&p_{01}&p_{01}&\\ b_{0}&b_{1}&b_{2}&b_{3}&\\ &b_{0}&b_{1}&b_{2}&b_{3}\\ &&\ddots&\ddots&\ddots&\ddots\\ &&&b_{0}&b_{1}&b_{2}&b_{3}\\ &&&&b_{0}&b_{1}&b_{2}+b_{3}\\ &&&&&b_{0}&b_{1}+b_{2}+b_{3}\end{array}\right]\text{.}$$

$${\mathbf{\pi }}^{(4)}{\mathbf{P}}_4={\mathbf{\pi }}^{(4)}\text{, }{\mathbf{\pi }}^{(4)}1=1\text{.}$$

$$\displaystyle\mathbf{P}_{4}=\left[\begin{array}[]{ccccccccc}a_{0}&b_{0}&0&0&\cdots\\ a_{1}&b_{1}&b_{0}&0&\cdots\\ a_{2}&b_{2}&b_{1}&b_{0}&\cdots\\ 0&b_{3}&b_{2}&b_{1}&\cdots\\ &0&b_{3}&b_{2}&b_{1}&\\ &&\ddots&\ddots&\ddots&\ddots\end{array}\right]\text{,}$$

$$P_{D_1}$$

$$P_{D_3}$$

$$P_{D_6}=p_{02}(1-{\mu }_6){\pi }_{M_6-1}^{(4)}+(p_{01}+2p_{02}){\pi }_{M_6}^{(4)}(1-{\mu }_6)\text{.}$$

$${\bar{Q}}_1$$

$${\bar{Q}}_4$$

$$T_i={\lambda }_i-P_{D_i}\text{, }\forall i\in \mathcal{K}\text{,}$$

$$D_i=\frac{{\bar{Q}}_i}{T_i}+\frac{1}{{\mu }_i}\text{, }\forall i\in \mathcal{K}\text{.}$$

$$D_{\text{sys}}=\alpha (D_1+D_2)+(1-\alpha )(D_3+D_4+D_5)+D_6\text{.}$$

$$\displaystyle\mathbf{P}_{6}=\left[\begin{array}[]{ccccccccccc}p_{00}&p_{01}&p_{02}&p_{03}&p_{04}&&&&&&\\ p_{10}&p_{11}&p_{12}&p_{13}&p_{14}&p_{15}&&&&&\\ b_{6}&b_{5}&b_{4}&b_{3}&b_{2}&b_{1}&b_{0}&&&&\\ &\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&&&\\ &&b_{6}&b_{5}&b_{4}&b_{3}&b_{2}&b_{1}&b_{0}&&\\ &&&b_{6}&b_{5}&b_{4}&b_{3}&b_{2}&b_{1}&b_{0}\\ &&&&b_{6}&b_{5}&b_{4}&b_{3}&b_{2}&b_{0}+b_{1}\\ &&&&&b_{6}&b_{5}&b_{4}&b_{3}&b_{2}+b_{1}+b_{0}\\ &&&&&&b_{6}&b_{5}&b_{4}&b_{3}+b_{2}+b_{1}+b_{0}\\ &&&&&&&b_{6}&b_{5}&b_{4}+b_{3}+b_{2}+b_{1}+b_{0}\\ \end{array}\right]\text{,}$$

$$P_{D_6|Q_6=M_6}$$

$$+2p_{03}{\pi }_{M_6}^{(4)}2\mu (1-\mu )+p_{03}{\pi }_{M_6}^{(4)}{\mu }^2+4p_{04}{\pi }_{M_6}^{(4)}(1-\mu {)}^2+3p_{04}{\pi }_{M_6}^{(4)}2(1-\mu )\mu$$

$$+2p_{04}{\pi }_{M_6}^{(4)}{\mu }^2\text{,}$$

$$P_{D_6|Q_6=M_6-1}$$

$$+3p_{04}{\pi }_{M_6-1}^{(4)}(1-\mu {)}^2+2p_{04}{\pi }_{M_6-1}^{(4)}2(1-\mu )\mu +p_{04}{\pi }_{M_6-1}^{(4)}{\mu }^2\text{,}$$

$$P_{D_6|Q_6=M_6-2}$$

$$P_{D_6|Q_6=M_6-3}$$

$$P_{D_6}=P_{D_6|Q_6=M_6-3}+P_{D_6|Q_6=M_6-2}+P_{D_6|Q_6=M_6-1}+P_{D_6|Q_6=M_6}\text{.}$$

$$D_{B{S}_1}$$

Full text

Future mobile networks supporting Internet of Things are expected to provide both high throughput and low latency to user-specific services. One way to overcome this challenge is to adopt Network Function Virtualization (NFV) and Multi-access Edge Computing (MEC). Besides latency constraints, these services may have strict function chaining requirements. In other words, each service has to be processed by a set of network functions in a specific order. The distribution of network functions over different hosts and more flexible routing caused by service function chaining raise new challenges for end-to-end performance analysis. In this paper, as a first step, we analyze an end-to-end communications system that consists of both MEC servers and a server at the core network hosting different types of virtual network functions. We develop a queueing model for the performance analysis of the system consisting of both processing and transmission flows. To approximate the behavior of the system, we decompose the system into subsystems and analyze them independently. By doing so we are able to provide approximate analytical expressions of the performance metrics such as system drop rate, end-to-end delay, and system throughput. Then, we show how to apply the similar method to an extended larger system and derive a stochastic model for systems with arbitrary number of servers at the edge. Simulation results show that our approximation model is accurate for the considered systems. We see in Section VI that the simulation and analytical results coincide. By evaluating the system under different scenarios, we provide insights for the decision making on traffic flow control and its impact on critical performance metrics.

I Introduction

The increasing demand of different kinds of network services and the requirements of 5G networks for low capital expenditure raise the need for the improvement of today’s networks in terms of both flexibility and scalability. In conventional communications networks, network functions (e.g., firewalls, transcoders, load balancers, etc.) are performed as dedicated hardware middleboxes. Although dedicated hardware can provide high performance, it causes high capital expenditure, low flexibility and scalability, and dependence on particular application. In future communications networks, these limitations will be addressed by Network Function Virtualization (NFV) [2, 3, 4]. The idea of NFV is to decouple the network functions from dedicated hardware equipment. More specifically, general purpose servers can host one or more types of network functions. The network can be flexible and it can deploy proper network functions according to the demands of various traffic types while reducing the capital expenditure. Furthermore, mission-critical mobile applications, e.g., augmented reality, connected vehicles, eHealth, will provide services that require ultra-low latency [5], [6]. To satisfy such requirements, Multi-access Edge Computing (MEC) has been proposed as a key solution [5]. The idea of MEC is to locate more computation resources closer to the users, e.g., at the base stations. NFV together with MEC are considered key technologies for 5G wireless sytems. However, computation capabilities and available resources of MEC servers are still limited compared to the high-end servers in the cloud. Therefore, it is interesting to further investigate the cooperation between the edge and the core, and the cooperation among MEC servers.

I-A Related work

Recently, the study of the performance of networks in Virtual Network Function (VNF)/ Software Defined Network (SDN) environment has attracted a lot of attention, [7, 8, 9, 10, 11, 12, 13, 14]. Ye et al.[7] analyze the end-to-end delay for embedded VNF chains. They consider two types of services that traverse different VNF chains and provide the delay analysis for each different chain. Miao et al.[9] provide an analytical model based on Stochastic Network Caclulus (SNC) to provide upper and lower delay bounds of a VNF chain. In their analysis, they consider both the case of bursty and non-bursty traffic. Along similar lines, Duan [8] analyzes an end-to-end delay performance of service function chaining for particular services and given resources. Authors in [10, 11, 12, 13, 14], apply tools from queueing theory to evaluate the performance of systems in SDN environment. In particular, Jarschel et al.[12] study the OpenFlow architecture, where the switch is modeled as an M/M/1 queue and the controller as a feedback system of the delay-loss type M/M/1/S queue. Similarly Goto et al.[13] analyze a simple OpenFlow-based switch in SDN environment, however, they distinguish traffic from the controller and exogenous traffic. Furthermore, a reasonable amount of works consider the modeling of connected VNF as a sequence of queues where the goal is to guarantee the stability of the system and some particular network or service requirements. A well known and widely used mathematical tool for stabilizing dynamic systems is the Lyapunov optimization theory. The authors in [15, 16, 17, 18, 19, 20] develop dynamic algorithms by applying Lyapunov optimization theory in order to guarantee system stability and fulfill additional service requirements. In particular, in [15, 18] flows traverse VNF chains and each node decides the resource allocation to each VNF and routes the flow to the next node. Gu et al.[20] develop a dynamic distributed algorithm that controls the flow and rate at each node. Their objective is to achieve fairness between the services and maximize the network utility while providing system stability. Chen et al.[19] consider the problem where VNF are installed in Virtual Machines. Each VM can be located in the same or different data center. In this work, each VM decides which functions to install or uninstall and which services to serve. A dynamic algorithm is developed that works in a distributed manner and takes online decisions. Considering a more static and known environment, researchers investigate the VNF placement and resource allocation problems [21, 22, 23, 24, 25]. In these works, the authors formulate the VNF placement problem as Mixed Integer Linear Problem (MILP) under the assumption of known traffic demand. In addition, approximation algorithms have been developed in order to provide a solution to the VNF placement problem, see for example [26, 27].

Besides the VNF/SDN technologies, MEC technology promises a significant improvement of the networks, especially in terms of latency reduction. Task offloading from the mobile devices to MECs has attracted considerable attention. For example, the authors in [28, 29, 30, 31] address the trade-off between the power consumption and the task processing delay. The authors in [32], propose an optimization model to perform joint slicing of mobile network and edge computation resources. As VNF and MEC are two strongly connected technologies, the research evolution of the one affects the evolution of the other. To the best of our knowledge, there are no works providing end-to-end network performance analysis for a network with MEC deployed and operated within the NFV environment.

I-B Contributions

In this paper, we investigate a VNF-facilitated end-to-end network where the cooperation between the edge and the core network, and the cooperation among MECs are considered. As a first step, we model and analyze a simple end-to-end communications system which consists of two MEC servers at the edge network and one at the core network hosting different types of VNFs. We develop a methodology by applying tools from queueing theory in order to derive approximate analytical expressions of the performance metrics of interest such as end-to-end delay, drop rate, and throughput. The methodology is based on the decomposion of the system into smaller subsystems easier to analyze. The information from analysis of each subsystem can be used in order to study the performance of the entire system. Since, in reality, there are no buffers with infinite capacity, we consider finite size buffers in the servers of our system model. However, the capabilities of the server at the core are considered to be much higher than those of the MEC servers. In order to observe the behavior of the system asymptotically, we also investigate, as a subcase, the scenario where the buffers at core server is infinite. Based on the methodology derived for the simple model, we extend the analysis to a more complex system with a larger number of MEC servers. As shown in the simulation results, the proposed stochastic model is accurate and indicates a robust behavior even for larger number of servers. Furthermore, we provide results that show the trade-off between the throughput, drop rate, and end-to-end delay. In this work we provide approximate analytical expressions for the performance metrics of our interest in order to provide insights of how to design a system to satisfy particular network requirements.

II System Model

The aim of this work is to study the performance of an end-to-end communications system that consists of different types of VNFs. We first provide the analysis for a simple end-to-end communications system and then, we provide the analysis for more general topologies.

We consider an end-to-end communications system consisting of a mobile device, two MEC servers, and one server located in the core network as depicted in Fig. 1. A task traverses a service chain of two consecutive VNFs: VNF $1$ and VNF $2$. In the chosen system, a MEC server, called Server $1$, is co-located with the base station and hosts one copy of VNF $1$ as the primary MEC server. A secondary MEC server, called Server $2$, is located nearby and also hosts a copy of VNF $1$. Server $2$, for example, may be located at a peer edge host with spare capacity or at a central office location within a metropolitan area network. In addition, Server $3$ in the core network hosts VNF $2$ and has more advanced computation capabilities than Servers $1$ and $2$. The analysis of such system can be used to study the deployment of AR-related applications. For example, VNF function Domain called Name Service (DMS) can be installed in MEC to map user requests to their corresponding nodes providing contents, while the compute-intensive VNFs such as AR or 3D gaming can be deployed in core networks.

We assume a slotted time system. At each time slot, the device transmits a task in form of a packet to a base station over a wireless channel. Because of the presence of fading in the wireless channel, transmissions may face errors. Thus, we assume that a task is successfully transmitted to the base station with a probability $p$ that captures fading, attenuation, noise, etc. The device attempts for a new task transmission only if the previous task is successfully received at the base station. The received tasks need to be distributed between the queue for local processing and the queue for transmission to the secondary MEC server.

Thus, there are two possible routes to pass through the service chain. A flow controller at the base station decides randomly the routing for each task. With probability $\alpha$ the task is processed by Server $1$ first to be processed by VNF $1$, and then forwarded to Server $3$ to be processed by VNF $2$. With probability $1-\alpha$ the task is forwarded to Server $2$, to be processed by VNF $1$, and then forwarded to Server $3$ for being processed by VNF $2$.

Each task that arrives at a server first waits in a queue for being processed by a VNF. Then, after the processing, it is stored in the transmission queue, waiting to be forwarded and processed by the next VNF. Let $Q_{i}$ denote the $i$-th queue, where $i\in\mathcal{K}$, and $\mathcal{K}$ is the set of the queues in the system. Note that the queues follow an early departure-late arrival model: at the beginning of the slot the departure takes place and a new arrival can enter the queue at the end of the slot. The queues for task transmission are $Q_{2}$, $Q_{3}$, and $Q_{5}$, and the queues for task processing are $Q_{1}$, $Q_{4}$, $Q_{6}$. The arrival rates for queues $Q_{1}$ and $Q_{3}$ are $p\alpha$ and $p(1-\alpha)$, respectively. We denote by $\mu_{i}$, $i\in\mathcal{K}$, the service rates of the queues. We assume that the service times are geometrically distributed. Furthermore, given that $Q_{1}$, $Q_{3}$, and $Q_{4}$ are non empty, the arrival rates of $Q_{2}$, $Q_{4}$, and $Q_{5}$ are equivalent to the service rates of $Q_{1}$, $Q_{3}$, and $Q_{4}$ (i.e., $\mu_{1}$, $\mu_{3}$, and $\mu_{4}$) respectively.

Furthermore, the queues are assumed to have finite buffer. Let $M_{i}$ denote the buffer size of each queue $i\in\mathcal{K}$. If a queue is full and no task departs at the same time that a new one arrives, the new task is dropped and removed from the system. Besides the case presented in Fig. 1, we will consider the case of a scaled-up system in Section V. Note that the described system model can be considered as an isolated network slice under the network slicing and MEC technologies.

III Performance Analysis

In this section, we perform the modeling and the performance analysis that allow us to derive the critical performance metrics. We model the considered queueing system utilizing Discrete Time Markov Chain (DTMC). Modeling the whole system as one Markov chain can drive in a quite complicated system difficult to be analyzed in terms of closed-form expressions. Thus, in order to simplify the analysis, we decompose the system into different subsystems. Since the computation capabilities of the server in the core are much higher than those of the MEC servers, we consider, in the analysis, the core server as one independent system. In order to develop a model that couples MEC Server $1$ and MEC Server $2$, we consider as one independent system the subsystem which consists of $Q_{3}$ and $Q_{4}$. $Q_{5}$ is considered as one independent system. We choose to decompose the system as described above in order to study the decomposition strategies where derive: i) a model for the interacting queues within one server, ii) a model for two queues of different servers that are interacting, and iii) an individual model for each queue within the same server. To summarize, we consider the following four subsystems: 1) $Q_{1}$ and $Q_{2}$, 2) $Q_{3}$ and $Q_{4}$ 3) $Q_{5}$, and 4) $Q_{6}$. The performance metrics for the whole system are approximated with the analytical expressions derived from the subsystems. The accuracy of the approximation is validated through simulations in Section VI.

III-A Subsystems 1 and 2: Two queues in tandem

The two queues in tandem, $Q_{1}$ and $Q_{2}$, are considered as one subsystem. The arrival rate for $Q_{1}$ is: $\lambda_{1}=p\alpha$. The Markov chain $\left\{(X_{n}\text{, }Y_{n})\right\}$ is described by

[TABLE]

where $X_{n}$ and $Y_{n}$ denote the states (in terms of queue length) of $Q_{1}$ and $Q_{2}$ at the $n$-th time slot, respectively; $i$ and $j$ are referred to as the level $i$ and phase $j$, respectively. In order to facilitate the presentation, we first analyze a simple example with buffer size $M_{1}=M_{2}=2$. The Markov chain of this example is shown in Fig. 2. However, the analysis presented below is quite general and independent of the specific buffer size. The Markov chain is a Quasi-Birth-and-Death (QBD) DTMC [33]. Note that the QBD only goes a maximum level up or down, the transition matrix has a block partitioned form:

[TABLE]

For simplicity, given a probability of an event, denoted by $p$, we denote the probability of its complementary event by $\bar{p}\triangleq 1-p$. The block matrices of $\mathbf{P}_{1}$ are shown below

[TABLE]

[TABLE]

We construct the Markov chain of this particular example. Let us define $b_{1}=\overline{\lambda}_{1}\bar{\mu}_{2}$, $b_{2}=\overline{\lambda}_{1}\mu_{2}$, $c_{1}=\lambda_{1}\bar{\mu}_{2}$, $c_{2}=\lambda_{1}\mu_{2}$, $a_{1}^{(0)}=\lambda_{1}\bar{\mu}_{1}\bar{\mu}_{2}$, $a_{2}^{(0)}=\lambda_{1}\bar{\mu}_{1}\mu_{2}$, $a_{0}^{(1)}=\lambda_{1}\mu_{1}\bar{\mu}_{2}$, $a_{1}^{(1)}=\lambda_{1}\mu_{1}\mu_{2}+\overline{\lambda}_{1}\bar{\mu}_{1}\bar{\mu}_{2}$, $a_{2}^{(1)}=\overline{\lambda}_{1}\bar{\mu}_{1}\mu_{2}$, and $a_{1}^{(2)}=\overline{\lambda}_{1}\mu_{1}\bar{\mu}_{2}$, $a_{2}^{(2)}=\overline{\lambda}_{1}\mu_{1}\mu_{2}$. Then, we construct the transition matrix of the Markov chain and observe a particular structure of the matrix because of the properties of a QBD DTMC.

By utilizing the properties of a QBD DTMC, we can analyze such systems with arbitrary buffer size in the following steps. First we define the following matrices [33]

[TABLE]

Then, the block matrices of the transition matrix are calculated as [33]

[TABLE]

[TABLE]

Following the steps described above, we can construct the transition matrix of Subsystem $1$ for arbitrary finite buffer sizes. Our goal is to derive the steady state distribution of the Markov chain defined above. We can apply direct methods in order to find the steady state distribution [33, Chapter 4]. Note that there are several efficient algorithms that can be used for this purpose, e.g., logarithmic reduction method [33]. We denote the steady state distribution of Subsystem $1$ by a row vector defined as $\bm{\pi}^{(1)}=\left[\pi_{0,0}^{(1)},\pi_{0,1}^{(1)},\ldots,\pi_{0,M_{2}}^{(1)},\pi_{1,0}^{(1)},\ldots,\pi_{M1,M2}^{(1)}\right]\text{.}$ We find $\bm{\pi}^{(1)}$ by solving the following linear system of equations $\bm{\pi}^{(1)}\mathbf{P}_{1}=\bm{\pi}^{(1)}\text{, }\bm{\pi}^{(1)}\bm{1}=1\text{,}$ where $\bm{1}$ denotes the column vector of ones. Hereafter we use $\bm{\mathbf{\pi}}^{(n)}$ to denote the steady state distribution vector of the $n$-th subsystem for $n=1,2,3,4$.

Furthermore, the arrival rate of $Q_{2}$ depends on the service rate of $Q_{1}$. However, the arrival rate of $Q_{2}$ is equal to $\mu_{1}$ if and only if $Q_{1}$ is non-empty. Therefore, we define the arrival rate of $Q_{2}$ as $\lambda_{2}=\text{Pr}\left\{Q_{1}>0\right\}\mu_{1}=\left(\sum_{j=0}^{M_{2}}\sum_{i=1}^{M_{1}}\pi_{i,j}^{(1)}\right)\mu_{1}\text{.}$ Similarly, we can construct the transition matrix $\mathbf{P}_{2}$ and the steady state distribution $\bm{\pi}^{(2)}$ for the second subsystem consisting of $Q_{3}$ and $Q_{4}$. The arrival rates of $Q_{3}$ and $Q_{4}$ are

[TABLE]

III-B Subsystem 3: $Q_{5}$

We consider $Q_{5}$ as an independent subsystem. We denote by $M_{5}$ the buffer size of the queue. We first define the arrival rate of $Q_{5}$ as $\lambda_{5}=\text{Pr}\left\{Q_{4}>0\right\}\mu_{4}=\left(\sum\limits_{i=0}^{M_{3}}\sum\limits_{j=1}^{M_{4}}\pi_{i,j}^{(2)}\right)\mu_{4}\text{.}$ The Markov chain of this system is shown in Fig. 3.

The transition matrix of this Markov chain is described below

[TABLE]

We denote the steady state distribution of Subsystem 3 by $\bm{\pi}^{(3)}=\left[\pi_{0}^{(3)}\text{, }\pi_{1}^{(3)}\text{,}\ldots\text{, }\pi_{M_{5}}^{(3)}\right]$. To derive $\mathbf{\pi}^{(3)}$, we solve the following linear system of equations, $\bm{\pi}^{(3)}\mathbf{P}_{3}=\bm{\pi}^{(3)}\text{, }\bm{\pi}^{(3)}\mathbf{1}=1\text{.}$

Using balance equation, we obtain

[TABLE]

III-C Subsystem $4$: $Q_{6}$ with finite buffer size

The arrival rate for $Q_{6}$ depends on the service rate of $Q_{2}$ and $Q_{5}$, and the probability the queues to be non-empty. Note that the departures from $Q_{2}$ and $Q_{5}$ can be considered independent stochastic processes. The arrival rates for $Q_{6}$ that occur due to $Q_{2}$ and $Q_{5}$ are $\lambda_{6,2}=\text{Pr}\left\{Q_{2}>0\right\}\mu_{2}\text{, }\text{and }\lambda_{6,5}=\text{Pr}\left\{Q_{5}>0\right\}\mu_{5}\text{, respectively.}$ The arrival rate of $Q_{6}$ is given by $\lambda_{6}=\lambda_{6,2}+\lambda_{6,5}\text{.}$ We model the system as a Markov chain as shown in Fig. 4, where

[TABLE]

The transition matrix that describes the Markov chain above is shown below

[TABLE]

In order to construct the transition matrix $\mathbf{P}_{6}$, we need only to calculate the probabilities described above. We denote the steady state distribution of Subsystem $4$ by $\bm{\pi}^{(4)}=\left[\pi_{0}^{(4)}\text{, }\pi_{1}^{(4)}\text{, }\ldots\right]$.

To derive $\bm{\pi}^{(4)}$, we solve the following linear system of equations

[TABLE]

We observe from (45) that $\bm{\pi}^{(4)}$ is the eigenvector of the transition matrix for $\lambda=1$. Therefore, we can use eigenvalue decomposition (EVD) in order to find the eigenvectors of the matrix. After applying EVD, the transition matrix can be written as $\mathbf{P}^{(4)}=\mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^{-1}\text{,}$ where $\mathbf{Q}=\left[\mathbf{q}_{1},\ldots\mathbf{q}_{n}\right]$ are the eigenvectors, and matrix $\mathbf{\Lambda}$ contains the eigenvalues in its diagonal. The eigenvector that corresponds to the the eigenvalue that has value equal to one, is the steady state distribution.

III-D Subsystem 4: $Q_{6}$ with infinite buffer size

In this subsection, we study the case where the buffer in the core server has infinite size. Since the capabilities of the server in the core can be much higher than those of the MEC server and therefore, the buffer has a very large capacity, we are interested in studying the stability conditions of such systems. We model the system as one Markov chain. The transition matrix that describes the Markov chain is shown below

[TABLE]

where $a_{0}=p_{00}$, $a_{1}=p_{01}$, $a_{2}=p_{02}$. The transition matrix is a lower Hessenberg matrix. The general expression for the equilibrium equations for the states is given by the $i$-th term in the following equation: $\pi_{i}^{(4)}=a_{i}\pi_{0}^{(4)}+\sum\limits_{j=1}^{i+1}b_{i-j}\pi_{j}^{(4)}\text{.}$ For the DTMC with infinite state space, we apply $z$-tansform approach to solve the state equations. The $z$-transform for the state transition probabilities $a_{i}$ and $b_{i}$ are $A(z)=\sum\limits_{i=0}^{2}a_{i}z^{-i}\text{ and }B(z)=\sum\limits_{i=0}^{3}b_{i}z^{-i}\text{,}$ respectively. The $z$-transform for the steady state distribution vector $\bm{\pi}^{(4)}$ is $\Pi(z)=\sum\limits_{i=0}^{\infty}\pi_{i}^{(4)}z^{-i}=\pi_{0}^{(4)}\frac{z^{-1}A(z)-B(z)}{z^{-1}-B(z)}\text{.}$ The solution for $\pi_{i}^{(4)}$ is given by $\pi_{0}^{(4)}=\frac{1+B^{\prime}(1)}{1+B^{\prime}(1)-A^{\prime}(1)}\text{, }\pi_{i}^{(4)}=c_{i}+\sum\limits_{j=1}^{m}r_{j}(p_{j})^{(i-1)}\text{, }i>0\text{,}$ where $r$, $p$, and $c$ are the residues, poles, and directs terms, respectively. Since $Q_{6}$ has infinite buffer size, we need to characterize the conditions under which the queue is stable. Stability is important since it implies finite queueing delay. A definition of queue stability is shown below [34].

Definition 1**.**

Denote by $Q_{i}(t)$ the length of queue $i$ at the beginning of time slot $t$. The queue is said to be stable if $\lim\limits_{t\rightarrow\infty}\Pr\left\{Q_{i}(t)<x\right\}=F(x)$ and $\lim\limits_{x\rightarrow\infty}F(x)=1$.

The corollary consequence of the previous definition is the Loynes’ theorem [35] that states: if the arrival and service processes of a queue are strictly jointly stationary and the average arrival rate is less than the average service rate, then the queue is stable. Therefore, $Q_{6}$ is stable if and only if the following inequality holds: $\lambda_{6}<\mu_{6}\text{.}$

Remark**.**

Since the stationary distributions of the presented Markov chains have been calculated, we can also obtain the distributions of the queue sizes. Thus, we can write the probability that a queue size goes beyond a congestion limit. For example, consider $C$ as the congestion limit for $Q_{5}$, the congestion violation probability is calculated as $\Pr\left\{Q_{5}>C\right\}=\sum\limits_{i=C+1}^{M_{5}}\pi^{(3)}_{i}$.

IV Key Performance Metrics

In this section, we provide analytical expressions of the performance metrics of our interests. First, we calculate the system drop rate and average number of tasks in the system that are necessary metrics for analyzing the throughput and delay of the system. We utilize the results of the previous section in order to obtain the corresponding expressions.

IV-A Drop rate and average queue length

The probabilities to have a dropped task at each time slot for $Q_{1}-Q_{5}$ are shown respectively in below

[TABLE]

[TABLE]

For the case of $Q_{6}$, drops can occur when $Q_{6}$ is on $(M_{6})^{\text{th}}$ or $(M_{6}-1)^{\text{th}}$ state. The drop rate of $Q_{6}$ is shown below

[TABLE]

The average length of each queue is given by

[TABLE]

Therefore, the system drop rate and the average number of tasks in the system can be described as $P_{D}=\sum\limits_{i\in\mathcal{K}}P_{D_{i}}\text{ and }\bar{Q}=\sum\limits_{i\in\mathcal{K}}\bar{Q}_{i}\text{,}$respectively.

IV-B Delay and throughput analysis

We denote by $T_{i}$ the throughput of each queue $i$. Since we consider finite buffers, and consequently packet drops, the throughput of each queue as well as system throughput depend on both the arrival and drop rate. Therefore, the throughput for $Q_{1}-Q_{6}$ is calculated as

[TABLE]

that is equivalent to the effective arrival rate of each queue.

In order to derive the per packet average delay expression, we first derive the expression of packet average delay for each queue that includes both the queueing and transmission delay. We denote by $D_{i}$, the average per packet delay for each queue. We utilize the Little’s theorem [36] and we derive the corresponding expressions as shown below

[TABLE]

Finally, the per packet average delay is calculated as

[TABLE]

V Scaled-Up System

In this section, we explain in details how the analysis that is shown above can be used in general cases of multiple user devices communicating to multiple base stations that are connected to the core server. We provide the analysis and methodology for two base stations. The same methodology can be applied for a system with arbitrary number of base stations. However, we choose the case of two base stations in order to facilitate the presentation. We analyze the system that is shown in Fig. 5, by following the same methodology in order to derive the approximate analytical expressions: we decompose the system into subsystems and analyze each subsystem independently. Furthermore, we consider that the server in the core has two CPUs. Therefore, the number of departures is up to two tasks per time slot. The number of arrivals can be up to four tasks per time slots. We decompose the system into seven individual subsystems: 1) $Q_{1}$ and $Q_{2}$, 2) $Q_{3}$ and $Q_{4}$, 3) $Q_{5}$, 4) $Q_{6}$, 5) $Q_{7}$ and $Q_{8}$, 6) $Q_{9}$ and $Q_{10}$, and 7) $Q_{11}$. For the subsystems $1-3$ and $5-7$, we apply the analysis that is shown in Sections III-IV. For Subsystem 4, we provide the structure of the transition matrix in below. First, we define the arrival rate of $Q_{6}$ that depends on the throughput of $Q_{2}$, $Q_{5}$, $Q_{8}$, and $Q_{11}$. The arrival rates of $Q_{6}$ that occur due to $Q_{2}$, $Q_{5}$, $Q_{8}$, and $Q_{11}$ are $\lambda_{6,2}=\text{Pr}\left\{Q_{2}>0\right\}\mu_{2}\text{, }$ $\lambda_{6,5}=\text{Pr}\left\{Q_{5}>0\right\}\mu_{5}\text{,}$ $\lambda_{6,8}=\text{Pr}\left\{Q_{8}>0\right\}\mu_{8}\text{, and }\lambda_{6,11}=\text{Pr}\left\{Q_{11}>0\right\}\mu_{11}\text{,}$ respectively. Furthermore, we assume that the service time are identically geometrically distributed. Therefore, $\mu_{6,1}=\mu_{6,2}$ and we set $\mu_{6,1}=\mu_{6,2}=\mu$. Let $X$ be a random variable that represents the number of departures at each time slot. Given that the queue is non-empty, $X$ follows a binomial distribution with success probability $\mu$. The Probability Mass Function (PMF) is given by $\text{Pr}\left\{X=k\right\}={n\choose k}\mu^{k}(1-\mu)^{n-k}\text{,}$ where $n$ is the number of cores, and $k$ is the number of departures in one slot. In our case, $n=2$, and $k\leq 2$. Therefore, we can easily calculate the probabilities to have $k$ departures at each time slot. We model this subsystem as one Markov chain. The corresponding transition matrix is shown below

[TABLE]

where the elements of the matrix are calculated in Appendix A.

Let $i$ and $j$ be the $i$-th row and $j$-th column of the transition matrix in (67), respectively. We observe that the transition matrix has repeated elements from the third row to the $(M_{6}-A_{\text{max}}+1)^{\text{th}}$ row, where $A_{\text{max}}$ is the maximum number of tasks that can arrive in $Q_{6}$ at each time slot. Therefore for the system above, $A_{\text{max}}=4$. Thus, after calculating the transition probabilities for the first three rows, it is easy to construct the transition matrix without calculating each element independently.

In order to derive analytical expressions for key performance metrics, we need to calculate the steady state distribution of the Markov chain of $Q_{6}$. We apply EVD to derive the steady state distribution and follow the same methodology as in Section III.

After obtaining the steady state distribution, we can calculate the drop rate of $Q_{6}$. In this case, at each time slot, more than one task can be dropped. The reason is because more than one task can arrive at each time slot. Drops occur when $Q_{6}$ has less than four vacant positions. This corresponds to the $(M_{6}-3)^{\text{th}}$ state. To summarize, drops occur when $Q_{6}$ is on the $(M_{6}-3)^{\text{th}}$, $(M_{6}-2)^{\text{th}}$, $(M_{6}-1)^{\text{th}}$, or $(M_{6})^{\text{th}}$ state. Below, we calculate the drop rate given that $Q_{6}$ is on a particular state

[TABLE]

Finally, the total drop rate of $Q_{6}$ is shown below

[TABLE]

We derive the throughput and delay for each queue by using the equations in (54), (55). The average delay that arrives in the first and the second base station are described below

[TABLE]

respectively. Finally, the per packet average delay is calculated as

[TABLE]

VI Numerical and Simulation Results

In this section, we evaluate the performance of our approximated model by comparing the analytical and simulation results. First, we study the performance of the system with one base station by showing how different parameters and routing decisions affect the system performance in terms of delay, throughput, and drop rate. Second, we study the performance of a system with two base stations. Our goal is to evaluate the performance of our approximation model for larger systems and provide insights on the operation of the considered setup. We developed a MATLAB-based behavioural simulator and each case run for $10^{6}$ timeslots.

VI-A The one base station case

VI-A1 Effect of the routing decision on the system performance - analytical vs simulation results

In Fig. 6 and Fig. 7, we provide results that show how different routing decisions affect the performance of the system for different settings of parameters $\mu_{i},i=1,2,\ldots,5$ . We also compare the simulation results with those of analytical in order to evaluate the performance of our approximated model. In the first case depicted in Fig. 6a, the capabilities of the MEC Server 1 in terms of service rates are lower than those of the MEC server 2. We observe that the optimal values of throughput and delay are achieved for values of $\alpha$ that are close to $0.2$. The value of $\alpha$ that minimizes delay, it does not necessarily maximize the throughput. For larger values of drop rates, we may have shorter delay of the system but smaller throughput as well. In Fig. 6, we observe results for smaller capabilities of the MEC server $2$. In Fig. 7, the capabilities of the MEC servers are identical. We observe that the maximum throughput is achieved for $\alpha=0.5$. However, the optimal $\alpha$ in terms of end-to-end delay is different. The value of $\alpha$ that minimizes the delay is around $0.65$. The flow controller forwards larger amount of traffic to the first flow. This strategy minimizes the delay because the tasks have to traverse a smaller amount of queues and therefore, face shorter waiting times. However, this has an impact on the system performance in terms of drop rate and therefore, the system throughput.

VI-A2 Effect of $\mu_{1}$ and $\mu_{3}$ on the system throughput systems with small buffer size

In this subsection, we provide results for the performance of the system in terms of throughput. In this case, our objective is to maximize the throughput in small buffered systems. In Fig. 8, we provide results that show the optimal routing decisions and the corresponding achieved throughput. The service rates of $Q_{2}$, $Q_{4}$ and $Q_{5}$ are fixed. We find the optimal routing decisions by applying exhaustive search for different values of $\mu_{1}$ and $\mu_{3}$. In Fig. 8a, we show the optimal routing decisions for different values of $\mu_{3}$ and $\mu_{1}$. We observe that when $\mu_{1}<\mu_{3}$, the optimal value of $\alpha$ is close to $0.18$. Therefore, almost $80\%$ of the traffic is processed by the MEC server $2$. On the other hand, when $\mu_{1}>\mu_{3}$, the largest part of the traffic is handled by the MEC server $1$. For the cases of $\mu_{1}=\mu_{3}$, the traffic is splitted between the two flows and the router achieves a balance between them. In Fig. 8b, we provide results that show the optimal achieved throughput for different values of $\mu_{1}$ and $\mu_{3}$. Note that the maximum achievable throughput is equal to the arrival rate, i.e., $0.8$. However, the maximum throughput is not achieved for this setup because of limited buffer capacity. We observe that even for large values of the service rates, the throughput is close to the arrival rate. This indicates that buffers with higher capacities are required in order to achieve a higher maximum throughput.

VI-A3 Effect of $\mu_{1}$ and $\mu_{3}$ on the system delay in systems with large buffer size

In this subsection, we provide results that show the performance of the system in terms of delay. In this case, we minimize the delay in systems with large buffer size by selecting the optimal values of $\alpha$ as shown in Fig. 9. In addition, the optimal values of $\alpha$ are depicted in Fig. 9a. We observe that when $\mu_{1}=\mu_{3}$, the optimal routing decision is to split the traffic into the two flows equally, i.e., $\alpha=0.5$. Another interesting observation is that for even small values of $\mu_{1}$ and large values of $\mu_{3}$, the router decides to split the traffic between the two flows. There are two reasons that explain this phenomenon. The first is that even if the capabilities of one MEC server are low, they should be still utilized in order to increase the performance. The second is that, for the case of minimization of the delay, the router decides to send part of the traffic to the first flow because of the less number of queues facing shorter delay.

VI-A4 Throughput - Delay - Drop rate trade-off

In this subsection, we provide results that show the performance region of the system in terms of throughput, delay, and system drop rate. In Fig. 10a, the performance region is shown. We generate the region by selecting different parameters of the system. In particular, we vary the service rates and the capacities of $Q_{1}-Q_{5}$, and we take the corresponding points as shown in the diagram. Each horizontal line is created by fixing the service rates and varying the capacity of the buffers. The color represents the system drop rate for each different case. We observe that the capacity of the buffers does not affect the system performance in terms of throughput and drop rate, when the service rates are small. On the other hand, small variations of the service rates have a significant impact on the performance metrics. In Fig. 10b, part of the region is provided. For specific requirements of the system, for example, low latency, high reliability, or high throughput, we can provide the proper setup by utilizing the information given by the performance region figures.

VI-B Scaled-up system

In this subsection, we study the performance of a scaled-up system with two base stations, BS1 and BS2. In addition, we evaluate the performance of the proposed approximation model by comparing the analytical and simulations results. The routing decisions for all the cases are equal to $0.5$.

In Fig. 11, the results show the trade-off between system throughput and system delay as the traffic of both base stations increases. It is easy to see that the analytical results are close to the simulation ones which validates the accuracy of the proposed model. In Fig. 12, we provide results that show the performance of each base station for different traffic volumes of BS1. In particular, in Fig. 12a, we observe the trade-off between the throughput and delay of BS1. We see that as the arrival rate increases, the throughput of BS1 increases. However, the throughput of BS1 is not equal to the arrival rate for large values of the arrival rate because of the packet drops.

In Fig. 12b, we provide results that show the trade-off between throughput and delay of BS2 as the arrival rate of BS1 increases. Our goal is to show how the traffic volume of BS1 affects the performance of BS2. It is clear that when the arrival rate of BS1 is larger than the arrival rate of BS2, the performance of BS2 decreases significantly.

In Fig. 13, we provide results that show how the traffic volume of BS2 affects the traffic of BS1. In this case, the service rates are higher than the previous case of Fig. 12. We observe that under this setup, the one system affects the other only in terms of delay. The maximum throughput of BS1 (equal to the arrival rate) is achieved for all the values of $p_{2}$. However, the effect on delay is more significant.

VII Summary

In this work, we consider an exemplary network topology with two MEC servers, a high-end server at core network, and VNF chains embedded in the servers. We model the network and provide a theoretical study on the system performance in terms of system drop rate and average number of the tasks in the system that can be useful for more general set-ups. We provide both experimental and theoretical results in order to evaluate the performance of our approximated model and as it is shown our derived model can approximate the system with high accuracy. Numerical results show that, we are able to offer some useful insights on the design of such systems or resource allocation at each server. Furthermore, we investigate numerically the routing policy that minimizes the system drop rate for different set-ups of the system. This work can be considered as an initial, but significant, step for analyzing and optimizing end-to-end delay, and throughput of such networks. The developed analysis, can provide guidelines for delay-aware routing and resource allocation schemes in similar systems.

Acknowledgment

The authors would like to thank Manuel Stein for numerous fruitful discussions and his valuable suggestions.

Appendix A

We model the queue in the core server of the system as a Markov chain. The Markov chain is depicted In Fig. 14b. In order to facilitate the presentation, we introduce a notation which allows us to design the transitions between states in a more compact way. Below, we calculate the transition probabilities of the Markov chain.

[TABLE]

[TABLE]