---
abstract: 'The radio access network (RAN) is regarded as one of the potential proposals for massive Internet of Things (mIoT), where the random access channel (RACH) procedure should be exploited for IoT devices to access to the RAN. However, modelling of the dynamic process of RACH of mIoT devices is challenging. To address this challenge, we first revisit the frame and minislot structure of the RAN. Then, we correlate the RACH request of an IoT device with its queue status and analyze the queue evolution process. Based on the analysis result, we derive the closed-form expression of the RA success probability of the device. Besides, considering the agreement on converging different services onto a shared infrastructure, we further investigate the RAN slicing for mIoT and bursty ultra-reliable and low latency communications (URLLC) service multiplexing. Specifically, we formulate the RAN slicing problem as an optimization one aiming at optimally orchestrating RAN resources for mIoT slices and bursty URLLC slices to maximize the RA success probability and energy-efficiently satisfy bursty URLLC slices’ quality-of-service (QoS) requirements. A slice resource optimization (SRO) algorithm exploiting relaxation and approximation with provable tightness and error bound is then proposed to mitigate the optimization problem.'
author:
- 'Peng Yang, Xing Xi, Tony Q. S. Quek, , Jingxuan Chen, Xianbin Cao, , Dapeng Wu, [^1]'
bibliography:
- 'Network\_slicing.bib'
title: 'RAN Slicing for Massive IoT and Bursty URLLC Service Multiplexing: Analysis and Optimization'
---

Massive IoT, random access channel, bursty URLLC, RAN slicing

Introduction
============

the explosive growth of the Internet of Things (IoT), massive IoT (mIoT) devices, the number of which is predicted to reach 20.8 billion by 2020, will access to the wireless networks for implementing advanced applications, such as e-health, public safety, smart traffic, virtual navigation/management, remote maintenance and control, and environment monitoring. To address the IoT market, the third-generation partnership project (3GPP) has identified mIoT as one of the three main use cases of 5G and has already initiated several task groups to standardize several solutions including extended coverage GSM (EC-GSM), LTE for machine-type communication (LTE-M), and narrowband IoT (NB-IoT) [@Cellular-Ericsson; @LTE-Nokia].

For establishing massive connections among the wireless networks and mIoT devices, the investigation of reliable and efficient access mechanisms should be prioritized. In accomplishing the massive connections, when an active IoT device wants to transmit signal in the uplink, it randomly chooses a random access (RA) preamble from an RA preamble pool and transmits it through an RA channel (RACH). If more than one device tries to access to a base station (BS) simultaneously, then interference occurs at the RRH. During the past few years, a rich body of works on RA mechanisms has been developed [@xing2019novel; @jang2016non; @jiang2019reinforcement; @zhang2019dnn; @jiang2018analyzing; @grau2019preamble; @jang2019recursive; @jin2017recursive; @jang2019versatile] to mitigate interference and improve the RA success probability or reduce the access delay of an IoT device.

Most of the studies [@xing2019novel; @jang2016non; @jiang2019reinforcement; @zhang2019dnn; @jiang2018analyzing; @grau2019preamble; @jang2019recursive; @jin2017recursive; @jang2019versatile], however, assumed that the whole network resources were reserved for the IoT service and did not investigate the case of the coexistence of IoT service and many other services such as enhanced mobile broadband (eMBB) and ultra-reliable and low latency communications (URLLC). The research of the coexistence of IoT service and other services is essential as future networks are convinced to converge variety of services with different latency, reliability, and throughput requirements onto a shared physical infrastructure rather than deploying individual network solution for each service [@alliance20155g]. What is more, owing to the shared characteristic of network resources, some conclusions obtained in the case of providing sole IoT service may become inapplicable if multiple types of services are required to be supported by the networks.

Network slicing is considered as a promising technology in future networks for providing scalability and flexibility in allocating network resources to various services. Recently, many network slicing frameworks have been developed to provide performance guarantees to IoT or massive machine-type communications (mMTC) service, eMBB service, and URLLC service [@popovski20185g; @budhiraja2019tactile; @ksentini2017toward; @rost2017network; @wen2018robustness; @buyakar2018poster; @guo2019enabling].

Different from previous works, this paper investigates the mIoT and bursty URLLC service multiplexing via slicing the radio access network (RAN). This study is highly challenging because i) performance requirements of a massive number of IoT devices should be satisfied. Yet, the typical 5G cellular IoT, NB-IoT can admit only 50,000 devices per cell [@3GPP15Cellular]; ii) RAN slicing operation (e.g., creating, activating, and releasing slices) has to be conducted in a timescale of minutes to hours to keep pace with the upper layer slicing. However, the wireless channel generally changes in a timescale of millisecond to seconds. Results of the RAN slicing operation are desired to be achieved based on the time-varying channel. Thus, the RAN slicing should tackle a two timescale issue [@tang2019service].

These challenges motivate us to investigate the RAN slicing for mIoT and bursty URLLC service provision to maximize the utility of mIoT slices and that of bursty URLLC slices. The main contributions of this paper can be summarized as the following:

-   We revisit the frame and minislot structure for mIoT transmission to accommodate more RA requests from a massive number of IoT devices.

-   We adopt a queueing model to track the IoT packet arrival, accumulate and departure processes and analyze the queue evolution process by employing probability and stochastic geometry theories. Based on the analysis result, we derive the closed-form expression of the RA success probability of a randomly chosen IoT device.

-   We define mIoT slice utility and bursty URLLC slice utility and formulate the RAN slicing for mIoT and bursty URLLC service multiplexing as a resource optimization problem. The objective of the optimization problem is to maximize the total mIoT and URLLC slice utilities, subject to limited physical resource constraints. The mitigation of this problem is difficult due to the existence of indeterministic objective function and thorny non-convex constraints and the requirement of tackling a two timescale issue as well.

-   To mitigate this thorny optimization problem, we propose a slice resource optimization (SRO) algorithm. In this algorithm, we first exploit a sample average approximate (SAA) technique and an alternating direction method of multipliers (ADMM) to tackle the indeterministic objective function and the two timescale issue. Then, a semidefinite relaxation (SDR) scheme joint with a Taylor expansion scheme are leveraged to approximate the non-convex problem as a convex one. The tightness of the SDR scheme and the error bound of the Taylor expansion are also analyzed.

The remaining of this paper is organized as follows. We review the prior arts in Section II. In Section III, we describe our system model and formulate the service multiplexing problem in Section IV. The problem-mitigating algorithm is presented in Sections V and VI. In Section VII, we give the simulation results and conclude this paper in Section VIII.

Prior arts
==========

Recently, many researches have been conducted to increase RA success probabilities and/or reduce the access delay of mIoT devices. They can be generally classified into two groups: traffic detection and estimation based algorithms and algorithms without traffic detection and estimation.

The fundamental idea of the traffic detection and estimation based algorithms is to design an RA algorithm based on the detected and/or estimated users’ activity and traffic congestion situation and so on. For example, to reduce the access delay, a grant-free non-orthogonal RA system relying on the accurate user activity detection and channel estimation was proposed in [@zhang2019dnn]. A traffic-aware spatiotemporal model for the contention-based RA analysis is conducted for mIoT networks in [@jiang2018analyzing]. With the spatiotemporal model, a hybrid power ramping and back-off RA scheme was then developed to improve the RA success probability. Besides, an extended pseudo-Bayesian backlog estimation scheme was exploited in [@jang2019versatile] to estimate the number of backlogged nodes to attempt access. A versatile access control mechanism was then designed to reduce the access delay based on the estimation results.

For algorithms without detection and estimation, they design RA schemes without detecting users’ activity or estimating the statistical characteristic of traffic. For instance, the work in [@xing2019novel] proposed to improve the RA success probability of an IoT device by exploiting a distributed queue mechanism and then proposed an access resource grouping mechanism to reduce the access delay caused by the queuing process of the distributed queue mechanism. To increase RA success probability, the work in [@jang2016non] proposed to increase the number of preambles at the first step of the RA procedure by utilizing a spatial group mechanism and improve resource utilization through non-orthogonally allocating uplink channel resources at the second step of the RA procedure. Additionally, without knowing the statistical characteristic of traffic, a reinforcement learning-based algorithm was proposed in [@jiang2019reinforcement] to determine the uplink resource configuration for RA such that the average number of served IoT devices was maximized while ensuring a high RA success probability.

Except for the IoT service, future networks are envisioned to simultaneously support different services and applications with significantly different requirements on reliability, latency and bandwidth. As a result, researchers are now paying more attention to the service multiplexing of IoT/mMTC and many other services such as eMBB and URLLC. For example, instead of slicing the RAN via orthogonal resource allocation among different services, the work in [@popovski20185g; @budhiraja2019tactile] studied the potential advantages of allowing for non-orthogonal RAN resources sharing in uplink communications from a collection of mMTC, eMBB, and URLLC devices to the same BS. The work in [@ksentini2017toward] developed a two-level scheduling process to allocate dynamically dedicated bandwidth to each network slice according to workload demand and slices’ quality of service (QoS) requirement such that flexible resource allocation could be implemented. The work in [@rost2017network] proposed to maintain slice-specific radio resource control elements with which the RAN protocol stacks and different slices were configured. Besides, the work in [@wen2018robustness] aimed to optimize the virtual network functions and infrastructure resources such as the system bandwidth to implement slice recovery and reconfiguration for mMTC and eMBB service provision. The work in [@buyakar2018poster] proposed to maintain slice isolation between mMTC and eMBB slices and meet the performance requirements of these slices through limiting and dynamically updating the amount of resources allocated to each slice and monitoring the resource usage of each slice. After representing the slice performance requirements as the required amount of resources per deadline interval, an idea of earliest-deadline and first-scheduling was exploited in [@guo2019enabling] to allocate radio resources to mMTC, eMBB, and URLLC slices effectively.

System model
============

We consider a coordinated-multipoint-enabled RAN slicing system for mIoT and bursty URLLC multiplexing service provision. From the viewpoint of the infrastructure composition, the system mainly includes one baseband unit (BBU) pool and multiple remote radio heads (RRHs) that connect to the BBU via fronthaul links. From the perspective of network slicing, two types of inter-slices, i.e., mIoT slices and URLLC slices, are exploited in this system with $\mathcal{S}^I$ and $\mathcal{S}^u$ representing the sets of mIoT slices and URLLC slices, respectively. We focus on the modelling of uplink IoT data transmission in mIoT slices and the modelling of downlink URLLC data transmission in URLLC slices. The IoT devices are spatially distributed in ${\mathbb R}^2$ according to an independent homogeneous Poisson point process (PPP) $\Phi_s = \{u_{i,s}; s \in \mathcal{S}^I, i = 1, 2, \ldots\}$ with intensity $\lambda_s^I$, where $u_{i,s}$ denotes the location of the $i$-th IoT device in the $s$-th mIoT slice. There are also $N^u$ URLLC devices that are randomly and evenly distributed in ${\mathbb R}^2$. The RRHs are spatially distributed in ${\mathbb R}^2$ according to an independent PPP $\Phi_R = \{v_j; j = 1, 2, \ldots\}$ with intensity $\lambda_R$, where $v_j$ represents the location of the $j$-th RRH. The number and locations of IoT devices and RRHs will be fixed once deployed. Besides, each RRH is equipped with $K$ antennas, and each device is equipped with a single antenna. In IoT network slices, each IoT device is assumed to connect to its geographically closest RRH [@jiang2018analyzing]; thus, the cell area of each RRH constitutes a Voronoi tessellation. In URLLC network slices, RRHs cooperate to transmit signals to a URLLC device to improve its signal-to-noise ratio (SNR). A flexible frequency division multiple access (FDMA) technique is utilized to achieve the inter-slice and intra-slice interference isolation [@tang2019service].

The system time is discretized and partitioned into time slots and minislots with a time slot consisting of $T$ minislots. On the one hand, at the beginning of each time slot, a RAN slicing coordinator [@How2019yang] will decide whether to accept or reject received network slice requests which will be defined in the following subsections. Once a slice request is accepted, a network slice management will be responsible for activating or creating a virtual slice that is well resource-configured to satisfied the QoS requirements of devices in the slice [@How2019yang]. The slice configuration process is time costly and will generally be conducted in a timescale of minutes to hours [@tang2019service]. On the other hand, at the beginning of each minislot, each active IoT device may try to connect to its associated RRH, and RRHs will generate cooperated beamformers based on sensed channel coefficients.

mIoT slice model
----------------

By referring to the concept of a network slice [@rost2017network], especially from the viewpoint of the QoS requirement of a slice, we can define a mIoT slice request as follows.

\[def:IoT\_slice\_definition\] A mIoT slice request is defined as a tuple $\{\lambda_s^I, \gamma_s^{th}, N_{a,s}\}$ for any slice $s \in {\mathcal S}^I$, where $\gamma_s^{th}$ is the requirement of data transfer rate from an IoT device in $s$ to its associated RRH, $N_{a,s}$ denotes the number of accumulated packets in a queue of an IoT device in $s$.

In this paper, all mIoT slice requests are always accepted by the RAN slicing coordinator. IoT devices with the same data transfer rate are assigned to the similar slice. For an IoT device in $s$, if it has the opportunity to send its endogenous arrival packets to the corresponding RRH, then it will randomly select a preamble (e.g., orthogonal Zadoff¨CChu sequences) from a BBU-maintained preamble pool and transmit the preamble to the RRH at the data rate $\gamma_s^{th}$. Just like the literature [@soorki2017stochastic; @gharbieh2017spatiotemporal], if the RRH can successfully decode the preamble, then a connection between the IoT device and the RRH are considered to be set up although the whole connection establishment process usually follows an RA four-step procedure [@grau2019preamble]. In other word, the RA success probability is regarded as the probability of successfully transmitting a preamble in this paper. Next, we will analyze an IoT device queue evolution model, with the analysis of which the RA success probability of the IoT device will be derived.

### Queue evolution model

The queue evolution process consists of the packet arrival process, packet accumulate process, and packet departure process.

During minislot $t$, a Poisson distribution with intensity (or arrival rate) $\epsilon_{w,s}(t)$ is exploited to model the random, mutually independent endogenous packet arrivals in an IoT device in slice $s$. Then during minislot $t$ with a duration $\tau$, the arrival intensity of new packets can be expressed as $\mu_{w,s}(t) = \epsilon_{w,s}(t)\tau$. Once arrived, new packets will not be sent out immediately in general and will enter a queue, which is modelled as an $M/M/k$ queue with unlimited capacity, to wait for their scheduling. In the $M/M/k$ queue, packets will be scheduled according to the first-come, first-served (FCFS) basis. The unlimited queue capacity indicates that the age of information [@kaul2011minimizing; @kaul2012real] of new arrivals will not be considered, and packets will not be dropped before sending out. Besides, owing to the RA behavior of a slotted-ALOHA protocol, new arrivals during $t$ will only be counted at minislot $t + 1$. Thus, the accumulated number of packets $N_{a,s}(t)$ of a randomly selected IoT device in slice $s$ at $t$ is determined by the accumulated number of packets and the number of new arrivals at $t - 1$ and whether the preamble of the device can be successfully decoded by its associated RRH. Table \[table:queue\_evolution\] shows the evolution of accumulated packets in an IoT device. In this table, $x_s = \gamma_s^{th}/L$ packets at the head of the queue will be popped out if the corresponding RA succeeds, where $L$ denotes the IoT packet length; otherwise, they will be kept in the queue and wait for the opportunity of re-transmission at the next minislot. The operation $[x]^+ = \max(x,0)$.

       Value                              Success                                          Failure
  ---------------- ------------------------------------------------------ -----------------------------------------
    $N_{a,s}(1)$                            $0$                                              $0$
    $N_{a,s}(2)$                 $[N_{w,s}{(1)}$ - $x_s]^+$                            $N_{w,s}{(1)}$
   $N_{a,s}{(3)}$       $[N_{a,s}{(2)}$ + $N_{w,s}{(2)}$ - $x_s]^+$            $N_{a,s}{(2)}$ + $N_{w,s}{(2)}$
      $\ldots$                            $\ldots$                                        $\ldots$
   $N_{a,s}{(t)}$   $ [N_{a,s}(t$ - $1)$ + $N_{w,s}(t$ - $1)$ - $x_s]^+$   $N_{a,s}(t$ - $1)$ + $N_{w,s}(t$ - $1)$

  : Accumulated packets evolution in an IoT device[]{data-label="table:queue_evolution"}

