Botao Yang
email: botao.yang19@student.xjtlu.edu.cn

Ye Liu
email: ye.liu02@xjtlu.edu.cn

Chung Shue Chen
email: chung_shue.chen@nokia-bell-labs.com

Uplink Scheduling in a NOMA-enabled Single-Cell Wireless Network Using Simulated Annealing

Keywords: Non-orthogonal multiple access (NOMA), resource scheduling, maximum completion time, serial collaborative optimization, simulated annealing (SA)

By allowing multiple users to transmit using the same frequency band at the same time, non-orthogonal multiple access (NOMA) can support more users as compared to orthogonal multiple access (OMA) given a fixed amount of time-frequency resources. In this paper, we study the resource allocation problem in the uplink of a single-cell network when NOMA is enabled, where the maximum completion time of serving all connected users is to be minimized. While such an objective function minimizes the required time to serve uplink users, the resource allocation problem is NP-hard. We propose a serial collaborative optimization framework based on simulated annealing (SA) to search for the optimal user pairing and scheduling solution. Simulation studies show that the proposed algorithm for NOMA scheduling can reduce the maximum completion time by more than 30% when compared against OMA scheduling and random NOMA user clustering.

I. INTRODUCTION

In 4G/5G wireless cellular systems, orthogonal frequencydivision multiple access (OFDMA) is used as the major radio access technology. OFDMA is an orthogonal multiple access (OMA) scheme that avoids intra-cell interference. Despite its ease of implementation, OMA suffers from low spectral efficiency when some bandwidth resources are allocated to users with poor channel conditions [START_REF] Islam | Powerdomain non-orthogonal multiple access (NOMA) in 5G systems: Potentials and challenges[END_REF]. With the rapid growth of demand for mobile data network capacity and transmission speed, OMA based cellular systems are unable to keep up with the applications in the next-generation radio access technology [START_REF] Wei | On the performance gain of NOMA over OMA in uplink communication systems[END_REF], , and non-orthogonal multiple access (NOMA) is considered to be a promising technique to overcome the limitations of OMA and improve system performance [START_REF] Ding | Application of non-orthogonal multiple access in LTE and 5G networks[END_REF]. In contrast to OMA schemes, NOMA allows multiple users to be superimposed on the same time-frequency resource block, where successive interference cancellation (SIC) can be used to separate signals from different users. As a result, the bandwidth allocated to the users with poor channel conditions can still be used by the users with strong channel conditions, which can significantly improve the system's spectral efficiency [START_REF] Tabassum | Uplink vs. downlink NOMA in cellular networks: Challenges and research directions[END_REF].

With higher sum capacity and enhanced spectrum efficiency, NOMA provides support for more users' communication needs at the same time, thus having the potential to improve quality-of-service (QoS) for a variety of communication applications. For delay insensitive services such as email, rate constraints are relatively relaxed for individual users and throughput (i.e., the sum of data rates delivered to all terminals in a network) is the most broadly used criterion to evaluate system performance. Over the past few years, this topic has been already thoroughly explored by academia and industry [START_REF] Ding | Impact of user pairing on 5G nonorthogonal multiple-access downlink transmissions[END_REF]- [START_REF] Mishra | Downlink connection density maximization for NB-IoT networks using NOMA with perfect and partial CSI[END_REF]. On the other hand, the explosive growth of the Internet in recent years has also created many applications that are sensitive to packet delays from end-to-end. In these cases, each user must transmit a certain amount of information within a certain period and maximizing the throughput is no longer an appropriate strategy [START_REF] Altintast | On a packet scheduling mechanism for supporting delay sensitive applications on high speed networks[END_REF]. Given a number of requests from multiple devices to transmit data packets of varying lengths, efficient scheduling scheme is required to avoid long queuing at the buffers of congested links.

In terms of delay optimization, the maximum completion time (i.e., the makespan, denoted by C max ) is a widely adopted metric. Recently, completion time reduction in NOMAenabled systems has received attention in the research community. Specifically, the uplink NOMA scheduling problem of minimizing C max is proven to be NP-hard [START_REF] Mohsenivatani | Completion-time-driven scheduling for uplink NOMAenabled wireless networks[END_REF], leading to high computational complexity in conventional iterative methods. To achieve completion time minimization, [START_REF] Wang | Completion time minimization in NOMA systems: Learning for combinatorial optimization[END_REF] proposed a learning-based hybrid algorithm to jointly optimize time slot and transmit power allocation. However, this work assumes that there is only one frequency channel available, which does not reflect most practical communication scenarios.

To solve this intractable problem within affordable time and complexity, heuristic algorithms with performance guarantees are preferred. Simulated annealing (SA) [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] is a metaheuristic approach which seeks for the global optimum of a given optimization problem. Compared to conventional iterative methods, the advantages of SA include its easy implementation and possibility of reaching the global optimal without being trapped in a local optimum [START_REF] Granville | Simulated annealing: A proof of convergence[END_REF]. SA scheme can provide satisfactory result with a relatively low number of iterations [START_REF] Szu | Nonconvex optimization by fast simulated annealing[END_REF], [START_REF] Chen | Joint optimization of radio resources in small and macro cell networks[END_REF], which makes it suitable for hard optimization problem.

Aiming at minimizing the maximum completion time of serving uplink users in a single-cell NOMA-enabled network, we formulate the scheduling problem as a combinatorial optimization problem with a discrete strategy set. The correlation between user clustering (or, equivalently, user grouping) and frequency channel allocation is taken by a serial collaborative optimization. In the proposed scheme, the optimizations of Fig. 1. User grouping in an uplink NOMA-based system, where the base station is located at the middle point. user clustering and channel assignment scheduling are carried out in a separate manner, which significantly decreases the computational complexity. Both steps are optimized by SA but with independent parameter conditions. Numerical results show that our proposed scheme reduces more than 30% of the maximum completion time against heuristic OMA on average in the considered simulation settings.

The rest of the paper is organized as follows. In Section II, we describe the system model and problem formulation. In Section III, we present the proposed solution and optimization algorithm. In Section IV, we show the simulation result and comparison. Finally, we conclude the paper in Section V.

II. SYSTEM MODEL AND PROBLEM FORMULATION