With the evolution of accumulated packets, we can define the non-empty probability of the queue of an IoT device in $s$ as the following.

\[def:non\_empty\_prob\] At minislot $t$, for a randomly selected IoT device in slice $s \in \mathcal{S}^I$, the probability that its queue is not empty can be defined as $$\label{eq:non_empty_prob}
    P_{ne,s}{(t)} = {\mathbb P}\{N_{a,s}(t) > 0\}, \forall s \in {\mathcal{S}^I}$$

(\[eq:non\_empty\_prob\]) implicitly reflects that new arrival packets at $t$ will not be sent out immediately. According to the evolution of $N_{a,s}(t)$, it can be observed that $P_{ne,s}(t)$ is determined by the probability distribution of $N_{a,s}(t-1)$ and the RA success probability. Since these probabilities and their correlations are unknown, the derivation of the explicit expression of $P_{ne,s}(t)$ is difficult.

Next, we describe the packet departure process combined with a frame and minislot structure for mIoT packets transmission. As mentioned above, partly because of the limitation on the frame and minislot structure, NB-IoT and LTE-M can only admit 50,000 devices. For NB-IoT, only one physical resource block (PRB) with a bandwidth of $180$ KHz in the frequency domain is allocated for IoT transmission, and each physical channel occupies the whole PRB. For LTE-M, although the physical channels are time and frequency multiplexed, it only reserves six in-band PRBs with a total bandwidth of $1.08$ MHz in the frequency domain for IoT data transmission. Therefore, the frame and minislot structure for mIoT transmission should be revisited if more RA requests from IoT devices want to be accepted.

Fig. \[fig:fig\_frame\_minislot\_structure\] depicts a frame and minislot structure for mIoT transmission in each mIoT slice[^2]. In this structure, both the frequency division multiplexing scheme and code division multiplexing scheme are leveraged to admit more IoT devices in the way of alleviating the mutual device interference. Particularly, the frequency division multiplexing scheme alleviates signal interference through orthogonal frequency allocation, and the code division multiplexing scheme mitigates the co-channel signal interference via reducing the cross-correlation of simultaneous transmissions. The combination of the two schemes may significantly mitigate interference experienced at an RRH. In this way, the QoS requirements of more IoT devices may be satisfied, and the RAN slicing system may support more IoT devices. For a mIoT slice $s \in \mathcal{S}^I$, each subframe includes $F_s$ orthogonal uplink physical RA channels (PRACHs). A single tone mode with a tone spacing of size of $a$ MHz is adopted for each uplink PRACH, which indicates that each PRACH occupies a PRB. At the beginning of each minislot, an active IoT device, i.e., the device’s queue is non-empty, will randomly choose a preamble from a set of non-dedicated RA preambles of size $\xi$ and transmit the preamble through a randomly selected PRACH. For each preamble, it has an equal probability ${\frac{1}{\xi}}$ to be chosen by each IoT device. Similarly, each PRACH has an equal probability $\frac{1}{F_s}$ to be selected. Thus, the average number of IoT devices in mIoT slice $s \in \mathcal{S}^I$ choosing the same PRACH and the same preamble is $\frac{\lambda_s^I}{\xi F_s}$. Notably, a greater $\xi F_s$ may significantly reduce signal interference experienced at each RRH.

Then, the following question should be tackled: *how many PRBs should be reserved for mIoT transmission?* To improve the resource utilization, the resource allocated to mIoT should be determined according to the requirements of mIoT and other coexistence services. It motivates us to optimize the resources orchestrated for the mIoT service except for analyzing the RACH procedure of IoT devices. The optimization procedure will be discussed in detail in the next section.

![The frame and minislot structure. ’R’ and ’D’ denote the resource block reserved for preamble and IoT data transmission. PBCH, PSS and SSS represent the PRBs for physical broadcast channel, primary synchronization signal and secondary synchronization signal transmission, respectively.[]{data-label="fig:fig_frame_minislot_structure"}](frame_slot_structure.eps){width="2.9in"}

### Access control scheme

In a mIoT network slice, as the slotted-ALOHA protocol allows all active IoT devices to request for RA at the beginning of each minislot without checking the status of channels, IoT devices may simultaneously transmit preambles. It may incur severe slice congestion that may lower the RA success probabilities of IoT devices and degrade the system performance. Access control has been considered as an efficient proposal of alleviating congestion, and many access control schemes have been proposed [@jiang2018analyzing; @Study163GPP]. In this paper, we aim at illustrating the performance difference between a network slicing system without access control and with access control. Therefore, we adopt the following two schemes [@jiang2018analyzing]:

-   **Unrestricted scheme:** each active IoT device requests the RACH at the beginning of minislot $t$ without access restriction. If mIoT slices are not crowded or in a light-crowded condition, then this scheme may quickly flush queues of IoT devices. However, if a heavy-crowded condition is encountered, then this scheme may result in a high packet queueing delay.

-   **Access class barring (ACB) scheme:** at the beginning of $t$, each active IoT device throws a random number $q \in [0, 1]$ and can request the RACH only if $q < P_{ACB}$, where $P_{ACB}$ is an ACB factor determined by RRHs based on the slice congestion condition. The ACB scheme can relieve slice congestion to some extent by reducing RACH requests of active IoT devices.

With the introduced access control schemes, we can define the non-restriction probability of a randomly selected IoT device in $s$ as follows.

\[mydef:unrestricted\_prob\] At minislot $t$, for a randomly selected IoT device in slice $s \in \mathcal{S}^I$, the probability that its RACH request is not restricted is defined as $$\label{eq:non_restriction_prob}
    P_{nr,s}{(t)} = {\mathbb P}\{{\rm Unrestricted} \text{ } {\rm RACH} \text{ } {\rm requests}\}, \forall s \in {\mathcal{S}^I}$$

For all $s \in \mathcal{S}^I$ at any minislot $t$, we have $P_{nr,s}(t) = 1$ for the unrestricted scheme and $P_{nr,s}(t) = P_{ACB}$ for the ACB scheme.

### Analysis of RA success probability

For an RRH, two significant reasons may lead to an error preamble decoding i) the achieved preamble transfer rate at the RRH is less than a preset threshold; ii) the RRH simultaneously decodes at least two similar co-channel preambles, and thus preamble collision occurs. The research of the mitigation of preamble collision has been well conducted in [@jang2019versatile; @jiang2018collision]. Just like [@jiang2018random], we focus on the exploration of enabling successful single preamble transmission that is discussed in detail as follows.

We utilize a power-law path-loss model to calculate the path-loss between an IoT device and its RRH in mIoT slices and utilize a truncated channel inversion power control scheme to eliminate the ’near-far’ effect. In the power-law path-loss model, the IoT device transmit power decays at the rate of $r^{-{\varphi}}$ with $r$ representing the propagation distance and $\varphi$ denoting the path-loss exponent. In the power control scheme, IoT devices associated with the same RRH compensate for the path-loss to maintain the average received signal power at the RRH equal to a threshold $\rho_o$. Without loss of generality, the cutoff threshold $\rho_o$ is set to be the same for all RRHs. Owing to the channel deep fading, severe co-channel interference, and insufficient transmit power, an IoT device may experience uplink preamble transmission outage. The following definition describes the definition of the probability that a randomly selected IoT device can successfully transmit a chosen preamble to its corresponding RRH.

\[mydef:preamble\_trans\_suc\_prob\] At minislot $t$, for a randomly selected active IoT device in slice $s \in \mathcal{S}^I$, its RA success probability is defined as $$\label{eq:preamble_trans_suc_prob}
    P_{s}{(t)} = {\mathbb P}\{r_s(t) \ge \gamma_s^{th}\}, \forall s \in {\mathcal{S}^I}$$ where $r_s(t) = a{\log _2}(1 + SIN{R_s}(t))$ denotes the achieved preamble transfer rate at the IoT device’s associated RRH and $SIN{R_s}(t)$ is the signal-to-interference-plus-noise ratio (SINR).

Then, for any active IoT device in $s$, its QoS requirement is given by $$\label{eq:mMTC_QoS}
    P_s(t)  \ge {\pi _s}, \forall s \in {\mathcal{S}^I}$$ where $\pi_s$ denotes a threshold of the required RA success probability.

This definition shows that the QoS requirement of each active IoT device in $s$ should be satisfied if the slice request of $s$ is accepted. The definition also states that $P_{s}{(t)}$ is correlated with the non-empty probability $P_{ne,s}{(t)}$. Recall that the RA success probability of an IoT device impacts its non-empty probability, we can know that the RA success probability and the non-empty probability are intertwined. Additionally, $SIN{R_s}(t)$ is a function of complicated co-channel interference. Thus, it is hard to obtain the closed-form expression of $P_{s}{(t)}$.

Without any loss in generality, we perform the analysis of RA success probability on an RRH located at the origin. According to Slivnyak’s theorem [@haenggi2012stochastic], the analysis holds for a generic RRH located at a generic location. For a randomly selected IoT device with non-empty queue in $s \in \mathcal{S}^I$, the theoretical preamble transfer rate experienced at the RRH located at the origin can take the form $$\label{eq:preamble_trans_suc_prob}
    r_s(t) = a \log_2 \left (1 + \frac{{\rho_o {h_o}}}{{{\sigma ^2} + {{\cal I}_{s}}(t)}} \right ), \forall s \in {\mathcal{S}^I}$$ where $\sigma^2$ represents the noise power, ${{\cal I}_{s}}(t)$ denotes signal interference received at the RRH, the useful signal power equals to $\rho_o h_o$ due to the truncated channel inversion power control[^3] [@elsawy2014stochastic] with $h_o$ denoting the channel power gain between the IoT device and the RRH. It is noteworthy that the channel power gain experienced at a generic RRH is related to the spatial locations of both the RRH and its associated IoT devices. Nevertheless, we drop the spatial indices for notation lightening. Besides, just like [@elsawy2014stochastic], all channel gains are assumed to be known and be independent of each other, independent of the spatial locations, symmetric and are identically distributed (i.i.d.). Considering both the particular IoT device deployment environment and the convenience of theoretical analysis, the Rayleigh fading is assumed, and the channel power gain $h_o$ is assumed to be exponentially distributed with unit mean.

Based on the following five facts, we next present the analytical expression of signal interference

-   **Fact 1**: the average signal received from any single IoT device belonging to inter-cells is strictly less than $\rho_o$.

-   **Fact 2**: the average interference signal received from any single interfering IoT device associated with the origin RRH strictly equals to $\rho_o$.

-   **Fact 3**: IoT devices choosing the same co-channel preamble as the randomly selected IoT device may become an interfering IoT device.

-   **Fact 4**: at each minislot, IoT devices with non-empty queue may become interfering IoT devices.

-   **Fact 5**: IoT devices in difference slices may not mutually interfere.

Note that Fact 1 and Fact 2 are direct consequences of the device-RRH association policy and power control scheme. Fact 5 holds due to the exploration of intra-slice isolation. Therefore, the aggregate interference received at the origin RRH can take the following form $$\label{eq:I_intra}
\begin{array}{l}
    {\cal I}_{s}(t) = \sum\limits_{{u_{m,s}} \in \Phi_s \backslash \{ o\} } {{\mathbbm 1}({p_m}||{d_m}|{|^{ - \varphi }} = {\rho _o}){\mathbbm  1}(N_{a,s}{(t)} > 0)} \times \\
    \qquad \quad \text{ } {{\mathbbm 1}({f_m} = {f_o}){\rho _o}{h_m}}, \forall s \in {\mathcal{S}^I}
\end{array}$$ where $o$ represents the randomly selected IoT device associated with the RRH at the origin, $p_m$ represents the transmit power of the $m$-th IoT device, $||d_m||$ is the distance between the $m$-th IoT device and the origin RRH, $f_o$ denotes the preamble and channel chosen by the randomly selected IoT device, $f_o = f_m$ indicates that the randomly selected IoT device and the $m$-th IoT device select the same preamble and channel. ${\mathbbm 1}(\cdot)$ is the indicator function that equals to one if the statement ${\mathbbm 1}(\cdot)$ is ture; otherwise, it equals to zero. Just like [@zhang2015resource], in (\[eq:I\_intra\]), co-channel inter-cell interference is assumed as a part of thermal noise mainly because of the severe wall penetration loss.

Then, for the randomly selected IoT device in $s \in {\mathcal{S}^I}$, we can rewrite (\[eq:preamble\_trans\_suc\_prob\]) as the following form with (\[eq:I\_intra\]) $$\label{eq:QoS_analysis}
    \begin{array}{l}
P_s(t) = {\mathbb P}\{ SINR_s(t) \ge \theta _s^{th}\} \\
\qquad \text{ } = {\mathbb P}\{ {h_o} \ge \frac{{\theta _s^{th}}}{{{\rho _o}}}({\sigma ^2} + {{\cal I}_{s}(t)})\} \\
\qquad \text{ } \mathop  = \limits^{(a)} {{\mathbb E}}\left[ {\exp \left\{ { - \frac{{\theta _s^{th}}}{{{\rho _o}}}({\sigma ^2} + {{\cal I}_{s}(t)})} \right\}} \right]\\
\qquad \text{ } =  \exp \left\{ { - \frac{{\theta _s^{th}}}{{{\rho _o}}}{\sigma ^2}} \right\}{\cal L}_{{{\mathcal I}_{s}(t)}}\left(\frac{\theta_s^{th}}{\rho_o}\right), \forall s \in {\mathcal{S}^I}
\end{array}$$ where $\theta_s^{th} = {2^{\gamma _s^{th}/a}} - 1$. (a) follows from the full probability law over ${\cal I}_{s}(t)$, and ${\cal L}_{{{\cal I}_{s}(t)}} (\cdot)$ denotes the Laplace transform (LT) of the probability density function (PDF) of the random variable ${\cal I}_{s}(t)$. Note that the notation ${\cal L}_{{{\cal I}_{s}(t)}} (\cdot)$ is a terminology that is a slight abuse of subscript ${{\cal I}_{s}(t)}$.

The following lemma characterizes the LT of aggregate interference $\mathcal{I}_s(t)$.

\[lem:LT\_interference\_expression\] For the origin RRH, the LT of its received aggregate interference from active IoT devices associated with it is given by $$\label{eq:LT_interference_expression}
        {\cal L}_{{{\cal I}_{s}(t)}} \left(\varpi_s\right) = {\frac{{1 + \varpi_s{\rho _o}}}{{{{\left( {1 + \alpha_s \varpi_s{\rho _o}/\left( {1 + \varpi_s{\rho _o}} \right)} \right)}^{3.5}}}} - \frac{{1 + \varpi_s{\rho _o}}}{{{{\left( {1 + \alpha_s } \right)}^{3.5}}}}}$$ where $\varpi_s = \frac{\theta_s^{th}}{\rho_o}$, $\alpha_s  = \frac{{{P_{nr,s}}(t){P_{ne,s}}(t){\lambda _{{s}}^I}}}{{3.5{\lambda _R}{\xi F_s}}}$, for all $s \in {\mathcal{S}^I}$.

Please refer to Appendix A.

With the conclusion in Lemma \[lem:LT\_interference\_expression\], we can then obtain the mathematical expression of the RA success probability of a randomly selected IoT device at $t$ in the following corollary.

\[lem:QoS\_prob\_analysis\] For a randomly selected IoT device in a mIoT slice $s \in {\cal S}^I$, its RA success probability at minislot $t$ is given by $$\label{eq:QoS_prob_analysis}
       {P_s}(t) =  {\frac{{(1 + \varpi_s{\rho _o})}{{e^{-\varpi_s{\sigma ^2}}}}}{{{{\left( {1 + \alpha_s \varpi_s{\rho _o}/\left( {1 + \varpi_s{\rho _o}} \right)} \right)}^{3.5}}}} - \frac{{(1 + \varpi_s{\rho _o})}{{e^{-\varpi_s{\sigma ^2}}}}}{{{{\left( {1 + \alpha_s } \right)}^{3.5}}}}}$$

By substituting (\[eq:LT\_interference\_expression\]) into (\[eq:QoS\_analysis\]), we can obtain (\[eq:QoS\_prob\_analysis\]).

Although Corollary \[lem:QoS\_prob\_analysis\] presents a mathematical expression of $P_s(t)$, the expression is not in the closed-form as it is a function of $P_{ne,s}(t)$ the closed-form expression of which is not obtained. Next, we derive the closed-form expression of $P_{ne,s}(t)$.

### Analysis of non-empty probability

According to the definition of non-empty probability, $P_{ne,s}(t)$ is correlated with the number of accumulated packets $N_{a,s}(t)$ of the randomly selected IoT device in mIoT slice $s$. Thus, we theoretically analyze the non-empty probability of the randomly selected IoT device as the following.

As the number of the accumulated packets in the queue of a randomly selected IoT device in slice $s \in \mathcal{S}^I$ at the $1^{\rm st}$ minislot is empty, its non-empty probability $P_{ne,s}^{1}$ at the $1^{\rm st}$ minislot can take the form $$\label{eq:non_empty_prob_1}
    P_{ne,s}^{1} = {\mathbb P}\{N_{a,s}^1 > 0\} = 0, \forall s \in {\mathcal{S}^I}$$ where we write $x^t$ instead of $x(t)$ to lighten the notation.

The following lemma presents the closed-form expression of the non-empty probability of a randomly selected IoT device served by the origin RRH when minislot $t > 1$.

\[lem:non\_empty\_prob\_lemma\] The number of accumulated packets of a randomly selected IoT device served by the origin RRH at minislot $t > 1$ may be approximately Poisson distributed. Therefore, for any mIoT slice $s \in \mathcal{S}^I$, we approximate the number of accumulated packets $N_{a,s}^{t}$ at minislot $t$ as a Poisson distribution with intensity $\mu_{a,s}^{t}$, which is given by $$\label{eq:mu_accumulated_packets}
        \mu _{a,s}^t = \left [\mu _{w,s}^{t - 1} + \mu _{a,s}^{t - 1} - P_{s}^{t - 1}\left( {1 - {e^{ - \mu _{w,s}^{t - 1} - \mu _{a,s}^{t - 1}}}} \right) \right ]^+, \forall s \in {\mathcal{S}^I}$$

Then, the non-empty probability of the device at minislot $t$ can be written as $$\label{eq:non_empty_prob_m}
        P_{ne,s}^t = 1 - {e^{ - \mu _{a,s}^t}}, \forall s \in {\mathcal{S}^I}$$

Please refer to Appendix B.

Combine with (\[eq:QoS\_prob\_analysis\]) and (\[eq:mu\_accumulated\_packets\]), the closed-form expression of $P_s(t)$ ($s \in \mathcal{S}^I$) can be obtained.

Bursty URLLC slice model
------------------------

Similar to the definition of a mIoT slice request, a bursty URLLC slice request can be defined as the following.

\[mydef:URLLC\_slice\_request\] A bursty URLLC slice request is defined as four tuples $\{I_s^u, D_s, \alpha, \beta\}$ for slice $s \in {\mathcal S}^u$, where $I_s^u$ denotes the number of URLLC devices in $s$, $D_s$ denotes the transmission latency requirement of each URLLC device in $s$, $\alpha$ and $\beta$ represent the packet blocking probability and the codeword error decoding probability of each URLLC device, respectively.

In this definition, URLLC devices are grouped into $|{\mathcal S}^u|$ clusters according to the transmission latency requirement of each device. Owing to the low latency requirement URLLC packets should be immediately scheduled upon arrival; thus, all URLLC slice requests will always be accepted by the RAN coordinator. Except for the low packet error decoding probability that has been emphasized for URLLC transmission in a plenty of works [@tang2019service], this paper attempts to orchestrate slice resources to reduce the packet blocking probability for bursty URLLC transmission. This is because the bursty characteristic of URLLC traffic [@How2019yang] may lead to the packet blocking in URLLC slices, which may significantly reduce the reliability of URLLC transmission. Therefore, the indicators $\alpha$ and $\beta$ are involved to reflect the ultra-reliable requirement of URLLC transmission jointly.

Then, we address the following question: *how to orchestrate slice resources for reducing packet blocking probability and codeword error decoding probability?*

### Reduction of packet blocking probability

As mentioned above, the bursty feature of URLLC traffic is the crucial factor that leads to the URLLC packet blocking for URLLC transmission. Therefore, we next model the bursty URLLC traffic based on which we discuss how to orchestrate slice resources to alleviate the impact of bursty URLLC traffic.

During minislot $t$, an independent homogeneous Poisson distribution with intensity ${\bm \lambda} = \{\lambda_s; s \in {\mathcal{S}^u}\}$ is utilized to model the number of bursty URLLC packets aggregated at RRHs, where $\lambda_s$ denotes the intensity of new arrivals destined to devices belonging to URLLC slice $s$.

Once arrived, new URLLC arrivals will enter a queue maintained by an RRH to be immediately scheduled. An $M/M/W^u$ queueing system with limited bandwidth $W^u$ is exploited to model the queue. Without loss of any generality, we assume that each RRH maintains the same queue due to the exploration of cooperated transmission. In the queue, a packet destined to URLLC device $i \in {\mathcal{I}_s^u}$, $s \in {\mathcal{S}^u}$ will be allocated with a block of system bandwidth $\omega_{i,s}^u(t)$ for a period of time $d_s \le D_s$ at $t$. Owing to stochastic variations in the bursty packet arrival process, the limited bandwidth may not be enough to serve new arrivals occasionally. As such, URLLC packet blocking may happen. To reduce the probability of URLLC packet blocking, the redesign of URLLC frame and minislot structure may be required.

At minislot $t$, let $P_b(\bm \omega^u(t), \bm \lambda, \bm d, W^u(t))$ denote the packet blocking probability experienced at an RRH, where ${{\bm \omega}^u(t) = \{\omega_{1,1}^u(t), \ldots, \omega_{i,s}^u(t), \ldots, \omega_{I_{|{\mathcal S}^u|}^u,|{\mathcal S}^u|}^u(t)\}}$ and ${\bm d} = \{d_1, \ldots, d_s \ldots, d_{|{\mathcal S}^u|}\}$. The Theorem 1 in [@anand2018resource] provides us with a clue of redesigning the URLLC frame and minislot structure in the time-frequency plane for bursty URLLC traffic transmission. This theorem indicates that if we narrow the PRB of the URLLC frame in the frequency domain while widening it in the time domain, then the number of concurrent transmissions will be increased. As a result, the packet blocking probability is reduced.

Therefore, for a URLLC packet destined to device $i \in {\mathcal I}_s^u$, $s \in \mathcal{S}^u$, we should scale up $d_s$ and choose $d_s$ and $\omega_{i,s}^u(t)$ at $t$ using the following equation $$\label{eq:omega_s}
   d_{{s}}  = D_s \ {\rm and} \ \omega_{{i,s}}^{u}(t) = \frac{b_{i,s}^u(t){r_{{i,s}}^{u}(t)}}{{\kappa D_s}}, \forall i \in {\mathcal I}_s^u, s \in {\mathcal S}^u$$ where $r_{i,s}^u(t)$ denotes channel uses for transmitting a URLLC packet, $\kappa$ is a constant reflecting the number of channel uses per unit time per unit bandwidth of FDMA frame structure and numerology, $b_{i,s}^u(t)$ is an indicator variable reflecting whether the QoS requirement of device $i$ in slice $s$ can be satisfied at $t$. As network resources are limited and shared by all network slices, not all URLLC devices can be guaranteed to be served at every minislot. Certainly, we can adjust the slice priority weight that will be introduced in the following section to guide the resource orchestration for enforcing the entire URLLC devices coverage.

Based on the result in (\[eq:omega\_s\]) and the conclusion of the Lemma 3.2 in our previous work [@How2019yang], we can derive the minimum upper bound of bandwidth orchestrated for URLLC slices in the following lemma.

\[lem:upper\_bound\_urllc\_bandwidth\] At minislot $t$, for a given $M/M/W^u$ queue with packet arrival intensity $\bm \lambda$ and a family of packet transmit rates $\{{\kappa/r_{i,s}^u(t)}\}$, let $W^u({\bm r}(t))$ denote the minimum upper bound of bandwidth orchestrated for URLLC slices such that $P_Q^{M/M/W^u}\le \varsigma$ and $P_b({\bm \omega}^u(t), {\bm \lambda}, {\bm D}$, $W^u({\bm r}(t)))$ is of the order of $\alpha$, where $P_Q^{M/M/W^u}$ represents the queueing probability, and ${\bm D} = \{D_1, \ldots, D_{|{\mathcal S}^u|}\}$. If $\varsigma > \alpha$, then we have $$\label{eq:URLLC_bandwidth}
        \begin{array}{l}
            W^u(\bm r(t)) \approx \sum\limits_{s \in {\mathcal S}^u } {\sum\limits_{i \in {\mathcal I}_s^u} {{\lambda _{s}}b_{i,s}^u(t)\frac{{{{r_{i,s}^u(t)}}}}{\kappa}} } + \\
            \frac{{{\alpha} - \varsigma {\alpha}}}{{\varsigma  - {\alpha}}}\sqrt {\frac{\left({\sum\limits_{s \in {{\mathcal S}^u}} {\sum\limits_{i \in {\mathcal I}_s^u}{b_{i,s}^u(t)\lambda _s^2D_s^2}} }\right)\left( {\sum\limits_{s \in {\mathcal S}^u} {\sum\limits_{i \in {\mathcal I}_s^u} {{\lambda _{s}}b_{i,s}^u(t)\frac{{r_{i,s}^u(t)}^2}{{{\kappa ^2}{D_s}}}} }}\right)}{{\mathop {\min }\limits_{s \in {{\mathcal S}^u}} \{ {\lambda _s}{D_s}\} }}}
        \end{array}$$

We omit the proof here as the similar proof can be found in the proof section of Lemma 3.2 in [@How2019yang].

### Reduction of codeword error decoding probability

The crucial factor that impacts the codeword error decoding probability is the network capacity. Next, we discuss the relationship between the network capacity and codeword error decoding probability.

For any URLLC slice $s \in {\cal S}^u$, during minislot $t$, let $x_{i,s}^u(t)$ be the original data symbol destined to a URLLC device $i \in {\mathcal{I}_s^u}$ with ${\mathbb E}[|x_{i,s}^u(t)|^2] = 1$, $\bm g_{ij,s}(t) \in {\mathbb C}^K$ be the transmit beamformer pointing at the device $i$ from the $j$-th RRH and $\bm h_{ij,s}(t) \in \mathbb{C}^K$ be the channel coefficient between the $i$-th URLLC device and the $j$-th RRH. The channel coefficient may change over minislots. However, it is assumed to be i.i.d. over each minislot and remain unchanged during each minislot. Then, the received signal at device $i$ in $s$ during minislot $t$ is given by $$\label{eq:URLLC_receiving_signal}
{\hat x_{i,s}^u(t) = }{\sum\limits_{j \in {\cal J}} {{\bm h}_{ij,s}^{\rm{H}}} (t){{\bm g}_{ij,s}}(t)x_{i,s}^u(t) + {\sigma _{i,s}}(t)}, \forall i \in {\mathcal I}_s^u, s \in {\mathcal S}^u$$ where the first term is the useful signal for $i$ and $\sigma_{i,s}(t) \sim {\cal CN}(0, \sigma_{i,s}^2) $ is the additive white Gaussian noise (AWGN) experienced at $i$. Similar to [@tang2019service], interference signal is not involved in (\[eq:URLLC\_receiving\_signal\]) due to the utilization of a flexible FDMA mechanism. Then the SNR received at device $i$ in $s$ at minislot $t$ can be written as $$\label{eq:URLLC_SINR}
    SINR_{i,s}^u(t) = \frac{{|\sum\nolimits_{j \in {\cal J}} {{{\bm h}_{ij,s}^{\rm H}}{{(t)}}{{\bm g}_{ij,s}}(t)} {|^2}}}{{\phi \sigma _{i,s}^2}}, \forall i \in {\mathcal I}_s^u, s \in {\mathcal S}^u$$ where $\phi > 1$ is an SNR loss coefficient. The perception of perfect channel status information (CSI) or accurate channel coefficients requires the information exchange between an RRH and its associated device before data transmission, the process of which is generally time consuming. URLLC packets, however, have a stringent latency requirement. As a result, perfect CSI or accurate channel parameters may be unavailable for URLLC transmission, which may incur the SNR loss. The coefficient $\phi$ is then utilized to characterize the SNR loss [@hou2018burstiness].

Shannon capacity formula is created under a crucial assumption of transmitting a block with long enough blocklength. However, URLLC packets are typically very short to satisfy the ultra-low latency requirement. Thus, the famous Shannon capacity formula cannot be utilized to model the URLLC transmission data rate and capture the corresponding codeword error decoding probability. For URLLC transmission, we resort to the capacity analysis for a finite blocklength channel coding regime derived in [@yang2014quasi]. For any device $i \in {\mathcal{I}_s^u}$, $s \in \mathcal{S}^u$, the number of transmitted information bits $L_{i,s}^u(t)$ at minislot $t$ using $r_{i,s}^u(t)$ channel uses in AWGN channel can be accurately correlated with the codeword error decoding probability $\beta$ according to the following equation $$\label{eq:URLLC_bit_length}
    \begin{array}{l}
        L_{i,s}^u(t)  \approx  r_{i,s}^u(t)C(SNR_{i,s}^u(t)) - \\
\quad     {Q^{ - 1}}(\beta)\sqrt {r_{i,s}^u(t)V(SNR_{i,s}^u(t))}, \forall i \in {\mathcal I}_s^u, s \in {\mathcal S}^u
    \end{array}$$ where $C(SNR_{i,s}^u(t)) = \log_2(1 + SNR_{i,s}^u(t) )$ is the AWGN channel capacity under infinite blocklength assumption, $V(SNR_{i,s}^u(t)) = \ln^2 2\left( {1 - \frac{1}{{{{(1 + SNR_{i,s}^u(t))}^2}}}} \right)$ is the channel dispersion, $Q(\cdot)$ is the $Q$-function. It is noteworthy that a URLLC packet will usually be coded before transmission and the generated codeword will be transmitted in the air interface such that the transmission reliability can be improved.

The complicated expression of $V(SNR_{i,s}^u(t))$ in (\[eq:URLLC\_bit\_length\]) significantly hinders the theoretical analysis of network resources orchestrated for URLLC slices. Fortunately, as $V(SNR_{i,s}^u(t))$ is maximized by $\ln^2 2$, the closed-form expression of the minimum upper bound of $r_{i,s}^u(t)$ ($i \in {\mathcal I}_s^u, s \in {\mathcal S}^u $) with a codeword error decoding probability $\beta$ can be given by [@How2019yang] $$\label{eq:URLLC_channel_use}
    \begin{array}{l}