Consider the uplink of a single-cell NOMA network, where K users request to transmit data packets of varying lengths to the base station (BS), as shown in Fig. 1. The frequency resource is equally divided into F frequency channels, where each channel has a bandwidth of W Hz. The BS is equipped with a single antenna and all the K users are single-antenna devices. Each user has a maximum available transmission power of P max . Also, each user k transmits a single data packet of L k bits in length with power p k .

The wireless links experience independent and identically distributed block Rayleigh fading and additive white Gaussian noise (AWGN), where the power spectral density of the AWGN is denoted as N 0 . We consider power-domain NOMA, where SIC at the receiver is used to decode signals from users that are sent in the same time-frequency block. To minimize the complexity of SIC, we assume that each frequency channel can be simultaneously occupied by at most two users. Hence, at most K/2 NOMA clusters can be formed, where each NOMA cluster contains two users that concurrently transmit to the serving BS via the same frequency channel. Without loss of generality, K is assumed an even number, for otherwise we can add one dummy user with 0 bit to transmit.

The channel gain between user k and the BS via channel f ∈ F is denoted by h f k , where F is the set of frequency channels and the size of F, i.e., |F|, is equal to F . Specifically,

h f k ≜ g f k / d α k , where g f k ∼ CN (0, 1)
is the Rayleigh fading coefficient, d k is the distance between user k and the BS, and α is the path loss exponent. For simplicity, we assume that for any specific user k, the channel coefficient remains the same regardless of the frequency channel selected [START_REF] Mohsenivatani | Completion-time-driven scheduling for uplink NOMAenabled wireless networks[END_REF], i.e. 1 ,

|h 1 k | 2 = |h 2 k | 2 = • • • = |h F k | 2 = |h k | 2 . ( 1
)
Users are sorted and indexed in ascending order according to the channel conditions, i.e.,

|h 1 | 2 ≤ |h 2 | 2 ≤ • • • ≤ |h K | 2 . ( 2
)
It is assumed that that perfect channel state information (CSI) is available at BS and the channel coefficients are quasi-static. When two users i and j transmit in NOMA mode, we form a NOMA cluster (i, j) and perform scheduling based on this NOMA cluster. The BS receives a linear superposition of signals from these two users, i.e.,

y f = h i √ p i x i + h j √ p j x j + w, (3)
where y f is the received signal on the f -th frequency channel, x k is the normalized unit-power signal transmitted by user k such that E[x 2 k ] = 1, and w denotes the noise received per channel with power N 0 W . Suppose that user i is the strong user with a higher SNR, i.e., p i |h i | 2 > p j |h j | 2 . According to the principle of SIC, the strong user i is decoded first, treating user j's signal as interference. After the strong user's signal is correctly decoded, the BS can regenerate the strong user's signal and remove that from the received composite signal. Then, the BS can decode the weak user's signal without cochannel interference from the other user. For the first decoded user (i.e., strong user i), the achievable data rate is given by

R SIC,j i = W log 2 1 + p i |h i | 2 p j |h j | 2 + N 0 W = W log 2 1 + p i γ i p j γ j + 1 , (4)
where

γ i ≜ |h i | 2 /(N 0 W ).
Next, the second decoded user (weak user j) can enjoy the interference-free transmission so the achievable rate of user j equals to R OMA j , i.e.,

R OMA j = W log 2 1 + p j |h j | 2 N 0 W = W log 2 (1 + p j γ j ). (5)
With the above, the required transmission time of user j on any frequency channel f in OMA is given by

τ j = L j R OMA j . (6)
Two users within a NOMA cluster start transmitting at the same time. As illustrated in Fig. 2, the time required to finish the transmission of the two-user NOMA cluster (i, j) depends on which user finishes sending its packet first, i.e.,

u i,j =    τ j , if Li R SIC,j i < τ j (case 1), τ j + Li-R SIC,j i τj R OMA i , otherwise (case 2). (7)
Note that in the above case 1, the required transmission time of the NOMA cluster is the same as that of the weak user j, since user j may have a long packet to send and/or the OMA data rate of user j is low. In case 2, the transmission rate of the strong user i has the opportunity to be from R SIC,j i to R OMA i when the weak user j finishes transmission. While some signaling may be needed to realize such a data-rate increase, here we use the model in case 2 of ( 7) to study the potential gain using NOMA in ideal situations, ignoring any potential overhead that may be necessary in practice.

It is evident from ( 4) and ( 5) that within a NOMA cluster (i, j), by letting the strong user i transmit at a power level p i = P max , the processing time of user i (τ i in OMA or, when using NOMA, Li R SIC,j i in case 1 or u i,j in case 2) will be minimized. While p j affects the processing time of users i and j and would be subject to optimization, the value of p j and the user clustering decision are intertwined. Specifically, the user clustering decision needs to take the users' power and channel gain profiles into account, such that finding a suitable weak user to form a cluster with a strong user is similar to adjusting the transmission power of the weak user in an existing NOMA cluster. Thus, we let p k = P max , ∀k, in the rest of the paper.

In the following, we introduce the two essential elements of the makespan-minimization problem with NOMA.

A. Cluster Establishment

Denote l ∈ L as a particular clustering scheme that groups the K users into θ(l) clusters, where L is the set of all feasible clustering methods and each cluster may contain either one or two users, i.e., K/2 ≤ θ(l) ≤ K. For any pair of users i and j, one observes that

τ i + τ j ≥ u i,j , (8)
where the proof of ( 8) is detailed in Appendix A. Equation [START_REF] Mishra | Downlink connection density maximization for NB-IoT networks using NOMA with perfect and partial CSI[END_REF] suggests that if users i and j are transmitting on the same frequency channel, then using NOMA to serve the two users will never prolong the completion time on that channel as compared to OMA. While it is possible to reduce the maximum completion time of the network if u i,j takes the form of case 2 and users i and j are scheduled on different channels2 , allowing a cluster to have a single user enlarges |L| and may increase the complexity of solving the NOMA scheduling problem. Encouraged by the benefit of NOMA as shown in [START_REF] Mishra | Downlink connection density maximization for NB-IoT networks using NOMA with perfect and partial CSI[END_REF], we assume that when using NOMA, each cluster contains two users and K/2 clusters are formed. Given a user grouping scheme, the processing time of every cluster can be determined by [START_REF] Jang | Throughput maximization for uplink IoT NOMA systems with efficient resource allocation[END_REF]. We can index the clusters by m ∈ {1, 2, • • • , K/2} and sort them according to the descending order of the time required to finish the transmission for each cluster, denoted by u m , such that

u 1 ≥ u 2 ≥ • • • ≥ u K/2 . (9)

B. Spectrum and Time Resource Allocation

After the user clustering is determined, the remaining scheduling problem belongs to the class of parallel machine scheduling that is NP-hard [START_REF] Conway | Theory of Scheduling[END_REF], where each NOMA cluster requires a single transmission on one of the F identical but independent frequency channels.

Denote s ∈ S as a particular scheduling scheme that allocates the two-user clusters to the F frequency channels, where S is the set of all feasible scheduling schemes. A scheduling scheme s can be represented by a size-K 2 tuple whose elements indicate the destination channel for each cluster. Specifically, the m-th cluster is labeled to transmit using the s(m)-th channel, i.e.,

1 ≤ s(m) ≤ F, m ∈ {1, 2, • • • , K/2}. ( 10
)
We impose the following constraints when performing the scheduling:

1) A channel, when scheduled to carry uplink data signals from users, will keep busy until all scheduled transmissions are completed. In other words, we ignore the impact of any packet header or signaling that may exist in a particular communication protocol. 2) No cluster is allowed to occupy more than one frequency channel at the same time. 3) Once the transmission for a user starts in a particular frequency channel, the scheduler is not allowed to interrupt it until its completion.

While removing the second and the third constraints may provide further space for minimizing the maximum completion time, a more complicated communication protocol and scheduling algorithm would be needed if those constraints are not imposed. We leave the study of uplink NOMA scheduling without the above constraints to a future work. Given the above settings, the time needed to complete the transmission at channel f is given by the sum of the completion times of the clusters that are allocated to the same channel, i.e.,

C f = m∈M f u m , (11)
where M f ≜ {m | s(m) = f } is the set of clusters that are allocated to channel f . Assume that the transmissions at all the channels start at time zero. The maximum completion time of the NOMA system is defined as the time required to finish transmission in all the channels, as depicted by Fig. 3, i.e.,

C max = max f ∈F {C f }. ( 12
)
Our goal is to find the optimal user clustering and channel schedules that minimize the maximum completion time of the

III. OPTIMIZATION AND SOLUTION

A. Complexity Analysis

The problem in ( 13) is of combinatorial nature with a finite discrete strategy set. The total number of user clustering schemes is given by

|L| = K! 2 K/2 (K/2)! = Π K/2 i=1 (K/2 + i) 2 K/2 > (K/4) K/2 , ( 14
)
where the inequality is true because in

Π K/2 i=1 (K/2+i) 2 K/2
, there are K/2 terms in the numerator that are no smaller than K/2. As for the total number of channel scheduling schemes |S|, note that the transmission order of the clusters in a channel does not affect the completion time. Therefore, |S| can be calculated by the Stirling number of the second kind, i.e., S(K/2, F ) [START_REF] Ilie | Computation of Stirling numbers and generalizations[END_REF]. A lower bound on S(K/2, F ) can be obtained as

S(K/2, F ) > F K/2-F , (15)
where Appendix B details the proof. As a result, the size of the strategy set of our problem is lower bounded as

|L × S| > (K/4) K/2 F K/2-F . (16)
Equation ( 16) suggests that the size of the strategy set for NOMA optimization scales up exponentially as the number of channels and serving users in the system increase. Therefore, exhaustive search for optimal solution is computationally expensive even for moderate K and F . One class of approximation algorithms is stochastic search with iterative improvement [START_REF] Hoos | Solving combinatorial auctions using stochastic local search[END_REF]. This local search method works by continuously exploring the neighborhood of a solution in order to find a better solution. When there is a gain, the current solution is replaced by the new one. However, the NP-hardness nature of a problem can limit the performance of conventional stochastic search because such algorithms can be trapped by local minima. Besides, purely random search is known to be slow for convergence and may not be practical for the application in our context.

B. Simulated Annealing 1) Preliminary:

The core mechanism of simulated annealing (SA) is based on the analogy of ideal crystal formation in thermodynamics. In the controlled annealing process, the disordered configuration of particles in high energy system gradually rearranges into regular crystalline solid state with low energy [START_REF] Mahady | Sum-rate maximization of NOMA systems under imperfect successive interference cancellation[END_REF]. SA algorithm mimics this metallurgical method to minimize the objective function by a series of guided probabilistic random search in the strategy set.

In general, a SA follows the procedures below: 1) Select a starting point x 0 and initialize from a high temperature setting. 2) Perturb the current placement to a new state. If the change is accepted, the current state is updated to the new one. This process is repeated in an iterative manner, until thermodynamic equilibrium is reached with current temperature.

3) The temperature is then lowered, before carrying out a new series of transformation (state transitions). 4) The above loop continues until a freezing point or preset termination condition is reached. Note that in an iteration, SA may not always select the next state that corresponds to a better solution as compared to the current state. In particular, the next state that exhibits as a worse solution than the current state can be chosen probabilistically [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. In the minimization problem, a next state that decreases the objective function will be accepted; on the other hand, a next state that increases the objective function also has a chance of being accepted with a probability determined by a given acceptance function.

Note that the variability of the candidate solution and randomness of new state acceptance decrease as the parameter called the temperature T drops (cooling). Specifically, 1) When T → ∞ (high-energy system), it behaves like a purely random walk that is always moving to a successor chosen uniformly at random from the entire strategy set. 2) When T → 0 (low-energy system), it degenerates into a deterministic iterative improvement scheme that refuses any moves towards worse objective values. Note that it can get stuck on a local minimum. 3) At an intermediate temperature, the SA algorithm intermittently accepts the transformations that increase the objective function. This provides for the system to be pulled out from a local minimum. 2) Sampling Method: Sampling method is the core of SA algorithm, which governs the behavior of iteratively updating the current state depending on the evaluation of candidate states. There are two mainstream sampling methods for the trial point selection and acceptance test [START_REF] Glass | A comparison of local search methods for flow shop scheduling[END_REF].

1) Gibbs Sampler Algorithm (GS): Denote N as the state set. At a temperature T , the sampler will evaluate all the possible candidate state x ∈ N and then selects a new state according to the following probability distribution:

p(x) = e -U (x)/T x ′ ∈N e -U (x ′ )/T (17)
where U (•) is the evaluation or cost function.

2) Metropolis-Hasting Algorithm (MH): This sampling algorithm randomly selects a candidate state x ′ from N and evaluates the change to the cost function:

∆U = U (x ′ ) -U (x). (18)
If the state transition leads to a decrease in the cost, it is unconditionally accepted. If the transition would cause an increase to the cost function, the new state will be accepted with the following probability, i.e.,

p accept (x ′ ) = e -∆U/T . ( 19
)
When |N | is large, the computational complexity of GS for each state update is far above the MH algorithm. Therefore, we use MH for our problem. We conduct multiple iterations at each temperature stage in MH in order to compensate its possible deficiency due to the locally incomplete sampling performed per iteration. 3) Cooling Schedule: Cooling scheme would affect the convergence speed and the quality of the SA optimization. Empirically, the temperature should decrease gradually by stages to enforce regularities of the system. It can be proven that SA converges asymptotically to the optimal solution. However, such convergence requires exponential time [START_REF] Granville | Simulated annealing: A proof of convergence[END_REF]. Hence, in practice, faster cooling schedules are used, with which SA does not guarantee to converge exactly to the global optimum and thus provides an approximate solution. According to the SA sampling mechanism, randomness of state selection and the tolerance for bad move would be decreased with temperature drop. Therefore, an appropriate cooling scheme should be able to initially search the space for global minima and later, as the system cools down, converge on a solution with proper accuracy [START_REF] Hajek | Cooling schedules for optimal annealing[END_REF]. Annealing too fast is likely to result in defects. On the other hand, annealing too slow will lead to a large number of iterations taken for convergence. Hence, to strike a balance between the convergence speed and optimality, we adopt the following geometric cooling scheme in our work:

T n+1 = βT n , T n = T 0 β n , 0 < β < 1, ( 20
)
where T 0 is the initial temperature, β is a tunable parameter which controls the temperature drop rate, and n is the index for the temperature cooling schedule.

The SA algorithm for our problem is shown in Algorithm 1. Specifically, with a fixed temperature value, a number of I iterations are conducted. The number of new state acceptance within each I iterations is recorded by η, and the frequency of new state acceptance, calculated as η • I -1 , is compared against a threshold value "minMove". If η • I -1 < minMove, it is considered that the steady state is reached (also known as the thermodynamic balance), and the SA algorithm terminates with the output x * . Otherwise, the temperature is reduced according to [START_REF] Hajek | Cooling schedules for optimal annealing[END_REF], and another I iterations will be performed using the new temperature. The SA algorithm will also terminate if a preset maximum number of iterations is reached.

C. Serial Collaborative Optimization

In our problem, every feasible solution consists of a user grouping scheme combined with a channel scheduling scheme, Algorithm 1 SA Optimization Algorithm (MH) 1: Initialization: 2: x ← x 0 starting state 3: T ← T 0 initial temperature 4: x * ← x the best solution found so far 5: repeat 6:

for i = 1 to I do end if end for 21: until η • I -1 < minMove 22: Output: x * which are correlated with each other, i.e.,

x = (l, s). (21)
Instead of treating user clustering and resource allocation as a joint union [START_REF] Sun | Jointly optimizing user clustering, power management, and wireless channel allocation for NOMA-based Internet of Things[END_REF], our proposed serial collaborative method decomposes the problem into two separate parts; see Fig. 4. Such a decomposition can significantly reduce the search space in each SA loop and, therefore, reduce the overall computational complexity. Note that both SA loops are implemented according to Algorithm 1 with different parameters. 1) SA Loop 1 for Cluster Establishment: In this step, we perform SA optimization for user grouping scheme to minimize the sum of the processing time of all clusters. A candidate solution can be generated by the following procedure:

1) Randomly select two clusters and split them to four users. 2) Reassemble them into two different clusters. Define N l as the neighborhood set of any clustering scheme l. The size of N l can be found as

|N l | = 2 × K/2 2 = K(K -2) 4 , (22)
where

K/2 2 = (K/2)(K/2-1) 2 .
The utility function is selected as the sum of the processing time of all clusters, i.e., 0.5K m=1 u m . While such a utility function is not equal to the objective function in [START_REF] Granville | Simulated annealing: A proof of convergence[END_REF], the utility function minimizes a measure of the processing time of all the clusters which will then impact the maximum completion time for any given schedule.

2) SA Loop 2 for Channel Allocation: With the user clustering obtained in SA loop 1, the processing time of every cluster is fixed. The remaining problem is how the tasks are distributed over the available frequency channels. In this step, we perform SA optimization for channel scheduling scheme to minimize the makespan (i.e., maximum completion time).

A candidate solution can be generated by randomly selecting two clusters in different frequency channels and swapping their position. The number of candidate solutions is maximal when clusters are equally distributed among the available frequency channels. Define N s as the set of neighborhood of a scheduling scheme s. Therefore, the size of N s has the following upper bound, i.e.,

|N s | ≤ F 2 × K 2F 2 ≤ F 2 K 2F + 1 2 = 1 8 (K + 2F ) 2 (1 - 1 F ).( 23
)
The utility function for SA loop 2 is selected as the maximum completion time, i.e., max f ∈F C f .

3) Initialization: In both SA loops, the state generation schemes allow that every single state can be reached from any other state in a number of moves. Therefore, the initial solution can be selected randomly to start the iterative search.

IV. SIMULATION AND NUMERICAL RESULTS

A. Simulation Experiment Design

To test performance, we produce the simulation results by randomly generating 1000 network topologies and then averaging the makespan from all the topologies. Specifically, in each topology, each user is randomly scattered within a bounded annular region centered at the BS, where the user-BS distance d k is between 100 meters and 1000 meters. The maximum allowed user transmission power is set to be P max = 23 dBm. The path loss exponent is set to α = 3. The packet size of each user is randomly chosen in ranges that are specified in Table I. Parameters related to the SA algorithm can be found in Table II.

In order to provide some benchmark results in comparison with the serial collaborative algorithm, the maximum completion time is also evaluated for random NOMA clustering (i.e., randomly select two users for each cluster) as well as conventional OMA scheme scheduled by the longest-processingtime-first (LPT) rule. The LPT rule is a heuristic method for the NP-hard identical machine scheduling problem [START_REF] Conway | Theory of Scheduling[END_REF]. The LPT rule works as follows:

1) Order the jobs (i.e., the users in this paper's context) by descending order of their processing time. 2) Schedule each job in this sequence into a machine (i.e., a frequency channel) in which the current flowtime (total processing time of the scheduled jobs) is the shortest. The LPT heuristic tries to place the shorter jobs towards the end of a schedule, where they can be used for balancing the loads. In the worst case, the makespan of LPT schedule is no longer than 4/3 times that of the optimal schedule [START_REF] Conway | Theory of Scheduling[END_REF], i.e.,

C max (LPT) C * max ≤ 4 3 - 1 3F , (24)
where C max (LPT) is the makespan obtained by LPT.

B. Numerical Results

We firstly compare the results of the proposed algorithm against brute-force search in small networks, where F = 2 and K ∈ {4F, 5F, 6F, 7F }. In these scenarios, the NOMA-SA algorithm always finds the optimal clustering scheme that minimizes the sum of processing time that is evaluated in SA loop 1. When comparing the makespan, Fig. 5 indicates that the NOMA-SA algorithm succeeds in finding solutions that perform very closely to global optima, and OMA-SA can increase the makespan by 22% as compared to NOMA-SA.

The makespan versus K for F = 4, 8, and 12 are depicted in Fig. 6, Fig. 7, and Fig. 8, respectively, where the packet size of each user varies randomly from 10 4 bits to 10 8 bits. Note that "RAND" in the figures refers to random user clustering, such that the SA loop 1 optimization is not performed but the SA loop 2 for scheduling optimization is performed. As shown by the figures, the NOMA scheme optimized by our SA-based serial collaborative method outperforms the conventional OMA scheme in terms of the makespan. Specifically, in Fig. 8, the proposed NOMA approach reduces the makespan by about 10% as compared to OMA-SA when K = 36. Moreover, given a fixed amount of frequency channels, the system performance gain of NOMA over OMA in terms of the maximum completion time increases significantly with the growing number of users. This can be demonstrated again by Fig. 8, where the benefit of NOMA against OMA-SA increases significantly from 10% to about 32% when K = 90. A possible reason is that with the increase of concurrent users, the channel profiles are more diverse, which allows more user grouping choices in NOMA that result in more efficient utilization of the limited spectral resources. The performance of NOMA with random user pairing is even worse than conventional OMA scheme. This result reveals the importance of proper user grouping on the system performance of NOMAenabled network. In addition, the performance gap between SA algorithm and LPT rule for OMA scheduling is about 2%, showing that the main performance gain of NOMA against OMA comes from the user clustering.

Fig. 9 shows the reduction in makespan by applying the proposed NOMA SA approach as compared to OMA-SA, where F = 4 and the packet size of each user varies within different ranges. We see that the benefit of NOMA against OMA increases as the number of users in the system increases. Moreover, the benefit of NOMA is more significant when the packet size variation is larger. Specifically, the reduction in makespan due to NOMA is about 16% when K/F = 16 and the packet size varies from 10 4 bits to 4 × 10 4 bits, but about 24% reduction in makespan is achieved when the packet size varies from 10 4 bits to 10 8 bits for the same K/F .

The number of iterations before the termination of Algorithm 1 for both SA loop 1 and SA loop 2 in Fig. 9 are recorded as follows. Specifically, when the packet size varies from 10 4 bits to 10 8 bits, SA loop 1 converges in a maximum of 15200 iterations, and SA loop 2 converges in a maximum of 17400 iterations. These two maximum iteration numbers increase to 26400 and 33200, respectively, when the packet size variation range is changed to [10 4 , 4×10 4 ]. Fig. 10 and Fig. 11 showcase the sum of processing time as evaluated in SA loop 1 and the maximum completion time as evaluated in SA loop 2 for a random network topology, respectively, where K = 100 and F = 10. We can see that the utility functions fluctuate but then converges as the number of iterations increases.

V. CONCLUSION

In this paper, we have proposed a SA-based serial collaborative optimization algorithm to minimize the maximum completion time of serving multiple uplink users that transmit packets of different lengths in a multi-channel NOMA system. We formulated the problem as a combinatorial optimization problem, analyzed its complexity and decomposed it into two separate parts. Simulation results have shown the convergence of the proposed algorithm. Moreover, the proposed algorithm is able to find clustering and scheduling solutions for NOMA that, on average, reduce the maximum completion time for more than 30% as compared to the OMA counterpart. APPENDIX A PROOF OF [START_REF] Mishra | Downlink connection density maximization for NB-IoT networks using NOMA with perfect and partial CSI[END_REF] As shown in [START_REF] Jang | Throughput maximization for uplink IoT NOMA systems with efficient resource allocation[END_REF], the expression of u i,j has two forms. When u i,j takes the form of case 1,

τ i + τ j -u i,j = τ i = L i R OMA i ≥ 0. (25)
When u i,j takes the form of case 2,

τ i + τ j -u i,j = τ i - L i R OMA i + R SIC,j i τ j R OMA i = R SIC,j i τ j R OMA i ≥ 0 (

Fig. 2 .

 2 Fig. 2. Illustration of the processing time of two users in a NOMA cluster.

Fig. 3 .

 3 Fig. 3. Multi-channel scheduling for NOMA clusters.

7 :η

 7 ← 0 reset counter for new state acceptance 8: Randomly generate x ′ in the neighborhood of x 9: if U (x ′ ) < U (x) then accept it 10: else 11: accept it with probability exp -U (x ′ )-U(x) T 12:

Fig. 4 .

 4 Fig. 4. Structure of SA-based serial collaborative optimization.

Fig. 5 .

 5 Fig. 5. Makespan versus K/F ratio (F = 2).

Fig. 6 .

 6 Fig. 6. Makespan versus the number of users (F = 4).

Fig. 7 .

 7 Fig. 7. Makespan versus the number of users (F = 8).

Fig. 8 .

 8 Fig. 8. Makespan versus the number of users (F = 12).

Fig. 9 .

 9 Fig. 9. Reduction on C max achieved by NOMA using SA as compared to OMA-SA versus K/F ratio under different allowable package length variations, where F = 4.

Fig. 10 .

 10 Fig. 10. Convergence result in SA loop 1.

Fig. 11 .

 11 Fig. 11. Convergence result in SA loop 2.

  26) because R SIC,j i , τ j , and R OMA i are all non-negative. The proof completes since (25) and (26) prove (8). APPENDIX B A LOWER BOUND ON THE STIRLING NUMBER OF THE SECOND KIND The number of ways of partitioning a set of n elements into k non-empty sets. Denoted by S(n, k), Stirling numbers of the second kind obey the following recurrence relation:S(n, k) = S(n -1, k -1) + kS(n -1, k). (27)We can obtain a lower bound for S(n, k) as follows:S(n, k) = S(n -1, k -1) + kS(n -1, k) > kS(n -1, k) = k(S(n -2, k -1) + kS(n -2, k)) > k 2 S(n -2, k) (28) > k n-k S(k, k) (29) = k n-k ,(30)where (28) follows because S(n -2, k -1) ≥ 0, and (29) is obtained by repeating the argument that leads to (28).

Table I :

 I Parameter settings for network topology generation User-to-BS distance 100m ≤ d k ≤ 1000m Max device transmit power P max = 23 dBm (0.2W ) , 10 8 ], [10 4 , 10 5 ], or [10 4 , 4 × 10 4 ]

	Noise spectral density	N 0 = -174 dB/Hz
	Path loss exponent	α = 3
	Bandwidth per channel	W = 180 kHz
	Data packet size (in bits)	[10 4

Table II :

 II Parameter settings for SA

	SA Parameters SA Loop 1 SA Loop 2
	T 0	150	100
	I	200	200
	minMove	5%	5%
	β	0.95	0.95

This work has been supported in part by Xi'an Jiaotong-Liverpool University Research Development Fund RDF-21-01-045.

Fixing the NOMA user clustering, the remaining scheduling problem is still NP-hard even with the assumption in[START_REF] Islam | Powerdomain non-orthogonal multiple access (NOMA) in 5G systems: Potentials and challenges[END_REF]; see Section II-B.

This can happen if u i,j > τ i , u i,j > τ j , and when scheduling users i and j using OMA, both user i and user j can finish earlier than the case when scheduling the NOMA cluster (i, j).


Xuan Wu
email: xuan.wu@ensam.eu

Virginie Goepp
email: virginie.goepp@insa-strasbourg.fr

Ali Siadat
email: ali.siadat@ensam.eu

-Morocco "New advances and challenges for sustainable and smart industries" FUNCTION ANALYSIS OF MANUFACTURING EXECUTION SYSTEMS (MES) IN CYBER PHYSICAL PRODUCTION SYSTEMS (CPPS) 1,2

Keywords: Cyber physical systems, Cyber physical production systems, Manufacturing execution systems, Functional modelling (IDEF0), 5C architecture, ANSI/ISA-95

Cyber Physical Systems (CPSs) play a crucial role in the Industry 4.0 paradigm. The application of CPSs in production environments gave rise to the term Cyber Physical Production Systems (CPPSs). The emergence of CPPSs transforms the automation pyramid into a decentralized structure. Manufacturing Execution Systems (MES) are being highly influenced by this transformation. Currently, there remains a dispute on the role of MES in CPPSs, specifically with respect to the functions they support in CPPSs. This paper aims to shed light on the role of MES in CPPSs by analyzing the functional models of the two. Firstly, functional models of CPPSs and MES are proposed in the form of IDEF0 diagrams. Then, data flows between the functional models of CPPSs and MES are described and the functions of the two are compared. We find that MES will continue to play a significant role in CPPSs, but it will take new forms. The next-generation of MES requires greater flexibility, dynamism and improved functionality.

INTRODUCTION

In the era of Industry 4.0, Cyber Physical Systems (CPSs) play an essential role because they make the fusion of the physical and the virtual worlds [START_REF] Kagermann | Securing the future of German manufacturing industry: Recommendations for implementing the strategic initiative industrie 4.0[END_REF]. The application of CPSs in production environments leads to the development of Cyber Physical Production Systems (CPPSs). The first detailed description of a CPPS was given by [START_REF] Monostori | Cyber-physical systems in manufacturing[END_REF], which has been widely cited and broadly accepted in recent years. As a complement to Monostori's description, Cardin (2019) added several missing points, including knowledge management, decision making and adaptability, and adapted it as following (the complementary part of Cardin's description has been highlighted in bold): "CPPSs are systems of systems of autonomous and cooperative elements connecting with each other in situation dependent ways, on and across all levels of production, from processes through machines up to production and logistics networks, enhancing decisionmaking processes in real-time, response to unforeseen conditions and evolution along time". [START_REF] Monostori | Cyber-physical systems in manufacturing[END_REF] claimed that CPPSs partly break with the traditional automation pyramid. Traditionally, the Manufacturing Execution Systems (MES) was positioned as a middle layer to bridge the shop-floor with the Enterprise Resource Planning (ERP) system (left side of Figure 1). Today, the automation pyramid will be transformed into a more decentralized structure (right side of Figure 1). However, this transformation leads to vague on the role of MES in CPPSs: Does this decentralized structure of CPPSs make MES superfluous? If not, what is the role of MES in CPPSs? And how does MES need to change to fulfill this role? Some scholars attempted to answer these questions. For example, [START_REF] Rossit | Industry 4.0: Smart Scheduling[END_REF] claimed that current MES, which take care of scheduling and dispatching work orders, will be absorbed by CPPSs. This minimizes response times, increase the flexibility of the entire system, which is useful for decision making. However, [START_REF] Arica | Status and future of manufacturing execution systems[END_REF] and [START_REF] Mantravadi | An Overview of Next-generation Manufacturing Execution Systems: How important is MES for Industry 4.0?[END_REF] outlined that with data collection, analysis, and communication functions across the value chains, MES should play a central role in Industry 4.0 manufacturing systems. ZVEI, (2017) working group also indicated that the significance of MES in Industry 4.0 will increase and not decrease. But the focus will shift away from simple execution management towards comprehensive coverage of all Manufacturing Operation Management (MOM) activities and away from the control of production towards optimization including the incorporation of events. Almada-Lobo, (2016) pointed out that MES is required to change to cope with the challenges created by Industry 4.0, especially considering the following four aspects: decentralization, vertical integration, connectivity and cloud computing and advanced analysis. We find that, in current academic research, there still remains a dispute on the role of MES in CPPSs, ranging from "no role because the hierarchical structure will be dissolved" and "fewer and fewer role because MES functions will be cannibalized by enhanced CPPSs and ERP functions" right up to "central role because MES will be the decision-making MOSIM '20 -November 12-14, 2020 -Agadir -Morocco center of an organization and create the optimal value chain". The existing studies are either based on literature reviews or empirical methods, but no one builds functional models to elaborate and compare the functions between CPPSs and MES. Therefore, this paper aims to shed light on the role of MES in CPPSs, by proposing functional models of CPPSs and MES in the form of IDEF0 diagrams and thus comparing their functions. The rest of the paper is organized as follows. In Section 2, functional models of CPPSs and MES are proposed using IDEF0 diagrams. In Section 3, the data flows between the functional models of CPPSs and MES are described and then the functions of the two are compared. Section 4 concludes this paper and indicates future work.