r_{i,s}^u(t) = \frac{{L_{i,s}^u(t)}}{{C(SINR_{i,s}^u(t))}} + \frac{{{{({Q^{ - 1}}(\beta))}^2}}}{{2{{(C(SINR_{i,s}^u(t)))}^2}}}\\
\quad + \frac{{{{({Q^{ - 1}}(\beta))}^2}}}{{2{{(C(SINR_{i,s}^u(t)))}^2}}}\sqrt {1 + \frac{{4L_{i,s}^u(t)C(SINR_{i,s}^u(t))}}{{{{({Q^{ - 1}}(\beta))}^2}}}}
\end{array}$$

Problem formulation
===================

This section aims to formulate the problem of RAN slicing for mIoT and bursty URLLC service multiplexing based on the above models.

In mIoT slices, each RRH may transmit feedback signal to its connected IoT devices for the connection establishment according to an RA four-step procedure [@grau2019preamble]. Meanwhile, in URLLC slices, each RRH may transmit URLLC packets to URLLC devices. As the transmit power $E_j$ of each RRH is limited, we have the following transmit power constraint $$\label{eq:RRH_energy}
\begin{array}{l}
   \sum\limits_{s \in {\cal S}} {{(1 + \alpha_g)\frac{\lambda_s^I}{\lambda_R}{{\hat E}_{j}^I}}}  + \\
   \qquad \sum\limits_{s \in {\cal S}^u} {\sum\limits_{i \in {\cal I}_s^u} {b_{i,s}^u(t){{\bm g}_{ij,s}^{\rm H}}{{(t)}}{{\bm g}_{ij,s}}(t)} }  \le {E_j}, \forall j \in \mathcal{J}
\end{array}$$ where ${\hat E}_{j}^I$ is assumed to be a constant and denotes the transmit power of the $j$-th RRH for connecting to its associated IoT devices over downlink, $\alpha_g$ represents a safety margin coefficient. As a PPP with intensity $\lambda_s^I$ is utilized to model the distribution of IoT devices, the actual number of IoT devices may be greater than $\lambda_s^I$ once deployed. As a result, the coefficient $\alpha_g$ is introduced to reserve transmit power for exceeded IoT devices.

In the RAN slicing system, as the total limited system bandwidth $W$ will be shared by mIoT slices and URLLC slices, we have the following bandwidth constraint $$\label{eq:total_bandwidth}
   \sum\limits_{s \in {\cal S}^I} {(1+\alpha_g)\omega_s(\bar t)}  + W^u(\bm r(t))  \le W$$ where $\omega_s(\bar t)$ denotes the bandwidth allocated to mIoT slice $s \in \mathcal{S}^I$ that is correlated with $F_s$ by means of $F_s = \left\lfloor {{{\omega _s}(\bar t)}/{{a}}} \right\rfloor$, and $\alpha_g \omega_s(\bar t)$ denotes a block of reserved bandwidth.

In (\[eq:total\_bandwidth\]), $F_s$ is an integer, and some integer variable recovery schemes [@tang2019systematic] can be leveraged to obtain the suboptimal $F_s$. However, considering the high computational complexity of optimizing an integer variable and the utilization of the scheme of spectrum safety margin, we directly relax the integer variable into a continuous one, i.e., let $F_s = {{{\omega _s}(\bar t)}/{{a}}}$. Without loss of any generality, we regard $\omega_s(\bar t)$ as an independent variable below. Besides, as at least one PRB should be allocated to mIoT slices, we have $$\label{eq:mMTC_bandwidth}
   \omega_s(\bar t) \ge a, \forall s \in {\cal S}^I$$

Owing to the exploration of mIoT and bursty URLLC service multiplexing, we should orchestrate network resources for all mIoT slices and URLLC slices to simultaneously maximize the utilities of mIoT slices and URLLC slices.

For a mIoT slice $s \in \mathcal{S}^I$, its primary goal is to offload as many data packets as possible from IoT devices. In this way, the number of accumulated packets in each IoT device should be kept at a low level. Considering that a great RA success probability of an IoT device will lead to a low number of accumulated packets in its queue, we define the utility of a mIoT slice as the following.

\[mydef:IoT\_slice\_utility\] Over a time slot of duration $T$, the mIoT slice utility is defined as the time-average of RA success probabilities of IoT devices in all mIoT slices, which is given by $$\label{eq:IoT_utility}
{{\bar U}^I} = \frac{1}{T}\sum\limits_{t = 1}^T {{U^I}(t)}  = \frac{1}{T}\sum\limits_{t = 1}^T {\tilde P(t) }$$ where $\tilde P(t) = \sum\limits_{s \in {{\mathcal S}^I}} {\frac{{{\lambda _s^I}{P_s}(t)}}{{\sum\nolimits_{s \in {{\mathcal S}^I}} {{\lambda _s^I}} }}}$ with the numerator $\lambda_s^I P_s (t)$ representing the expected sum of RA success probabilities of IoT devices in slice $s \in \mathcal{S}^I$ and the denominator $\sum\nolimits_{s\in\mathcal{S}^I}{\lambda_s^I}$ denoting a normalization coefficient.

In (\[eq:IoT\_utility\]), $\frac{\lambda_s^I}{\sum\nolimits_{s \in \mathcal{S}^I} {\lambda_s^I}}$ can be regarded as an intra-slice priority coefficient. A mIoT slice serving more IoT devices will be orchestrated with more network resources.

For a URLLC slice $s \in \mathcal{S}^u$, its primary objective is to maximize the slice gain that is reflected by the parameters in the slice request at a low cost. Therefore, we define an energy-efficient utility for URLLC slices, as presented below.

\[mydef:URLLC\_slice\_utility\] Over a time slot of duration $T$, the bursty URLLC slice utility is defined as the time-average energy efficiency for serving all URLLC devices, which is given by $$\label{URLLC_long_term_utility}
    \begin{array}{l}
        {{\bar U}^u} = \frac{1}{T}\sum\limits_{t = 1}^T {{U^u}(t)}  = \frac{1}{T}\sum\limits_{t = 1}^T {\sum\limits_{s \in {{\mathcal S}^u}} {U_s^u({D_s},{\bm g_{ij,s}}(t))} } \\
\qquad         = \frac{1}{T}\sum\limits_{t = 1}^T {\sum\limits_{s \in {{\mathcal S}^u}} {\sum\limits_{i \in {\mathcal I}_s^u} \frac{b_{i,s}^u(t)}{{1 - {e^{ - {D_s}}}}} } }  - \\
\qquad \quad        \frac{\eta }{T}\sum\limits_{t = 1}^T {\sum\limits_{s \in {{\mathcal S}^u}} {\sum\limits_{j \in {\mathcal J}} {\sum\limits_{i \in {\mathcal I}_s^u} {b_{i,s}^u(t)\bm g_{ij,s}^{\rm{H}}(t){\bm g_{ij,s}}(t)} } } }
    \end{array}$$ where $\eta$ is an energy efficiency coefficient reflecting a tradeoff between the URLLC slice gain and the RRH energy consumption.

Then, over a time slot of duration $T$, the RAN slicing problem for mIoT and URLLC service multiplexing can be formulated as follows

\[eq:original\_problem\] $$\begin{aligned}
{2}
& \mathop {\rm maximize}\limits_{{b_{i,s}^u(t)},{{\omega_s}}(\bar t),{{\bm g}_{ij,s}}(t)} {{\bar U}^I}  + {\tilde \rho} {{\bar U}^u} \\
& {\rm subject \text{ } to:} \nonumber \\
& b_{i,s}^u(t) \in \{ 0,1\}, \forall s \in {\cal S}^u, i \in {\cal I}_s^u \\
& \rm {constraints \text{ } (\ref{eq:mMTC_QoS}), (\ref{eq:RRH_energy})-(\ref{eq:mMTC_bandwidth}) \text{ } are \text{ } satisfied.}\end{aligned}$$

where $\tilde \rho$ is an inter-slice priority coefficient reflecting the priority of orchestrating network resources for mIoT slices and URLLC slices.

The mitigation of (\[eq:original\_problem\]) is quite challenging mainly because

-   **indeterministic objective function**: (\[eq:original\_problem\]) should be optimized at the beginning of the $1^{\rm st}$ minislot. The time-averaged objective function of (\[eq:original\_problem\]) can only be exactly computed according to the future channel information. Therefore, the value of the objective function is indeterministic at the beginning of the $1^{\rm st}$ minislot.

-   **two timescale issue**: the creation of a network slice is performed at a timescale of time slot. Thus, the variable $\omega_s(\bar t)$ should be determined at the beginning of the time slot $\bar t$ and kept unchanged over the whole time slot. The channel, however, is time-varying. As a result, the beamformer $\bm g_{ij,s}(t)$ should be optimized at each minislot. In summary, the variables in (\[eq:original\_problem\]) should be optimized at two different timescales.

-   **thorny optimization problem**: at each minislot $t$, the constraint (\[eq:mMTC\_QoS\]) is non-convex over $\omega_s(\bar t)$, and the constraints (\[eq:RRH\_energy\]), (\[eq:total\_bandwidth\]) are non-convex over $\bm g_{ij,s}(t)$, which together lead to a non-convex problem.

Problem solution with system generated channel
==============================================

This section aims to tackle these challenges by exploiting of an SAA technique [@kim2015guide], an ADMM method [@boyd2011distributed], a semidefinite relaxation scheme and a Taylor expansion scheme.

Sample average approximation and alternating direction method of multipliers
----------------------------------------------------------------------------

As mIoT slices and URLLC slices share the network resources, both $\bar U^I$ and $\bar U^u$ may be determined by channel coefficients experienced by URLLC slices. At each minislot $t$, due to the i.i.d. assumption on the channel coefficients of URLLC slices, we have $$\label{objfun_approx}
    \frac{1}{T}\sum\limits_{t = 1}^T {{U^I}(t)}  + \frac{1}{T}\sum\limits_{t = 1}^T {\tilde \rho {U^u}(t)}  \approx {E_{\hat {\bm h}}}\left[ {{{\hat U}^I} + \tilde \rho {{\hat U}^u}} \right]$$ where $\hat {\bm h}$ represents the channel samples of URLLC slices collected at the beginning of the time slot $\bar t$.

Given a collection of channel samples $\{\bm h_m\}$ with $\bm h_m = [\bm h_{11,1m};\ldots;\bm h_{1J,sm};\ldots;\bm h_{N^uJ,|\mathcal{S}^u|m}]$ and $m \in \mathcal{M}=\{1,\ldots,M\}$, Just like [@How2019yang], as constraints (\[eq:original\_problem\]b) and (\[eq:original\_problem\]c) construct a nonempty compact set, the conclusion of Proposition 5.1 in [@How2019yang] is applicable to this paper by exploiting the SAA technique. The conclusion indicates that if the number of channel samples $M$ is reasonably large, then $\frac{1}{M}\sum\limits_{m = 1}^M {U_m^I }  + \frac{\tilde \rho}{M}\sum\limits_{m = 1}^M {U_m^u }$ converges to ${E_{\bm {\hat h}}}\left[ {{{\hat U}^I} + {\tilde \rho {\hat U}^u}} \right]$ uniformly on the nonempty compact set almost surely. In other words, the SAA technique enables us to use the channel samples collected at the beginning of a time slot to approximate the unknown channel coefficients over the time slot. For notation lightening, we write $x_m$ instead of $x(m)$ that represents a variable corresponding to the channel sample $\bm h_m$.

Recall that the variable $\omega_s(\bar t)$ will be kept unchanged over the time slot $\bar t$ and the beamformer $\bm g_{ij,s}(t)$ should be calculated at each minislot $t$, we can further consider (\[eq:original\_problem\]) as a global consensus problem, which can be effectively mitigated by an ADMM method. In this problem, $\omega_s(\bar t)$ is a global consensus variable that should be maintained in consensus for all $\bm h_m$, and $\bm g_{ij,sm}$ that is calculated based on $\bm h_m$ is a local variable.

The fundamental principle of an ADMM method is to impose augmented penalty terms characterizing global consensus constraints on the objective function of an optimization problem. In this way, the local variables can be driven into the global consensus while still attempting to maximize the objective function. Let ${\bm G}_{i,sm} = {\bm g}_{i,sm} {\bm g}_{i,sm}^{\rm H} \in {\mathbb R}^{JK \times JK}$, ${\bm H}_{i,sm} = {\bm h}_{i,sm} {\bm h}_{i,sm}^{\rm H} \in {\mathbb R}^{JK \times JK}$, where ${\bm g}_{i,sm} = [{\bm g}_{i1,sm};\ldots;{\bm g}_{iJ,sm}] \in {\mathbb C}^{JK \times 1}$ and ${\bm h}_{i,sm} = [{\bm h}_{i_1,sm};\ldots;{\bm h}_{i_J,sm}] \in {\mathbb C}^{JK \times 1}$. By applying the matrix property ${{{\bm G}_{i,sm}} = {{\bm g}_{i,sm}}{\bm g}_{i,sm}^{\rm H} \Leftrightarrow {{\bm G}_{i,sm}} \succeq 0,} \text{ } {{\rm rank}({{\bm G}_{i,sm}}) \le 1}$ and utilizing the conclusions of SAA and ADMM, we can approximate (\[eq:original\_problem\]) as the following problem at the beginning of the time slot $\bar t$

\[eq:SAA\_admm\_problem\] $$\begin{aligned}
{2}
& \mathop {{\rm{minimize}}}\limits_{\{ \omega _{sm},\omega_s(\bar t),{b_{i,sm}^u},{\bm G_{i,sm}}\}\hfill}
\sum\limits_{m = 1}^M {\left[ { - \frac{{U_m^I}}{M} - \frac{{\tilde \rho U_m^u}}{M}} \right]}  + \nonumber \\
& \underbrace {\sum\limits_{m = 1}^M {\sum\limits_{s \in {{\mathcal S}^I}} {\left[ {{\psi _{sm}}\left( {{\omega _{sm}} - {\omega _s}(\bar t)} \right) + \frac{\mu }{2}\left\| {{\omega _{sm}} - \omega{_s}(\bar t)} \right\|_2^2} \right]} } }_{{\rm augmented} \text{ } {\rm penalty} \text{ } {\rm terms}} \\
& {\rm subject \text{ } to:} \nonumber \\
& P_{sm} \ge \pi_s, \forall s \in {\cal S}^I, m \in {\cal M} \\
& \sum\limits_{s \in {\cal S}^I} {{(1 + \alpha_g)\frac{\lambda_s^I}{\lambda_R}{\hat E}_{j}^I}}  + \sum\limits_{s \in {{\cal S}^u}} {\sum\limits_{i \in {\cal I}_s^u} {b_{i,sm}^u {\rm tr}({{\bm Z}_j}{{\bm G}_{i,sm}})} }  \le {E_j} \nonumber \\
& \qquad \qquad \forall j\in {\cal J},m\in {\cal M} \\
& \sum\limits_{s \in {{\cal S}^I}} {(1+\alpha_g) \omega _{s}(\bar t)}  + W^u(\bm r_m)  \le W, m \in {\cal M} \\
&{{\bm G}_{i,sm}} \succeq 0, \forall s \in {\cal S}^u, i \in {\cal I}_s^u, m \in {\cal M} \\
&{{\rm rank}({{\bm G}_{i,sm}}) \le 1}, \forall s \in {\cal S}^u, i \in {\cal I}_s^u, m \in {\cal M}  \\
& b_{i,sm}^u \in \{0,1\}, \forall s \in {\cal S}^u, i \in {\cal I}_s^u, m \in {\cal M}\end{aligned}$$

where $\psi_{sm}$ is the Lagrangian multiplier, $\mu$ is a penalty coefficient, ${\bm Z}_j$ is a square matrix with $J \times J$ blocks, and each block in ${\bm Z}_j$ is a $K \times K$ matrix. In ${\bm Z}_j$, the block in the $j$-th row and $j$-th column is a $K \times K$ identity matrix, and all other blocks are zero matrices.

(\[eq:original\_problem\]) is now reduced to a deterministic single timescale problem (\[eq:SAA\_admm\_problem\]). What is more, (\[eq:SAA\_admm\_problem\]) can be split into $M$ separate problems that can be optimized in parallel as its objective function is separable. Thus, the following ADMM-based framework from (\[eq:arg\_lagarangian\]) to (\[eq:psi\_update\]) can be exploited to mitigate (\[eq:SAA\_admm\_problem\])

\[eq:arg\_lagarangian\] $$\begin{aligned}
{2}
& \left\{ {\begin{array}{*{20}{l}}
{\omega _{sm}^{(k + 1)},b_{i,sm}^{u(k+1)}}\\
{\bm G_{i,sm}^{(k + 1)}}
\end{array}} \right\} = \mathop {{\rm{argmin}}}\limits_{\left\{ {\scriptstyle\omega _{sm},b_{i,sm}^{u},\hfill\atop
\scriptstyle{\bm G_{i,sm}}\hfill} \right\}} \overline {\cal L} (\omega _{sm},{{\bm G}_{i,sm}})  \\
& {\rm subject \text{ } to:} \nonumber \\
& {\rm for} \text{ } {\rm the} \text{ }  m-{\rm th} \text{ } {\rm sample},  (\ref{eq:SAA_admm_problem}b)-(\ref{eq:SAA_admm_problem}g) \text{ } \rm{are} \text{ } satisfied.\end{aligned}$$

$$\label{eq:omega_update}
\omega _s^{(k + 1)}(\bar t) = \frac{1}{M}\sum\limits_{m = 1}^M {\left( {\omega _{sm}^{(k + 1)} + \frac{1}{\mu }\psi _{sm}^{(k)}} \right)}, \text{ } \forall s\in {\cal S}^I$$

$$\label{eq:psi_update}
    \psi _{sm}^{(k + 1)} = \psi _{sm}^{(k)} + \mu \left( {\omega _{sm}^{(k + 1)} - \omega _s^{(k + 1)}(\bar t)} \right), \text{ } \forall s\in {\cal S}^I$$

where the augmented partial Lagrangian function $$\label{eq:average_Lagrangian_func}
\begin{array}{l}
    \bar{\mathcal{L}}(\omega _{sm},\bm {\bm G}_{i,sm}) =  {- \frac{{U_m^{I(k)}}}{M} - \frac{{\tilde \rho U_m^{u(k)}}}{M} +} \\
    { \sum\limits_{s \in {{\cal S}^I}} {\left[ {\psi _{sm}^{(k)}\left( {\omega _{sm} - \omega _s^{(k)}(\bar t)} \right) + \frac{\mu }{2}{{\left\| {\omega _{sm} - \omega _s^{(k)}(\bar t)} \right\|}_2^2}} \right]} }
    \end{array}$$

This ADMM-based framework can be executed on multiple processors. Each processor is responsible for optimizing (\[eq:arg\_lagarangian\]) and calculating (\[eq:psi\_update\]) with a global value as an input. (\[eq:omega\_update\]) is centrally updated in such a way that local variables converge to the global value, which is the solution of (\[eq:SAA\_admm\_problem\]). Unfortunately, (\[eq:arg\_lagarangian\]) is a mixed-integer non-convex optimization problem as there are zero-one variables, continuous variables and non-convex constraints in (\[eq:arg\_lagarangian\]). As a result, the optimization of (\[eq:arg\_lagarangian\]) is quite difficult. We next discuss how to handle this hard problem.

Alternative optimization
------------------------

In this subsection, we exploit a widely applied scheme, i.e., an alternative optimization scheme, to handle the mixed-integer non-convex optimization problem. Specifically, we first assume that continuous variables are known and attempt to mitigate a zero-one optimization problem. Given the zero-one variables, we then try to optimize a non-convex optimization problem. The process is alternatively conducted until convergence.

### URLLC device associations

Given continuous variables $\{{\bm G}_{i,sm}^{(k)}, \omega_{sm}^{(k)}\}$ at the $k$-th iteration, the association problem of URLLC devices in URLLC slices can take the following form

\[eq:arg\_lagarangian\_b\_isu\] $$\begin{aligned}
{2}
& \{ b_{i,sm}^{u(k+1)} \} = \mathop {{\rm{argmin}}}\limits_{\{b_{i,sm}^u\}} - \frac{{\tilde \rho U_m^{u(k)}}}{M}  \\
& {\rm subject \text{ } to:} \nonumber \\
& {\rm for} \text{ } m, \text{ } (\ref{eq:SAA_admm_problem}c), (\ref{eq:SAA_admm_problem}d),(\ref{eq:SAA_admm_problem}g) \text{ } \rm{are} \text{ } satisfied.\end{aligned}$$

This problem is non-linear and hard to be handled. In theory, an exhaustive algorithm can obtain the optimal solution of (\[eq:arg\_lagarangian\_b\_isu\]). The computation complexity of this algorithm is $O(2^{N^u})$ that may be impractical in implementation. Therefore, a greedy scheme of the computational complexity $O(N^u)$, which is summarized as the following, is proposed to obtain $\{b_{i,sm}^{u(k+1)}\}$

a)  initialize two device sets, i.e., candidate device set ${\mathcal{I}}^{u-} = {\mathcal I}^u$, association device set ${\mathcal I}^{u+} = \emptyset$.

b)  select the device that maximizes $\frac{{\tilde \rho U_m^{u(k)}}}{M}$ from ${\mathcal{I}}^{u-}$, remove it from ${\mathcal{I}}^{u-}$, and add it to ${\mathcal{I}}^{u+}$. Given ${\mathcal{I}}^{u+}$, check the feasibility of (\[eq:arg\_lagarangian\_b\_isu\]). If (\[eq:arg\_lagarangian\_b\_isu\]) is feasible, then accept the device; otherwise, remove the device from ${\mathcal{I}}^{u+}$. Continue till ${\mathcal{I}}^{u-} = \emptyset$.

### joint bandwidth and beamforming optimization

Given the obtained $b_{i,sm}^{u(k+1)}$, (\[eq:SAA\_admm\_problem\]) will be reduced to the following joint bandwidth and beamforming problem.

\[eq:bandwidth\_beamforming\_problem\] $$\begin{aligned}
{2}
& \left\{
{\omega _{sm}^{(k + 1)}},
{\bm G_{i,sm}^{(k + 1)}}
 \right\} = \mathop {{\rm{argmin}}}\limits_{\{\omega _{sm},{\bm G_{i,sm}}\}} \overline {\cal L} (\omega _{sm},{{\bm G}_{i,sm}})  \\
& {\rm subject \text{ } to:} \nonumber \\
& {\rm for} \text{ } m,  (\ref{eq:SAA_admm_problem}b)-(\ref{eq:SAA_admm_problem}f) \text{ } \rm{are} \text{ } satisfied.\end{aligned}$$

In (\[eq:bandwidth\_beamforming\_problem\]), the low-rank constraint (\[eq:SAA\_admm\_problem\]f) is non-convex, and its objective function is not convex and even not quasi-convex w.r.t. ${\omega_{sm}}$, the tackling of which is quite tricky.

To tackle the non-convex low-rank constraint (\[eq:SAA\_admm\_problem\]e), we resort to the semidefinite relaxation technique. The primary procedures of SDR are i) directly drop the low-rank constraint; ii) solve the optimization problem without the low-rank constraint to obtain the solution; iii) owing to the relaxation, the obtained solution cannot satisfy the low-rank constraint in general. If it is, its principal component is the optimal solution to the problem. If not, then some manipulations such as randomization/scale [@ma2010semidefinite] are needed to perform on the solution to impose the low-rank constraint.

For the tricky objective function, we are reminded of the art of dealing with a non-convex function, i.e., study the structure of the function if it is non-convex. A crucial observation is that $P_{sm}$ is quasi-concave w.r.t. $\omega_{sm}$ although the objective function is not quasi-convex w.r.t. $\omega_{sm}$. Therefore, we resort to the Taylor expansion to approximate the tricky objection function.

The following analysis is built under the following two facts

-   the value of the objective function of (\[eq:bandwidth\_beamforming\_problem\]) is mainly determined by that of ${\tilde P_{m}^{(k)}}$ (or $U_m^{I(k)}$);

-   for all $s \in \mathcal{S}^u$, $m \in {\mathcal{M}}$, the solution $\omega_{sm}$ maximizing ${\tilde P_{m}^{(k)}}$ must locate in the range of $[\omega_{sm}^{lb}, S_{sm}^{\star}]$ shown in Fig. \[fig:fig\_quasi\_convex\_structure\], where $\omega_{sm}^{lb}$ denotes the lower bound of $\omega_{sm}$ satisfying the constraint (\[eq:SAA\_admm\_problem\]b), $S_{sm}^{\star}$ is the $\omega_{sm}$ maximizing $P_{sm}$, and the notation $P_{sm}|_{\omega_{sm}}$ is utilized to explicitly indicate that $P_{sm}$ is a function of $\omega_{sm}$.

Fact 1 holds because the linear terms w.r.t. $\omega_{sm}$ will donate little to the objective function as the consensus constraint is active. Besides, the quadratic terms pull local values towards the consensus; thus, they will also donate little to the objective function. Fact 2 holds because the total bandwidth is limited and shared. For example, given a value $\omega_{sm, 2} \in [S_{sm}^{\star}+\delta_{\omega}, W]$ with $\delta_{\omega}$ being a small positive constant, there must exist a value $\omega_{sm,1} \in [\omega_{sm}^{lb}, S_{sm}^{\star}]$ such that $P_{sm}|_{\omega_{sm,1}} = P_{sm}|_{\omega_{sm,2}}$. Thus, a small $\omega_{sm}$ will be preferred as it indicates that more bandwidth can be allocated to URLLC slices to further improve the objective function.

For all $s \in \mathcal{S}^I$, it can be proved that $P_{sm}$ is concave in the interval $(a_1, a_2]$ by evaluating the second-order derivative of $P_{sm}$. Therefore, we can leverage the $2^{\rm nd}$ degree Taylor expansion to approximate $P_{sm}$ in this interval. Considering that $P_{sm}$ is convex in the interval $[\omega_{sm}^{lb}, a_1]$, the $1^{\rm st}$ degree Taylor expansion is always leveraged to obtain the lower bound of $P_{sm}$. However, this interval is usually rather narrow, and the value of $P_{sm}$ in this interval is much lower than the value of that in the interval $(a_1, a_2]$. What is more, the error bound of the $1^{\rm st}$ degree Taylor expansion is greater than that of the $2^{\rm nd}$ expansion. Therefore, we explore the $2^{\rm nd}$ degree Taylor expansion to approximate $P_{sm}$ in the interval $[\omega_{sm}^{lb}, a_2]$. Fig. \[fig\_Taylor\_expansion\] shows an example of the $2^{\rm nd}$ degree Taylor expansion of $P_{sm}$. Given a local point $\bm \omega_{m}^{(k,q)}$ at the $q$-th iteration, the Taylor expansion of $-\tilde P_{m}^{(k)}$ at the local point can be given by $$\label{eq:P_sm_Taylor_expansion}
    \begin{array}{*{20}{l}}
{ - {{\tilde P}_{m}^{(k)}} \approx  - \tilde P_{m}^{(k,q)} - \nabla \tilde P_m^{(k,q)}{{({{\bm \omega} _m} - {\bm \omega} _m^{(k,q)})}^{\rm T}} - }\\
{\frac{1}{2}({{\bm \omega} _m} - {\bm \omega} _m^{(k,q)})H({\bm \omega} _{m}^{(k,q)}){{({\bm \omega _m} - {\bm \omega} _m^{(k,q)})}^{\rm T}} - R_2(\bm \omega_{m})}
\end{array}$$ where ${\bm \omega}_m = [\omega_{1m}, \ldots, \omega_{|{\mathcal S}^I|m}]$, $\nabla \tilde P_m^{(k,q)}$ is the gradient of $\tilde P_m^{(k)}$ over ${\bm \omega_m}$ at the local point ${\bm \omega_m^{(k,q)}}$ with $$\label{eq:first_order_Taylor_expansion}
    \begin{array}{l}