FUNCTIONAL MODELS OF CPPS AND MES

IDEF0 (Integration DEFinition of function modelling) offers a functional modelling methodology for activity and information-flow analysis, focusing on the hierarchical decomposition of activities and the interaction flows between activities (IDEF0, 1993). In this section, using IDEF0 diagrams, the functional model of a CPPS based on the 5C architecture proposed in [START_REF] Lee | A cyberphysical systems architecture for Industry 4.0-based manufacturing systems[END_REF] and the functional model of an MES based on the ISA-95 standard (ISA-95 2005) are worked out. Level IV (cognition level) generates a thorough insight of the whole system, which supports expert users to make decisions. At Level V (configuration level), CPPSs can self-configure and self-adapt to return to normal behaviors. For each level, we can extract a major activity, which is "acquire raw data from the physical space", "convert raw data to meaningful data", "build the cyber space", "make decisions", and "self-configure and self-adapt", respectively.

Functional models of CPPSs

5C architecture

Functional model of CPPSs

According to the description of the 5C architecture, each level has a major activity. Therefore, the context diagram A0 "Implement a CPPS" can be decomposed into these five sub-activities, A1-A5, as shown in Figure 3, which shows the activities involved in the manufacturing process along with the information flows required to support these activities. A detailed description of the function modules and their relationship is given as follows.