\frac{{\partial P_{sm}^{(k)}}}{{\partial \omega _{sm}^{(k,q)}}} = \frac{{{\lambda _s^I}(1 + {\varpi _s}{\rho _o}){e^{ - {\varpi _s}{\sigma ^2}}}}}{{\sum\nolimits_{s \in {{\cal S}^I}} {{\lambda _s^I}} }} \times \\
\qquad \qquad \left[ {\frac{{3.5{y_{sm}}{z_s}\omega _{sm}^{2.5(k,q)}}}{{{{({y_{sm}}{z_s} + \omega _{sm}^{(k,q)})}^{4.5}}}} - \frac{{3.5{y_{sm}}\omega _{sm}^{2.5(k,q)}}}{{{{({y_{sm}} + \omega _{sm}^{(k,q)})}^{4.5}}}}} \right]
\end{array}$$ and $H({\bm \omega}_m^{(k,q)})$ is a Hessian matrix with $$\label{eq:second_order_Taylor_expansion}
\begin{array}{l}
\frac{{{\partial ^2}P_{sm}^{(k)}}}{{\partial \omega _{sm}^{2(k,q)}}} = \frac{{\lambda _s^I(1 + {\varpi _s}{\rho _o}){e^{ - {\varpi _s}{\sigma ^2}}}}}{{\sum\nolimits_{s \in {{\mathcal S}^I}} {\lambda _s^I} }}\left[ {\frac{{15.75y_{sm}^2z_s^2\omega _{sm}^{1.5(k,q)}}}{{{{({y_{sm}}{z_s} + \omega _{sm}^{(k,q)})}^{5.5}}}} - } \right.\\
\frac{{7{y_{sm}}{z_s}\omega _{sm}^{1.5(k,q)}}}{{{{({y_{sm}}{z_s} + \omega _{sm}^{(k,q)})}^{4.5}}}} + \frac{{7{y_{sm}}\omega _{sm}^{1.5(k,q)}}}{{{{({y_{sm}} + \omega _{sm}^{(k,q)})}^{4.5}}}} - \left. {\frac{{15.75y_{sm}^2\omega _{sm}^{1.5(k,q)}}}{{{{({y_{sm}} + \omega _{sm}^{(k,q)})}^{5.5}}}}} \right]
\end{array}$$ $$\label{eq:second_order_intersection_term}
  \frac{{{\partial ^2}P_{sm}^{(k)}}}{{\partial \omega _{sm}^{(k,q)}\partial \omega _{s'm}^{(k,q)}}} = 0, \forall s \ne s'$$ $y_{sm} = \frac{{a{P_{nr,sm}}{P_{ne,sm}}{\lambda_s^I}}}{{3.5{\lambda _R}}}$, ${z_s} = \frac{{\theta_s^{th}}}{{ 1 + \theta_s^{th} }}$. Besides, we write $\omega_{sm}^{2.5(k,q)}$ rather than ${\left( {\omega _{sm}^{(k,q)}} \right)^{2.5}}$ for lightening the notation.

\[lemma:error\_bound\] Let the function $\tilde P_{m}^{(k)}: \mathbb{R}^{|{\mathcal S}^I|} \to \mathbb{R}$ be three times differentiable in a given interval $[\omega_{sm}^{lb}, S_{sm}^{\star}]$ for all $s \in {\mathcal{S}^I}$, then the error bound of the $2^{\rm nd}$ degree Taylor expansion of $\tilde P_{m}^{(k)}$ at the local point $\bm \omega_{sm}^{(k,q)}$ with $\omega_{sm}^{(k,q)} \in [\omega_{sm}^{lb}, S_{sm}^{\star}]$ is given by $$\label{eq:error_bound_lemma}
    {R_2}({\bm \omega _m}) = \frac{1}{{3!}}{\left[ {\sum\limits_{s \in {{\cal S}^I}} {\left( {{\omega _{sm}} - \omega _{sm}^{(k,q)}} \right)\frac{\partial }{{\partial \omega _{sm}^{(k,q)}}}} } \right]^3} \tilde P_m^{(k)}{{|_{\bm \omega _m^{lb}}^{{\bm S_m^ \star }}}}$$ where $\tilde P_m^{(k)}{{|_{\bm \omega _m^{lb}}^{{\bm S_m^ \star }}}} = \max \left\{ {\tilde P_m^{(k)}{|_{\bm \omega _m^{lb}}},\tilde P_m^{(k)}{|_{\bm S_m^ \star }}} \right\}$, $\bm \omega_{m}^{lb} = [\omega_{1m}^{lb}, \ldots, \omega_{|\mathcal{S}^I|m}^{lb}]$ and $\bm S_{m}^{\star} = [S_{1m}^{\star}, \ldots, S_{|\mathcal{S}^I|m}^{\star}]$.

Please refer to Appendix C.

After conducting the $2^{\rm nd}$ degree Taylor approximation, the objective function becomes a convex function. Although the constraint (\[eq:SAA\_admm\_problem\]b) is $P_{sm}$ related, we need not to conduct the Taylor approximation on (\[eq:SAA\_admm\_problem\]b) as $P_{sm}$ is quasi-concave and unimodal. In fact, the probability constraint (\[eq:SAA\_admm\_problem\]b) is equivalent to the following inequality $$\label{eq:variable_range}
    \omega_{sm}^{lb} \le \omega_{sm} \le \omega_{sm}^{ub}, \forall s \in \mathcal{S}^I$$ where $\omega_{sm}^{ub} \le W$ represents the upper bound of $\omega_{sm}$ satisfying (\[eq:SAA\_admm\_problem\]b).

Next, a low-complexity bisection-search-based scheme, the main procedures of which are described below, is developed to obtain $\omega_{sm}^{lb}$, $S_{sm}^{\star}$, and $\omega_{sm}^{ub}$

a)  let the function $Q_{sm} = P_{sm} - \pi_s$. Perform the bisection search method [@yang2018three] on $Q_{sm} = 0$ to obtain $\omega_{sm}^{lb}$ and $\omega_{sm}^{ub}$ that are the two zero points of $Q_{sm}$.

b)  with the obtained $\omega_{sm}^{lb}$ and $\omega_{sm}^{ub}$, find the maximum value $S_{sm}^{\star}$ of $P_{sm}$ using the bisection search method again.

According to the above analysis, at the $q$-th iteration, we can rewrite (\[eq:bandwidth\_beamforming\_problem\]) as

\[eq:SCA\_bandwidth\_beamforming\_problem\] $$\begin{aligned}
{2}
& \left\{
{\omega _{sm}^{(k + 1,q+1)}},{\bm G_{i,sm}^{(k + 1,q+1)}}
 \right\} = \mathop {{\rm{argmin}}}\limits_{\{\omega _{sm},{\bm G_{i,sm}}\}} \bar {\cal L}^{(q)} (\omega _{sm},{\bm G_{i,sm}})  \\
& {\rm subject \text{ } to:} \nonumber \\
& {\rm for} \text{ } m,   (\ref{eq:SAA_admm_problem}c)-(\ref{eq:SAA_admm_problem}e),(\ref{eq:variable_range}) \text{ } \rm{are} \text{ } satisfied.\end{aligned}$$

where $$\label{eq:average_Lagrangian_func}
\begin{array}{l}
    \bar {\mathcal{L}}^{(q)}(\omega _{sm},\bm G_{i,sm}) =  -{ \frac{1}{M}{\tilde P_{m}^{(k)}} - \frac{{\tilde \rho U_m^{u(k)}}}{M} +} \\
    { \sum\limits_{s \in {{\cal S}^I}} {\left[ {\psi _{sm}^{(k,q)}\left( {\omega _{sm} - \omega _s^{(k,q)}(\bar t)} \right) + \frac{\mu }{2}{{\left\| {\omega _{sm} - \omega _s^{(k,q)}(\bar t)} \right\|}_2^2}} \right]} } \nonumber
    \end{array}$$

In (\[eq:SCA\_bandwidth\_beamforming\_problem\]), the objective function is convex, (\[eq:SAA\_admm\_problem\]c) is affine, and the constraint (\[eq:SAA\_admm\_problem\]d) can be proved to be convex w.r.t. both $\omega_{sm}$ and $\bm G_{i,sm}$ [@How2019yang]. Therefore, (\[eq:SCA\_bandwidth\_beamforming\_problem\]) is a convex problem that can be effectively mitigated by some standard convex optimization tools such as CVX [@CVX] and MOSEK [@MOSEK].

Then we can summarize the main steps of mitigating the problem (\[eq:SAA\_admm\_problem\]) in Algorithm \[alg1\].

**Initialization:** Randomly initialize $\bm G_{i,s}^{(0,0)}$, $\{\omega_{s}^{(0,0)}\}$, let $k_{\rm max}=250$, $q_{\rm max}=250$, $q = 0$, $k = 0$, and generate channel samples {$\bm H_{i,sm}$}. Given $\bm G_{i,sm}^{(k,q)}$, $\omega_{sm}^{(k,q)}$, call the greedy scheme to obtain $b_{i,sm}^{u(k,q+1)}$. Optimize (\[eq:SCA\_bandwidth\_beamforming\_problem\]) with obtained $b_{i,sm}^{u(k,q+1)}$ to achieve $\bm G_{i,sm}^{(k,q+1)}$ and $\omega_{sm}^{(k,q+1)}$. Update $q = q + 1$. Let $\omega_{sm}^{(k+1,q+1)} = \omega_{sm}^{(k,q+1)}$, update $\psi _{sm}^{(k + 1)}$ using (\[eq:psi\_update\]). Call (\[eq:omega\_update\]) to update $\omega_{s}^{(k+1)}(\bar t)$. Update $k = k + 1$.

\[lemma\_5\] For all $i\in\mathcal{I}_s^u$, $s \in {\mathcal{S}^u}$, and $m \in {\mathcal{M}}$, the obtained power matrix $\bm G_{i,sm}^{(k,q)}$ by Algorithm 1 at the $(k,q)$-th iteration satisfies the low-rank constraint, i.e., the SDR for the power matrix utilized in Algorithm 1 is tight.

Please refer to Appendix D.

Optimization of beamforming with system sensed channels
=======================================================

In Section V, we obtained a family of global consensus variables $\{\omega_{s}(\bar t)\}$ with the system generated channel samples. The time-varying actual channels may require the re-optimization of beamformers and device associations at each minislot. According to system sensed channels at each minislot, we next discuss how to calculate beamforms and device associations.

At each minislot $t$, given the global consensus variables $\{\omega_s(\bar t)\}$, the original problem (\[eq:original\_problem\]) will be reduced to the following problem

\[eq:mini\_time\_scale\_transformed\_problem\] $$\begin{aligned}
{2}
& \mathop {{\rm{maximize}}}\limits_{\{b_{i,s}^u(t), {\bm G}_{i,s}(t)\}} \text{ } {\tilde \rho} {{U}^u(t)}  \\
& {\rm subject \text{ } to:} \nonumber \\
& \rm {constraints \text{ } (\ref{eq:RRH_energy}), (\ref{eq:total_bandwidth}),(\ref{eq:original_problem}b) \text{ } are \text{ } satisfied.}\end{aligned}$$

In (\[eq:mini\_time\_scale\_transformed\_problem\]), the channels are system sensed ones at $t$. According to the convexity analysis in Section V, (\[eq:mini\_time\_scale\_transformed\_problem\]) is a mixed-integer non-convex programming problem with positive semidefinite matrices, which is hard to be mitigated. Therefore, the alternative optimization scheme presented in subsection V-B can be leveraged to achieve the solutions $b_{i,s}^u(t)$ and $\bm G_{i,s}(t)$ of (\[eq:mini\_time\_scale\_transformed\_problem\]). Lemma \[lemma\_5\] indicates that the achieved ${\rm rank}(\bm G_{i,s}(t)) \le 1$. Thus, we can obtain the beamformers $\bm g_{i,s}(t)$ by performing the eigendecomposition on $\bm G_{i,s}(t)$. To sum up, over a time slot $\bar t$, the slice resource optimization algorithm designed for the RAN slicing system can be summarized as follows.

**Initialization:** $\{{\bm H}_{i,s}(t)\}$, $\forall i \in {{\mathcal I}^u}$, $s \in {{\mathcal S}^u}$, and let $P_{s}^{1} \in [0,1]$, $\mu_{a,s}^{1}=0$, $\forall s \in {{\mathcal S}^I}$. Call Algorithm \[alg1\] to obtain $\{\omega_s(\bar t)\}$, for all $s\in {\mathcal S}^I$. Given $\{\omega_s(\bar t)\}$, mitigate (\[eq:mini\_time\_scale\_transformed\_problem\]) by exploiting the alternative optimization scheme to obtain beamformers $\{{\bm g}_{i,s}(t)\}$ and URLLC device associations $b_{i,s}^u(t)$ for all $i \in {\mathcal I}_s^u$, $s \in {\mathcal S}^u$.

Simulation results
==================

In this section, we aim at evaluating the proposed algorithm via extensive simulations.

Comparison algorithms and parameter setting
-------------------------------------------

We compare the following three algorithms to verify the effectiveness of the proposed algorithm and to explain the impact of access control schemes on the RAN system performance intuitively i) SRO algorithm that adopts the unrestricted access control scheme; ii) SRO-ACB$_{\rm I}$ algorithm that utilizes the ACB access control scheme with $P_{ACB} = 0.9$; iii) SRO-ACB$_{\rm II}$ algorithm that adopts the ACB access control scheme with $P_{ACB} = 0.5$.

The parameter setting is as follows: RRHs and IoT devices are deployed following independent PPPs in a one km$^2$ area. URLLC devices are randomly and uniformly distributed in this area. There are three mIoT slices and two URLLC slices in the RAN slicing system. For the mIoT slices, set the new endogenous packet arrivals rate $\epsilon_{w,s}(t) = [1.5, 1.0, 0.5]$ packets/minislot, $\pi_s = 0.5$, $\forall s,t$. Let the path-loss component $\varphi = 4$, $L = 2000$ bits, $\tau = 1$ unit, $\sigma^2 = -90$ dBm, $\rho_o = -90$ dBm, ${\hat E}_{j}^I = 0.03$ mW, $\lambda_R = 3$ RRHs/km$^2$, $\{\lambda_{s}^I\} = [18000, 18000, 18000]$ IoT devices/km$^2$, $\{\gamma_s^{th}\} = \{5.8, 4.35, 2.9\}$ Kbits/minislot. For the URLLC slices, the transmit antenna gain at each RRH is set to be $5$ dB, and a log-normal shadowing path-loss model is leveraged to simulate the path-loss between an RRH and a URLLC device with the log-normal shadowing parameter being $10$ dB. A path-loss is computed by $h({\rm dB}) = 128.1 + 37.6\log_{10}d$, where $d$ (in km) represents the distance between a device and an RRH. Let $L_{i,s}^u = 160$ bits, $\sigma_{i,s}^2 = -100$ dBm, $\lambda_{s} = \lambda = 0.1$ packets/minislot, $\forall i,s$, $\{I_s^u\} = \{3, 5\}$ devices, and $\{D_s\} = \{1, 2\}$ milliseconds, $E_j = 3$ W, $\forall j$ [@tang2019service]. Other simulation parameters are shown in Table \[table\_simulation\_parameters\].

    Para.       Value       Para.            Value              Para.               Value
  ---------- ----------- ------------ -------------------- --------------- -----------------------
     $J$         $3$         $K$              $2$           $\tilde \rho$            $1$
    $\eta$      $100$        $T$              $60$               $W$              $60$ MHz
     $M$        $100$       $\phi$           $1.5$               $a$             $0.18$ MHz
    $\xi$       $54$      $\alpha_g$         $0.05$           $\kappa$      $5.12 \times 10^{-4}$
   $\alpha$   $10^{-5}$    $\beta$     $2 \times 10^{-8}$    $\varsigma$     $2 \times 10^{-5}$

  : Simulation parameters[]{data-label="table_simulation_parameters"}

Performance evaluation
----------------------

To evaluate the comparison algorithms, the following performance indicators are utilized i) RA success probability $P_s(t)$ that is computed using (\[eq:QoS\_prob\_analysis\]); ii) expected queue length per IoT device at minislot $t$, $E[Q_s(t)] = \mu_{a,s}(t)$; iii) total slice utility $\bar U$ that is the objective function of (\[eq:original\_problem\]).

We first evaluate the convergence of the proposed SRO algorithm. Fig. \[fig:fig\_convergence\] illustrates the convergence of SRO with $\Delta_{\omega} = {\rm{ }}\sum\nolimits_{s \in {{\mathcal S}^I}} {\left| {\omega _s^{(k + 1)}(\bar t) - \omega _s^{(k)}(\bar t)} \right|} $. It shows that SRO can converge after several iterations.

![The convergence curve of the proposed SRO algorithm.[]{data-label="fig:fig_convergence"}](Convergence_curve.eps){width="2.9in"}

We next plot the tendency of the RA success probability $P_s(t)$ and the corresponding expected queue length $E[Q_s(t)]$ during a time slot in Fig. \[fig:fig\_queueLen\_Pst\]. Fig. \[fig:fig\_queueLen\_Pst\](a) and \[fig:fig\_queueLen\_Pst\](c) show the tendency of $P_s(t)$ and $E[Q_s(t)]$ in the case of $\{\gamma_s^{th}\} = \{1.8, 1.35, 0.9\}$. Fig. \[fig:fig\_queueLen\_Pst\](b) and \[fig:fig\_queueLen\_Pst\](d) depict the tendency of $P_s(t)$ and $E[Q_s(t)]$ in the case $\{\gamma_s^{th}\} = \{5.8, 4.35, 2.9\}$.

From Fig. \[fig:fig\_queueLen\_Pst\], we obtain the following interesting conclusions: the queue of each IoT device is not stable when the queue serving rate $\gamma_s^{th}$ is small. In this case, the average queue length monotonously increases over minislot $t$. On the contrary, the queue of each IoT device is periodically flushed when a great queue serving rate is configured.

![Trends of $P_s(t)$ and $E[Q_s(t)]$. (a) and (c) are results of the parameter setting $\{\gamma_s^{th}\} = \{1.8, 1.35, 0.9\}$ Kbits/minislot; (b) and (d) correspond to the parameter setting $\{\gamma_s^{th}\} = \{5.8, 4.35, 2.9\}$ Kbits/minislot.[]{data-label="fig:fig_queueLen_Pst"}](Ps_EQs_vs_t.eps){width="3.1in"}

Let the IoT device intensity $\bm \lambda^I = [900 n, 900 n, 900 n]$ with $n \in \{6, 8, \ldots, 26\}$. Under the existence of both mIoT and URLLC slices, we plot trends of the total slice utility $\bar U$ and bursty URLLC slice utility $\bar U^u$ w.r.t. $n$ in Fig. \[fig:fig\_utility\_vs\_lamdba\_IoT\] to understand the impact of the mIoT slices on the performance of all comparison algorithms. In this figure, $B = [b_{11}^u, \ldots, b_{31}^u, b_{12}^u, \ldots, b_{52}^u]$, $\omega^I=[\omega_{SRO}^I, \omega_{ACB_{\rm I}}^I, \omega_{ACB_{\rm II}}^I]$ MHz with $\omega_{SRO}^I$, $\omega_{ACB_{\rm I}}^I$ and $\omega_{ACB_{\rm II}}^I$ representing the bandwidth allocated to mIoT slices by executing SRO, SRO-ACB$_{\rm I}$, and SRO-ACB$_{\rm II}$ algorithms, respectively, and $\bar U^I = [\bar U_{SRO}^I, \bar U_{ACB_{\rm I}}^I, \bar U_{ACB_{\rm II}}^I]$ with $\bar U_{SRO}^I$ denoting the achieved mIoT slice utility of SRO.

\

The following observations can be obtained from Fig. \[fig:fig\_utility\_vs\_lamdba\_IoT\]: i) when $n < 16$, all algorithms almost obtain the same total slice utility, and the obtained utilities are robust to the average number of IoT devices; ii) when $16 \le n \le 26$, the conclusion changes. For the SRO algorithm, its achieved $\bar U$ decreases with an increasing $n$ due to increasing interference. A great $n$, however, does not cause a significant decrease in the total slice utilities obtained by SRO-ACB$_{\rm I}$ and SRO-ACB$_{\rm II}$. Thanks to the exploration of an access control scheme, both SRO-ACB$_{\rm I}$ and SRO-ACB$_{\rm II}$ can achieve greater $\bar U$ than SRO. For example, compared with SRO, SRO-ACB$_{\rm II}$ improves $\bar U$ by $6.65\%$ when $n = 24$; iii) when $n = 26$, which means that the total average number of IoT devices reaches $70,200$ devices, the RAN slicing system fails to create and manage mIoT slices as the QoS requirements of mIoT slices serving such a massive average number of devices cannot be simultaneously satisfied. In this case, all system resources are allocated to URLLC slices, and the maximum bursty URLLC slice utility is obtained; iv) as mIoT slices and URLLC slices share the system resources, an increasing $n$ results in a decreasing bursty URLLC slice utility $\bar U^u$; Besides, it is interesting to find that the two access-control-based algorithms may not outperform SRO in terms of obtaining $\bar U^u$. It indicates that URLLC slices do not benefit from access control schemes of mIoT slices when changing $n$; v) the RAN slicing system can always accommodate the QoS requirements of all URLLC devices.

Next, to understand the impact of URLLC slices on the performance of all comparison algorithms, we plot the trends of the total slice utilities and the mIoT slice utilities obtained by all comparison algorithms w.r.t. URLLC packet arrival rate $\lambda$ with $\lambda = \{0.1,0.5,1.0,\ldots,4.5,5.0\}$ packets per unit time in Fig. \[fig:fig\_utility\_vs\_lamdba\_URLLC\]. Similarly, the following notations are involved in this figure: $\omega^u = [\omega_{SRO}^u, \omega_{ACB_{\rm I}}^u, \omega_{ACB_{\rm II}}^u]$, $\bar U^u = [\bar U_{SRO}^u, \bar U_{ACB_{\rm I}}^u, \bar U_{ACB_{\rm II}}^u]$ with $\omega_{SRO}^u$ and $\bar U_{SRO}^u$ denoting the bandwidth allocated to URLLC slices and the URLLC slice utility obtained by running the SRO algorithm, respectively.

\

From Fig. \[fig:fig\_utility\_vs\_lamdba\_URLLC\], we can observe that: i) the obtained $\bar U$ of all algorithms decrease with $\lambda$ mainly due to the decrease of the bursty URLLC slice utility. Two algorithms adopting the access control scheme always achieve greater utilities $\bar U$ than SRO. For example, when $\lambda = 5$, compared with the SRO algorithm, the obtained $\bar U$ of SRO-ACB$_{\rm II}$ is increased by $29.41\%$; ii) for all algorithms, the computed bandwidth for URLLC slices increases with an increasing $\lambda$. However, their obtained URLLC slice utilities $\bar U^u$ are reduced owing to the increase of energy consumption; iii) SRO-ACB$_{\rm II}$ may achieve greater $\bar U$ than SRO-ACB$_{\rm I}$ as a greater $\bar U^I$ is obtained by reducing the number of interfering IoT devices; iv) the obtained mIoT slice utilities $\bar U^I$ of SRO-ACB$_{\rm I}$ and SRO-ACB$_{\rm II}$ are robust to the URLLC packet arrival rate. The obtained $\bar U^I$ of SRO decreases with an increasing $\lambda$; v) an important observation is that the $\bar U^I$ of the access-control-based SRO-ACB$_{\rm I}$ algorithm is $1.65$ times that of the SRO algorithm when $\lambda = 5$. It explicitly reflects that mIoT slices can still benefit from access control schemes even though $\lambda$ is changed.

Figs. \[fig:fig\_utility\_vs\_lamdba\_IoT\] and \[fig:fig\_utility\_vs\_lamdba\_URLLC\] illustrate the situation of a given total system bandwidth. We next change the total bandwidth $W$ and plot its impact on the obtained total slice utilities of all algorithms in Fig. \[fig:fig\_total\_utility\_vs\_bandwidth\].

![Trend of achieved total slice utility vs. system bandwidth.[]{data-label="fig:fig_total_utility_vs_bandwidth"}](total_slice_utility_vs_bandwidth.eps){width="3.1in"}

The following conclusions can be obtained from this figure i) when $W=45$ MHz, the QoS requirements of all IoT devices cannot be simultaneously satisfied. As a result, the total bandwidth is allocated to URLLC slices; ii) when $W$ locates in the range of $(45, 55]$ MHz, the achieved total slice utilities $\bar U$ of SRO and SRO-ACB$_{\rm I}$ increase with $W$. Owing to the utilization of the access control scheme, SRO-ACB$_{\rm I}$ and SRO-ACB$_{\rm II}$ obtain higher $\bar U$ than SRO. For example, compared with the SRO algorithm, the SRO-ACB$_{\rm II}$ algorithm improves the achieved total slice utility by $6.66\%$ when $W = 50$ MHz. iii) when $W > 55$ MHz, all algorithms cannot remarkably improve $\bar U$.

At last, we discuss other crucial parameters’ impact on the performance of the comparison algorithms. We reconfigure $\{\gamma_s^{th}\}$ of mIoT slices as $\gamma_1^{th} = 3.6 m$, $\gamma_2^{th} = 2.7 m$ and $\gamma_3^{th} = 1.8 m$ Kbits/minislot with $m \in \{1.5, 1.6, \ldots, 2.1\}$ and $\{D_s\}$ of URLLC slices as $D_1 = 0.00025 d$ second and $D_2 = 0.0005 d$ second with $d \in \{2, 3, \ldots, 10\}$. The impact of QoS requirements of network slices on the total slice utility is plotted in Figs. \[fig:fig\_total\_utility\_vs\_gammath\] and \[fig:fig\_total\_utility\_vs\_Ds\]. The impact of energy efficiency coefficient $\eta$ is plotted in Fig. \[fig:fig\_total\_utility\_vs\_eta\]. In this figure, we denote the energy consumption of RRHs of all algorithms by $E^u = [E_{SRO}^u, E_{ACB_{\rm I}}^u, E_{ACB_{\rm II}}^u]$ with $E_{SRO}^u = \sum\nolimits_{t = 1}^T {\sum\limits_{s \in {{\mathcal S}^u}} {\sum\limits_{i \in {\mathcal I}_s^u} {b_{i,s}^u{\rm tr}({\bm G_{i,s}})} } } $.

![Trend of achieved total slice utility vs. $m$.[]{data-label="fig:fig_total_utility_vs_gammath"}](total_slice_utility_vs_Csth.eps){width="3.1in"}

![Trend of achieved total slice utility vs. $d$.[]{data-label="fig:fig_total_utility_vs_Ds"}](total_slice_utility_vs_Ds.eps){width="3.1in"}

![Trend of achieved total slice utility vs. $\eta$.[]{data-label="fig:fig_total_utility_vs_eta"}](total_slice_utility_vs_eta.eps){width="3.1in"}

From these figures, the following observations can be achieved: i) the obtained utilities $\bar U$ of all algorithms decrease with an increasing $m$. This is because a great $m$ indicates that the accumulated IoT packets in the queue of each IoT device can be quickly emptied, and then a small $P_s(t)$ is obtained; ii) a great $D_s$ will reduce RRHs’ energy consumption. However, it also reduces the URLLC slice gain. Then, it may be hard to conclude the trend of $\bar U^u$ w.r.t. $D_s$ as the energy efficiency coefficient $\eta$ significantly affects the value of $\bar U^u$; iii) it is also uneasy to conclude the trend of $\bar U^u$ w.r.t. $\eta$. An increasing $\eta$ causes a decrease of RRHs’ energy consumption. Yet, the value of $\bar U^u$ is determined by the multiplier of $\eta$ and $E^u$; iv) the SRO-ACB$_{\rm II}$ algorithm may perform better than the SRO algorithm. However, the performance of the other access-control-based algorithm, SRO-ACB$_{\rm I}$, is slightly worse than the SRO algorithm. Besides, it cannot ensure that the $\bar U^I$ obtained by the access-control-based algorithms are always higher than that of SRO. At sometimes, access control schemes may drag down the utility of the mIoT service.

To sum up, in the case of service multiplexing, RA control schemes for alleviating signal interference and enhancing mIoT slice utility may be preferred for mIoT slices. However, considering both the CAPEX and the improvement of slice utility, RA control schemes should be carefully designed and employed because some RA control schemes may worsen the mIoT and even the total slice utilities.

Conclusion
==========

In this paper, we revisited the frame and minislot structure of a RAN slicing system to admit more IoT devices and proposed a queue evolution model to analyze the RACH of a randomly chosen IoT device. Based on the analysis result, we derived the closed-form expression of the RA success probabilities of the device with unrestricted access control scheme and ACB access control scheme. Next, we formulated the RAN slicing for mIoT and bursty URLLC service multiplexing as an optimization problem to optimally orchestrating RAN resources for mIoT slices and URLLC slices, and efficient mechanisms such as SAA and ADMM were exploited to mitigate the optimization problem. Simulation results showed that RA control schemes should be carefully designed and employed in the case of service multiplexing.

Proof of Lemma 1
----------------