A1-Acquire raw data from the physical space

This box describes the main activity in the Level I of the 5C architecture. The input is "operational commands", which is transformed into "raw data". The raw data refers to measurements data from plant instrumentation. The activity uses the mechanisms "data sources (such as sensors and controllers)", "communication protocols (such as MTConnect)" and "data acquisition methods", with the controls "production rules", "production capability" and "production routing".

A2-Convert raw data to meaningful data

This box describes the main activity in the Level II of the 5C architecture. The input that is linked to the prior output "raw data" of the activity A1 is transformed into "meaningful data". The meaningful data refers to the data at a higher semantic level, such as movement events of resources and statuses data of resources. Data analysis technologies as the mechanism are needed for getting meaningful data. In addition to having the same controls as the activity A1, another control is the data format.

MOSIM'20 -November 12-14, 2020 -Agadir -Morocco A3-Build the cyber space This box describes the main activity in the Level III of the 5C architecture. It has one input which is linked to the prior output "meaningful data" of the activity A2, one output "simulation data of status and behaviors of CPPSs", and the same controls as the activity A1. Simulation software and information systems are needed to support the activity and its function.

A4-Make decisions

This box describes the main activity in the Level IV of the 5C architecture. It has one input which is linked to the prior output "simulation data" of the activity A3, and one output "corrective and preventive decisions". Simulation software, information systems, proper presentation tools and expert users are needed to support the activity and this function. In addition to having the same controls as the activity A1, another control, the risk criterion, is needed to support the correct decisions to be taken.

A5-Self-configure and self-adapt

This box describes the main activity in the Level V of the 5C architecture. The input that is linked to the prior output "decisions" of the activity A5 is transformed into "adjusted production activities". It has the same controls as the activity A1. Physical entities in shop floor are needed to support the activity and its function.

Functional models of MES

ANSI/ISA-95 standard

The part 3 of the standard ISA-95 (2005) gives a detailed overview of all the activities and data flows of MOM in the area of production, inventory, maintenance and quality operational management. MOM could be regarded as the upper limited range of functions involved in MES applications. Generally speaking, the main focus of current MES applications is only on the production operations management. However, future MES functions will have to be developed towards MOM. Therefore, in this paper, MES refers to MOM, which covers all areas of operations management. One of the main advantages of ISA95 is the definition of the main functions and data flows within MOM. Therefore, for this study, MES function's classification is based on the ISA-95.

Functional model of MES

The ISA-95 described the generic activities in MOM, as shown in Figure 4. Each of the four areas of operations management (production, inventory, maintenance and quality) is composed of these eight activities. According to this, the context diagram B0 "Implement an MES" was decomposed into eight sub-activities, B1-B8, as depicted in Figure 5. Activity B1 (Manage definitions of work) manages the work information and develop work rules. Activity B2 (Manage resource) manages the information about the resources including machines, tools, labor, materials and MOSIM'20 -November 12-14, 2020 -Agadir -Morocco energy. Activity B3 (Detailed schedule) makes the schedule and determines the optimal use of resources to meet the schedule requirements. Activity B4 (Dispatch work) dispatches work lists to equipment and personnel. Activity B5 (Manage executions of work) directs the performance of work, as specified by the dispatch list. Activity B6 (Track work) summarizes and reports information about the personnel and equipment used, the material consumed, material produced and other relevant data such as costs and performance analysis results. Activity B7 (Collect data) gathers, compiles and manages data for specific work processes. Activity B8 (Analyze performance) analyzes and reports performance information of production systems.

The data flows between these activities are extracted from part 3 of the ISA-95. For example, one of the tasks in product definition management, a subset of "B1manage definitions of work", is extracted from the ISA-95: providing product production rules to personnel or other activities. From this task, we can conclude that the production rule is the output of B1. Like this, the data flows between these activities can be found in the ISA-95, so we will not discuss them here.

FUNCTION ANALYSIS OF MES IN CPPS