For the origin RRH, the LT of its aggregate interference from interfering IoT devices in $s \in {\mathcal{S}^I}$ can be derived as $$\label{eq:LT_intra_interference}
\begin{array}{*{20}{l}}
{{{\mathcal L}_{{{\mathcal I}_s}(t)}}({\varpi _s}) = {E_{{{\mathcal I}_s}(t)}}\left[ {{e^{ - {\varpi _s}{{\mathcal I}_s}(t)}}} \right]}\\
= {E_{{{\mathcal I}_s}(t)}}\left[ {\exp ( { - {\varpi _s}\sum\limits_{{u_{m,s}} \in {\Phi _s}\backslash \{ o\} } {{\mathbbm 1}({p_m}||{d_m}|{|^{ - \varphi }} = {\rho _o}) \times } } } \right.\\
\text{ } \left. {\left. {{\mathbbm 1}({N_{a,s}}(t) > 0){\mathbbm 1}({f_m} = {f_o}){\rho _o}{h_m}} \right)} \right]
\\
\mathop  = \limits^{(a)} {E_{{\Phi _s}}}\left[ {{\prod _{{u_{m,s}} \in {\Phi _s}\backslash \{ o\} }}{E_{{h_m}}}\left[ {\exp \left( { - {\varpi _s} \times } \right.} \right.} \right.\\
\text{ } \left. {\left. {{\mathbbm 1}({p_m}||{u_{m,s}}|{|^{ - \varphi }} = {\rho _o}){\mathbbm 1}({N_{a,s}}(t) > 0){\mathbbm 1}({f_m} = {f_o}){\rho _o}{h_m}} \right)} \right] \\
{\mathop  = \limits^{(b)} \sum\limits_{n = 0}^\infty  {P\{ |{Z_s}| = n\} \prod\limits_{{u_{m,s}} \in {Z_s}} {{E_{{h_m}}}\left[ {{e^{ - {\varpi _s}{\rho _o}{h_m}}}} \right]} } }\\
{\mathop  = \limits^{(c)} P\{ |{Z_s}| = 0\}  + \sum\limits_{n = 1}^\infty  {P\{ |{Z_s}| = n\} {{\left( {\frac{1}{{1 + {\varpi _s}{\rho _o}}}} \right)}^n}} }\\
{\mathop  = \limits^{(d)} {{\tilde P}_{{X_s}}}\{ {X_s} = 1\}  + \left\{ {\sum\limits_{n' = 0}^\infty  {{{\tilde P}_{{X_s}}}\{ {X_s} = n'\} {{\left( {\frac{1}{{1 + {\varpi _s}{\rho _o}}}} \right)}^{n'}}}  - } \right.}\\
\text{ } {\left. {\sum\limits_{n' = 0}^1 {{{\tilde P}_{{X_s}}}\{ {X_s} = n'\} {{\left( {\frac{1}{{1 + {\varpi _s}{\rho _o}}}} \right)}^{n'}}} } \right\}(1 + {\varpi _s}{\rho _o})}
\end{array}$$ where $\varpi_s = \frac{\theta_s^{th}}{\rho_o}$, $Z_s$ denotes the set of interfering IoT devices in mIoT slice $s$, $X_s$ represents the number of active IoT devices associated with the origin RRH in $s$. According to the conclusion of Lemma 1 in [@yu2013downlink], the probability mass function (PMF) ${\tilde P}_{X_s} \{ X_s = n'\}$ can be written as $$\label{eq:intra_interference_proba}
    {\tilde P}_{X_s}\{ X_s = n'\}  = \frac{{{{3.5}^{3.5}}\Gamma (n' + 3.5){{(\frac{{{P_{nr,s}}(t)P_{ne,s}(t){\lambda _{s}^I}}}{{{\lambda _R}\xi F_s}})}^{n' }}}}{{\Gamma (3.5)(n')!{{(\frac{{P_{nr,s}}(t){P_{ne,s}(t){\lambda _{s}^I}}}{{{\lambda _R} \xi F_s}} + 3.5)}^{n' + 3.5}}}}$$ with $\Gamma(\cdot)$ being the gamma function. Besides, in (\[eq:LT\_intra\_interference\]), (a) follows from the i.i.d. distribution of $h_m$ and its further independence from the Poisson point process $\Phi_s$; (b) follows from the expectation of a discrete random variable; (c) follows from the LT over $h_m$; (d) follows from the fact that the number of active IoT device in a specific Voronoi cell is one more than the number of active interfering IoT devices in this cell.

From (\[eq:intra\_interference\_proba\]), we can deduce that $X_s$ ($s\in \mathcal{S}^I$) is a gamma¨CPoisson random variable with $X_s \sim {\rm gamma-Poisson}(\alpha_s, 3.5)$ and $\alpha_s = \frac{{{P_{nr,s}}(t){P_{ne,s}}(t){\lambda _{s}^I}}}{{3.5{\lambda _R}\xi F_s}}$.

For a gamma¨CPoisson random variable $X_s \sim {\rm gamma-Poisson}(\alpha, \beta)$, the following expression holds: $E[e^{X_s}] = (1+\alpha-\alpha e)^{-\beta}$. Thus, we can rewrite (\[eq:LT\_intra\_interference\]) as (\[eq:LT\_interference\_expression\]). This completes the proof.

Proof of Lemma 2
----------------

As new endogenous packet arrivals in any IoT device at each minislot $t$ is modelled as a Poisson distribution, the departure process of packets can be regarded as an approximated thinning process of new arrivals, where the thinning factor is related to the RA success probability. The number of accumulated packets in the queue of any IoT device can then be approximated as a Poisson distribution with intensity $\mu_{a,s}^{t}$ ($s\in \mathcal
S^I$) after the thinning process in a specific minislot $t$ ($t > 1$) [@jiang2018random]. Thus, we can derive the expression of $\mu_{a,s}^{t}$ ($t > 1$) via combining with the following facts

-   [**Fact 1:**]{} the accumulated packets during the $t-1$-th minislot will contribute to the accumulated packets at the $m$-th minislot.

-   [**Fact 2:**]{} the arrival packets during the $t-1$-th minislot will also contribute to the accumulated packets in the queue of an IoT device at the $m$-th minislot.

-   [**Fact 3:**]{} an IoT device can send packets only if its preamble is successfully transmitted.

-   [**Fact 4:**]{} at the same minislot, the new packet arrival process and the packet accumulated process are independent.

Similar as the Theorem 2 in [@jiang2018random], we can infer that at the $2^{\rm nd}$ minislot, for all $s \in \mathcal{S}^I$, $\mu_{a,s}^{2}$ depends on the intensity of new packet arrivals $\mu_{w,s}^{1}$ and the probability $P_{s}^{1}$ of a randomly selected IoT device at the $1^{\rm st}$ minislot, which is given by $$\label{eq:mu_acc_2_1}
    \begin{array}{l}
\mu _{a,s}^2 = \mu _{w,s}^1 - x_s P_s^1\left( {1 - {e^{ - \mu _{w,s}^1}}} \right)
\end{array}$$

The detailed proof of (\[eq:mu\_acc\_2\_1\]) is omitted for brevity, and a similar proof can be found in the proof section of Theorem 2 in [@jiang2018random].

Considering that $\mu_{a,s}(t)$ is non-negative at each minislot $t$, we have $$\label{eq:mu_acc_2}
    \begin{array}{*{20}{l}}
{\mu _{a,s}^2 = }
\end{array}{\left[ {\mu _{w,s}^1 - {x_s}P_s^1\left( {1 - {e^{ - \mu _{w,s}^1}}} \right)} \right]^ + }$$

Then, according to the definition of non-empty probability and the Poisson approximation, the non-empty probability of a randomly selected IoT device in mIoT slice $s \in \mathcal{S}^I$ at the $2^{\rm nd}$ minislot can be approximated as $$\label{eq:non_empty_prob_2_expression}
    P_{ne,s}^2 = 1 - {e^{ - \mu _{a,s}^2}}$$

At the $3^{\rm rd}$ minislot, the intensity of accumulated data packets in the queue of a randomly selected IoT device can be derived as the following $$\label{eq:mu_acc_3}
    \begin{array}{l}
\mu _{a,s}^3 = P_s^2\left( {\sum\nolimits_{n = 1}^\infty  {( {{{[n - x_s]}^ + } \sum\nolimits_{z = 0}^n {{P_{N_{w,s}^2}}(z){P_{N_{a,s}^2}}(n - z)} } )} } \right)\\
\quad + (1 - P_s^2)\left( {\sum\nolimits_{n = 1}^\infty  {n  \sum\nolimits_{z = 0}^n {{P_{N_{w,s}^2}}(z){P_{N_{a,s}^2}}(n - z)} } } \right)\\
\mathop  = \limits^{(a)} P_s^2\left [ {\sum\nolimits_{n = 1}^\infty  {\sum\nolimits_{z = 0}^n {\frac{{{{\left( {\mu _{w,s}^2} \right)}^z}{e^{ - \mu _{w,s}^2}}}}{{z!}}\frac{{{{\left( {\mu _{a,s}^2} \right)}^{n - z}}{e^{ - \mu _{a,s}^2}}}}{{(n - z)!}}}  \times n}  - } \right.\\
\quad \left. {x_s\sum\nolimits_{n = 1}^\infty  {\sum\nolimits_{z = 0}^n {\frac{{{{\left( {\mu _{w,s}^2} \right)}^z}{e^{ - \mu _{w,s}^2}}}}{{z!}}\frac{{{{\left( {\mu _{a,s}^2} \right)}^{n - z}}{e^{ - \mu _{a,s}^2}}}}{{(n - z)!}}} } } \right ]^+ + \\
\quad (1 - P_s^2)\sum\nolimits_{n = 1}^\infty  {\sum\nolimits_{z = 0}^n {\frac{{{{\left( {\mu _{w,s}^2} \right)}^z}{e^{ - \mu _{w,s}^2}}}}{{z!}}\frac{{{{\left( {\mu _{a,s}^2} \right)}^{n - z}}{e^{ - \mu _{a,s}^2}}}}{{(n - z)!}}}  \times n} \\
\mathop  = \limits^{(b)} \left [ \mu _{w,s}^2 + \mu _{a,s}^2 - x_sP_s^2\left( {1 - {e^{ - \mu _{w,s}^2 - \mu _{a,s}^2}}} \right) \right ]^+
\end{array}$$ where $P_{N_{w,s}^2}$ and $P_{N_{a,s}^2}$ represent the PMFs of new arrival packets and accumulated packets at the $2^{\rm nd}$ minislot, respectively. Besides, (a) follows from the fact: for any two independent Poisson distribution $\Phi_{X_1}$ and $\Phi_{X_2}$, $P_{X_1,X_2}({X_1} + {X_2} = x) = \sum\nolimits_{y = 0}^x {P_{X_1}({X_1} = y)P_{X_2}(X_2 = x - y)} $; (b) holds as $\Phi_{X_1,X_2}$ is a two dimensional Poisson distribution with an intensity $\lambda_{X_1} + \lambda_{X_2}$, and $\sum\nolimits_{x = 1}^\infty  {P_{X_1, X_2}({X_1} + {X_2} = x)}  = 1 - P_{X_1,X_2}({X_1} + {X_2} = 0)$.

Similarly, we have $$\label{eq:non_empty_prob_3}
    P_{ne,s}^3 = 1 - {e^{ - \mu _{a,s}^3}}$$

When $t > 3$, since the accumulated packets evolution model of the queue of any IoT device is the similar as that at $t = 3$, we can directly extend the conclusion obtained at $t = 3$ to that at $t > 3$. Therefore, we can obtain the closed-form expression of $\mu_{a,s}^t$ for all $s \in \mathcal{S}^I$ at $t > 1$ with $$\label{eq:mu_acc_m_expression_proof}
    \mu _{a,s}^t = \left [ \mu _{w,s}^{t - 1} + \mu _{a,s}^{t - 1} - x_sP_s^{t - 1}\left( {1 - {e^{ - \mu _{w,s}^{t - 1} - \mu _{a,s}^{t - 1}}}} \right) \right ]^+$$ and $$\label{eq:non_empty_prob_m_expression_proof}
    P_{ne,s}^t = 1 - {e^{ - \mu _{a,s}^t}}$$ This completes the proof.

Proof of Lemma 4
----------------

The $2^{\rm nd}$ degree Taylor expansion of $\tilde P_{m}^{(k)}$ at the local point $\bm \omega_{m}^{(k,q)}$ is $$\tilde P_{2,m}^{(k)} = \sum\limits_{j = 0}^2 {\frac{1}{{j!}}{{\left[ {\sum\limits_{s \in {\mathcal S}^I} {\left( {{\omega _{sm}} - \omega _{sm}^{(k,q)}} \right)\frac{\partial }{{\partial \omega _{sm}^{(k,q)}}}} } \right]}^j}{\tilde P_m^{(k)}}|_{\bm \omega_m^{(k,q)}}}$$

The $3^{\rm rd}$ degree Taylor expansion of $\tilde P_{m}^{(k)}$ at $\bm \omega_{m}^{(k,q)}$ must be more accurate than $\tilde P_{2,m}^{(k)}$ with $$\tilde P_{3,m}^{(k)} = \tilde P_{2,m}^{(k)} + {\frac{1}{{3!}}{{\left[ {\sum\limits_{s \in {\mathcal S}^I} {\left( {{\omega _{sm}} - \omega _{sm}^{(k,q)}} \right)\frac{\partial }{{\partial \omega _{sm}^{(k,q)}}}} } \right]}^3}{\tilde P_m^{(k)}}|_{\bm \omega_m^{(k,q)}}}$$

Since the error of $\tilde P_{2,m}^{(k)}$ is no greater than the maximum difference between $\tilde P_{3,m}^{(k)}$ and $\tilde P_{2,m}^{(k)}$, we have $$\label{eq:error_R2}
    R_{2}(\bm \omega_m) = \max \{ {\frac{1}{{3!}}{{[ {\sum\limits_{s \in {\mathcal S}^I} {( {{\omega _{sm}} - \omega _{sm}^{(k,q)}} )\frac{\partial }{{\partial \omega _{sm}^{(k,q)}}}} } ]}^3}{\tilde P_m^{(k)}}|_{\bm \omega_m^{(k,q)}}} \}$$

In (\[eq:error\_R2\]), $\bm \omega_{m}^{(k,q)}$ is a constant vector, the $\max$ operation will not affect the constant vector and the vector $\bm \omega_{m}$. For any $s\in \mathcal{S}^I$, the maximum value obtainable by $\frac{{{\partial ^3}{\tilde P_m^{(k)}}|_{\bm \omega_m^{(k,q)}}}}{{\partial \omega _{sm}^{3(k,q)}}}$ will not exceed the greatest value of that derivative in the interval $[\omega_{sm}^{lb}, S_{sm}^{\star}]$. Additionally, the maximum value of $\frac{{{\partial ^3}{P_m}|_{\bm \omega_m^{(k,q)}}}}{{\partial \omega _{sm}^{3(k,q)}}}$ will generally occur at one of the endpoints of the interval $[\omega_{sm}^{lb}, S_{sm}^{\star}]$. Therefore, we obtain (\[eq:error\_bound\_lemma\]). This completes the proof.

Proof of Lemma 5
----------------

For all $i \in {\mathcal{I}^u}$, $s \in {\mathcal{S}^u}$, $m \in {\mathcal{M}}$, a feasible way of proving that ${\rm rank}(\bm G_{i,sm}) \le 1$ is to utilize the Lagrange method. However, owing to the complicated expression of $W^{u}(\bm r_m)$ w.r.t. $\bm G_{i,sm}$ it will be uneasy to do that. Fortunately, we find that the proof can be conducted if a family of auxiliary variables is introduced.

For the constraint (\[eq:SAA\_admm\_problem\]d), if we introduce the auxiliary variables $\{\nu_{i,sm}\}$ and let $$\label{eq:linear_snr_tau_ism}
    \frac{{\rm{tr}}({\bm H_{i,sm}}{{\bm G}_{i,sm}})}{\phi \sigma _{i,s}^2} \ge \nu _{i,sm},\forall i \in {\mathcal I}_s^u,s \in {{\mathcal S}^u},m \in {\mathcal{M}}$$ then (\[eq:SAA\_admm\_problem\]d) is equivalent to $$\label{eq:}
    \sum\limits_{s \in {{\cal S}^I}} {(1+\alpha_g)\omega _{sm}(\bar t)}  + W^u(\bm f_m)  \le W, \text{ } {\rm and} \text{ } {\rm (\ref{eq:linear_snr_tau_ism})}$$ where $\bm f_m = \{f_{i,sm};i \in {\mathcal{I}_s^u}, s \in {\mathcal{S}^u}\}$ and $$\begin{array}{*{20}{l}}
{{f_{i,sm}} = \frac{{L_{i,s}^u}}{{{{\log }_2}(1 + {\nu _{i,sm}})}} + \frac{{{{({Q^{ - 1}}(\beta))}^2}}}{{2\log _2^2(1 + {\nu _{i,sm}})}}}\\
{\qquad \quad  + \frac{{{{({Q^{ - 1}}(\beta))}^2}}}{{2\log _2^2(1 + {\nu _{i,sm}})}}\sqrt {1 + \frac{{4L_{i,s}^u{{\log }_2}(1 + {\nu _{i,sm}})}}{{{{({Q^{ - 1}}(\beta))}^2}}}} }
\end{array}$$

We omit the proof of the equivalence as a similar proof can be found in the proof section of constraints’ equivalence in [@How2019yang].

The partial Lagrangian function of (\[eq:SCA\_bandwidth\_beamforming\_problem\]) can be written as $$\label{eq:lagrangian_func}
\begin{array}{l}
L( \ldots ) = \sum\limits_{s \in {{\mathcal S}^u}} {\sum\limits_{i \in {\mathcal I}_s^u} {\left[ {\frac{{\tilde \rho \eta }}{M}{\rm{tr}}({\bm G_{i,sm}}) - {{\bar \mu }_{i,sm}}\frac{{{\rm{tr}}({\bm H_{i,sm}}{\bm G_{i,sm}})}}{{\phi \sigma _{i,s}^2}} + } \right.} } \\
\qquad \qquad \left. {\sum\limits_{j \in {\mathcal J}} {{{\bar \lambda }_{jm}}{\rm{tr}}(b_{i,sm}^{u(k,q)}{\bm Z_j}{\bm G_{i,sm}})}  - {{\bm {\bar X}}_{i,sm}}{\bm G_{i,sm}}} \right]
\end{array}$$ where $\bar \lambda_{jm}$, $\bar \mu_{i,sm}$, and $\bar {\bm X}_{i,sm}$ are Lagrangian multipliers corresponding to constraints (\[eq:SAA\_admm\_problem\]c), (\[eq:linear\_snr\_tau\_ism\]) and (\[eq:SAA\_admm\_problem\]e). Besides, only terms related to ${\bm G}_{i,sm}$ are included in this function for brevity. According to the Karush-Kuhn-Tucker (KKT) conditions, the necessary condition for obtaining the optimal matrix power at the $(k,q)$-th iteration ${\bm G}_{i,sm}^{(k,q)\star}$ is given by $$\label{eq:KKT_condition_2}
    \begin{array}{l}
\frac{{\partial L( \ldots )}}{{\partial {\bm G}_{i,sm}^{(k,q)\star}}} =  \frac{{\tilde \rho \eta }}{M}{\bm I_{i,sm}} + \frac{\bar \mu_{i,sm}{{\bm H_{i,sm}}}}{{\phi \sigma _{i,s}^2}} - \\
\qquad \qquad \text{ } \sum\limits_{j \in {\mathcal J}} {{{\bar \lambda }_{jm}}{b_{i,sm}^{u(k,q)}\bm Z_j}} - {{\bm X}_{i,sm}} = 0
\end{array}$$ where ${\bm I_{i,sm}} \in \mathbb{R}^{JK \times JK}$ is an identity matrix.

Then, we can conclude that ${\rm rank}({\bm X}_{i,sm}) \ge JK - 1$. The reasons are i) ${{\bar \lambda }_{jm}}$, $b_{i,sm}^{u(k,q)}$, and $\bar \mu _{i,sm}$ are nonnegative and the matrix $\bm I_{i,sm}$ is full rank; ii) ${\rm rank}(\bm H_{i,sm}) \le 1$.

Next, according to the complementary slackness condition, we have $$\label{eq:complementary_slack_condition_2}
    {{\bm X}_{i,sm}}{\bm G}_{i,sm}^{(k,q)\star} = 0$$

Based on (\[eq:complementary\_slack\_condition\_2\]) and the rank result of $\bm X_{i,sm}$, we can conclude that ${\rm rank}(\bm G_{i,sm}^{(k,q)\star}) \le 1$. This completes the proof.

[^1]: P. Yang and T. Q. S. Quek are with the Information Systems Technology and Design, Singapore University of Technology and Design, 487372 Singapore.

    X. Xi, J. Chen, and X. Cao are with the School of Electronic and Information Engineering, Beihang University, Beijing 100083, China, and also with the Key Laboratory of Advanced Technology, Near Space Information System (Beihang University), Ministry of Industry and Information Technology of China, Beijing 100083, China.

    D. Wu is with the Department of Electrical and Computer Engineering, University of Florida, Gainesville FL 32611 USA.

[^2]: We do not show all channels in this figure as the detailed research of the physical layer supporting the mIoT service is out of the scope of this paper.

[^3]: Owing to the truncated channel inversion power control, not all of the IoT devices in mIoT slices can communicate in the uplink when the cutoff threshold is relatively high [@elsawy2014stochastic]. However, this paper assumes that the transmit power of each IoT device is large enough such that the IoT device will not experience preamble outage resulting from the insufficient power.