Based on the established functional model of CPPSs and MES, this section describes the data flows between CPPSs and MES, which represents the functional interfaces between them. Then based on the data flow, this section analyses which functions of MES are redundant in CPPSs and which functions are central.

There are two-way data flows in MES. On the one hand, MES deals with the top-down data flow: the requirements provided by the organizational level must be transformed into detailed operational commands to the production areas. On the other hand, it manages the bottom-up data flow. MES collects shop floor data, analyzes it, and extracts useful information to provide key performance indicators (KPIs) to the company for making commercial decisions and improving production performance. Therefore, the data flows between CPPSs and MES are described as follows.

• Data flow 1: "equipment and process rules" from the output of the Activity B1  mechanisms of the Activity A1 to A5.

Equipment and process rules define the specific production/inventory/maintenance/quality test instructions that are executed in the shop floor based on the specific assigned tasks, such as programs for CNC machines.

• Data flow 2: "operational commands" from the output of the Activity B5  the input of the Activity A1.

Operational commands define the commands to start or complete elements of a work order in the shop floor, such as procedures for setting up machines or maintenance of machines.

• Data flow 3: "operational responses" from the output of Activity A3  the input of the Activity B5.

Operational responses define information received from the shop floor in response to commands, such as the completion or status of elements of work orders.

• Data flow 4: "equipment and process specific data, resource data, maintenance history, quality test data, inventory data " from the output of Activity A3  the input of the Activity B7.

Equipment and process specific data defines the data about the process being performed and the resources involved.

• Data flow 5: "Production, maintenance, quality and inventory related KPIs" from the output of Activity B8

 the input of Activity A4

The manufacturing operations management areas have different sets of performance indicators, which are used together to make decisions and to monitor the realization of enterprise business objectives.

• Data flow 6: "current resource availability" from the output of Activity A3  the input of Activity B2.

Resources include machines, tools, labor (with specific skill sets), materials, and energy.

• Data flow 7: "resource capability" from the output of Activity B2  the input of A3

Resource capability defines the committed available, or unattainable of resources. The information is based on the current statuses and future needs (as identified in Activity B3).

• Data flow 8: "production and movement events" from the output of Activity A3  the input of the Activity B6.

Production and movement events include the movements of the path of materials and updates to lot locations.

For the data flow 1, 3, 4, 6, 7 and 8, MES mainly serves as an information platform that contains all information in the production, inventory, maintenance and quality test operations management and we find that the functions of MES are very similar to those of digital twin in terms of work definitions management, resource management, work track and data collection. For example, in terms of the MES's function "data collection", digital twins store all data from the field and other information systems in real-time, as in the MES. Digital twin refers to the virtualization of the physical resources in the cyber world and its main purpose is to facilitate the decision-making process through real-time simulation [START_REF] Tao | Digital twin-driven product design, manufacturing and service with big data[END_REF]. Therefore, future MES can be used as a digital twin due to its ability to provide digital images and monitor manufacturing processes in real-time.

For the data flow 5, MES analyzes and reports performance information to expert users for decision-making.

In the Industry 4.0 context, the function of performance analysis in MES becomes more complex and important. With the development of information and communication technology (ICT), the amount of data collected from the shop floor increases tremendously. MES should always be required to present the most relevant data to the humans for decision making.

For the data flow 2, MES sends the operational commands to the physical space of CPPSs. The function of executions management is equally important because it is not only necessary to have data coming from the shop floor to MES, but also to be able to act from the MES to the shop floor, by offering operational commands on the production control system. In the CPPSs context, it should be possible in the future to decentralize operational decisions and delegate a part of decisions made previously in MES to lower levels of production systems, such as smart machines and smart products. This reduces the number of tasks of this function in terms of the operational level decisions. Since MES covers all the areas of MOM, it can collect all data related to production, quality, inventory and maintenance. This endows MES with the holistic insights for the manufacturing processes and, which can make commands before any disruptions.

CONCLUSION AND FUTURE WORK

The paper studies the role of MES in CPPSs by analyzing the functional models of CPPSs and MES. MES's functions are still required for CPPSs, but they need to MOSIM'20 -November 12-14, 2020 -Agadir -Morocco be presented in new and different forms, with some being incorporated into the CPPSs' functions, some expanding their horizons and some decentralizing to the local intelligent entities. For example, the MES's function "detailed schedule" will only have high level tasks, making rough planning of manufacturing orders, while flexible planning of partial orders will be made by local entities of CPPSs through taking into account local resource availabilities in real-time. Correspondingly, MES software will move from monolithic software systems towards functional modules that will be implemented in a highly distributed way. Therefore, the next-generation MES requires greater flexibility, dynamism and improved functionality than today's MES.

In the future work, validation in a real platform will be worked out to evaluate the quality of the proposed model. Besides, we will implement a "dynamic link" between CPPSs and MES from a functional perspective. In our view, the functional level is the best suited for this implementation as it enables to free of the underlying technical solutions. Indeed, the MES and CPPS interrelated functions cannot be engineered independently and then simply combined together. The idea is rather to coordinate the design of both systems to explore hidden functions and ensure that their interactions are mutually supportive. This "dynamic link" is twofold: one the one hand, if the configuration of CPPSs (e.g. the layout of CPPSs) changes, the organization and business processes will be changed automatically. One the other hand, if new functionalities of EISs are needed, CPPSs can configure themselves automatically to meet these new requirements. The digital twin could be a solution to implement such dynamic link as the change of state of the physical object will have impacts on the state of the cyber object automatically and vice versa in digital twins. [START_REF] Cimino | Review of digital twin applications in manufacturing[END_REF] pointed out the future possibility of integrating the digital twin of the field devices with the MES in order that decision making could be made in a fully reactive way through the digital twin to the MES.

Figure 1 :

 1 Figure 1: Decomposition of the automation hierarchy with distributed services (VDI/ VDE, 2013)

Figure 2 :

 2 Figure 2: 5C architecture for the implementation of CPSs (Lee et al., 2015)

Figure 3 :

 3 Figure 3: Activity A0 -Implement a CPPS.

Figure 4 :

 4 Figure 4: Generic activity model of manufacturing operations management (ISA-95, 2005)

ACKNOWLEDGMENTS

The paper of the first author is financially supported by China Scholarship Council under Grant 201706020154.