{'fraction_non_alphanumeric': 0.04400462962962963, 'fraction_numerical': 0.009328703703703704, 'mean_word_length': 4.577921239509361, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
IRS Assisted NOMA Aided Mobile Edge Computing with Queue Stability: Heterogeneous Multi-Agent Reinforcement Learning

By employing powerful edge servers for data processing, mobile edge computing (MEC) has been recognized as a promising technology to support emerging computation-intensive applications. Besides, non-orthogonal multiple access (NOMA)-aided MEC system can further enhance the spectral-efficiency with massive tasks offloading. However, with more dynamic devices brought online and the uncontrollable stochastic channel environment, it is even desirable to deploy appealing technique, i.e., intelligent reflecting surfaces (IRS), in the MEC system to flexibly tune the communication environment and improve the system energy efficiency. In this paper, we investigate the joint offloading, communication and computation resource allocation for IRS-assisted NOMA MEC system. We firstly formulate a mixed integer energy efficiency maximization problem with system queue stability constraint. We then propose the Lyapunov-function-based Mixed Integer Deep Deterministic Policy Gradient (LMIDDPG) algorithm which is based on the centralized reinforcement learning (RL) framework. To be specific, we design the mixed integer action space mapping which contains both continuous mapping and integer mapping. Moreover, the award function is defined as the upper-bound of the Lyapunov drift-plus-penalty function. To enable end devices (EDs) to choose actions independently at the execution stage, we further propose the Heterogeneous Multi-agent LMIDDPG (HMA-LMIDDPG) algorithm based on distributed RL framework with homogeneous EDs and heterogeneous base station (BS) as heterogeneous multi-agent. Numerical results show that our proposed algorithms can achieve superior energy efficiency performance to the benchmark algorithms while maintaining the queue stability. Specially, the distributed structure HMA-LMIDDPG can acquire more energy efficiency gain than centralized structure LMIDDPG.

To alleviate the heavy communication overheads caused by centralized scheme, distributed resource allocation frameworks have been investigated. The most typical of which is multi-agent reinforcement learning (MARL) [12]. The MARL has the following advantages: firstly, agents learn their strategies in the distributed manner by observing local environment. Secondly, agents can be designed to share experience with each other based on customized purposes. Thirdly, multi-agent is more robust when some agents fail, as the remaining agents can take over the tasks.
In [10], the computation offloading mechanism with resource allocation has been formulated as a stochastic game. The independent learners based multi-agent Q-learning algorithm has been proposed. In addition, the authors designed a distributed multi-agent deep reinforcement learning scheme to minimize the overall energy consumption of in the small cell networks [11]. To further decrease the computation complexity and communication overheads of the training process, a federated DRL scheme which only share the model parameters has also been proposed in [11].
Energy-efficient and spectral-efficient task offloading is the key to the success of MEC, which drives the deployment of NOMA in MEC offloading [13]. Explicitly, NOMA allows multiple wireless devices to occupy the same frequency and time slot resources to transmit data to the BS, and BS uses SIC technology to decode the received information. Based on the system optimization metrics, the existing NOMA-aided MEC works are generally divided into two categories: energy consumption [14], [15] and time delay [16], [17]. In [14], a hybrid NOMA-MEC offloading strategy, which combines conventional orthogonal multiple access and pure NOMA was proposed. Moreover, the multi-objective optimization problem has been formulated to minimize the energy consumption. Similarly, the authors also proposed an energy consumption optimization problem by optimizing the resource allocation and sub-channel assignment with the latency constraint in NOMA enabled MEC in [15]. To minimize the time delay in the NOMA-MEC system, the authors in [16] discussed that under the same power constraint, the deployment of NOMA can achieve substantial time delay reduction compare to time division multiple access (TDMA). To reduce the user average offloading delay, the optimal power allocation is obtained by convex programming, followed by iteratively semidefinite relaxation and convex-concave for NOMA user paring and offloading decision making. To minimize the maximal offloading latency of NOMA users, the authors exploited the successive convex approximation method to solve the complicated non-convex problem [17]. Apart from aforementioned research, there are several works investigated the RL framework for NOMA-MEC [18], [19]. An online algorithm based on DRL framework to learn the near-optimal solution by minimizing the total energy consumption of end devices in multitask NOMA-MEC system was proposed in [18]. Additionally, to minimize the long-term network computation cost, in [19], the authors proposed a cooperative multi-agent deep reinforcement learning framework to learn the decentralized policies at end users for task offloading decisions and local execution power allocation.
Recently, as a revolutionizing technology, IRS has been proved to be able to improve the performance of wireless communications with new degrees of freedom through highly controllable intelligent signal reflection [20]. A few works exploited the beneficial role of IRSs in MEC, and they can be divided into two categories based on different optimization metrics: energy consumption [21] and system throughput [22], [23]. Authors in [21] designed the resource allocation problem in the wireless powered IRS-MEC system by utilizing the alternating optimization technique and successive convex approximation method to minimize the system energy consumption. Differently, [22] focused on optimizing IRS-MEC system's sum computational bits and the Lagrange dual method and Karush-Kuhn-Tucker conditions were proposed. What's more, [23] also focused on maximizing the total computation bits maximization problem for IRS-enhanced wireless powered MEC networks and the iterative algorithm was proposed to solve the non-convex non-linear optimization problem.

B. Motivation and Contributions
As discussed above, with the development of wireless communications, the research in MEC has evolved from fundamental communication environment [7]- [11] to NOMA-aided [14]- [17] and IRS-assisted systems [21]- [23]. Despite the advantage of NOMA-aided MEC system with resource-efficient transmission, there are several limitations, such as the fluctuation of wireless channels, path loss, high user mobility, and so on that can influence the performance of the offloading communication. Although increasing the number of antennas at both transmitter and receiver can overcome such problems, the higher energy consumption occurs at the same time.
As an near-zero energy consumption technique, IRS has been identified to support NOMA communication with enhanced performance from four aspects: tuned channel gains, improved fair power allocation, enhanced coverage range, and high energy efficiency [5], which motivates us to explore the IRS-assisted NOMA-aided MEC system. For optimization metrics design, minimizing energy consumption [14], [15] or time delay [16], [17] are two widely explored problems in NOMA-enable MEC. Apart from saving energy consumption [21], system computation throughput [22], [23] has been considered for IRS-assisted MEC. However, it is necessary to design a customized optimization metric for IRS-assisted NOMA-aided MEC system that can include both energy consumption, system throughput and time latency at the same time.
Although the conventional optimization algorithms (i.e., Block coordinate descent (BCD), successive convex approximation, and so on) have been broadly proposed to solve different multiconvex or non-convex resource allocation problems in MEC, due to the computational difficulty and the requirements of the statistical knowledge of the MEC systems, it is impractical to deploy them in the highly dynamic environment. Variant RL frameworks become one of the efficient alternatives to solve the complex joint resource allocation problems without the prior statistical knowledge [7]- [11], [19], [24]. Normally, the resource allocation tasks in MEC are mixed integer problems. However, most of the existing RL-based algorithms are limited to the action space with pure discrete [7]- [11] or pure continuous [19] action in MEC systems, which motivates us to consider the mixed integer mapping in the algorithm design. Although [24] considered the multi-stage stochastic mixed integer non-linear programming problem and applied the centralized structure of RL framework to solve the decoupled problems, it is inevitable to consider the distributed manner MARL framework which allows end users to make decisions without outside information and can further save the heavy communication overheads. Moreover, the existing MARL frameworks on solving the MEC allocation problems mainly have homogeneous agents which have similar state and action spaces. However, there is no work considering the agents with different state and action spaces.
Motivated by the aforementioned literature review, we consider the IRS-assisted NOMA-aided MEC system and propose the energy efficiency with queue stability as the evaluation metric for optimal resource allocation strategy design. We further design both centralized RL framework and distributed MARL framework to solve the problem with mixed integer action space. By considering the different characteristics of the devices in MEC, the heterogeneous agents are designed in MARL framework. The main contributions are summarized as follows.
• We formulate the energy efficiency maximization problem with queue stability in IRSassisted NOMA-aided MEC system by optimizing the mixed integer action space, which consists of edge computing decision, local computing decision, transmit power allocation, local computing frequency, edge computing frequency and IRS phase shifts.
• We further derive the upper-bound of the Lyapunov drift-plus-penalty function based on our proposed problem, which can reflect the system throughput, system power consumption and the queue stability of the designed problem.  • To save the communication interaction overheads and preserve privacy at the execution stage, the distributed RL framework algorithm named as Heterogeneous Multi-agent LMIDDPG (HMA-LMIDDPG) is proposed. Specially, there are heterogeneous multiple agents which include homogeneous end devices (EDs) and heterogeneous BS.
• The proposed LMIDDPG and HMA-LMIDDPG can maintain the queue stability and outperform benchmark algorithms with higher energy efficiency. Moreover, the distributed framework HMA-LMIDDPG can achieve further performance gain compared to the centralized framework LMIDDPG.
The rest of this paper is organized as follows. The IRS-assisted NOMA-aided MEC system model is introduced in Section II. In Section III, the task offloading and resource allocation problem is formulated. The proposed algorithm LMIDDPG and HMA-LMIDDPG are presented in Section IV. Numerical results and conclusions are then provided in Section V and Section VI, separately.

II. SYSTEM MODEL
As illustrated in Fig. 1, we consider the IRS-assisted NOMA-aided MEC system, where a set of single antenna devices, denoted by U = {u 1 , ..., u k , ..., u K } with k ∈ K = {1, 2, ..., K} are served by one BS and one IRS. The EDs can be different kinds of sensors and mobile devices that are energy-constrained and computation-limited. The system is assumed to operate at a time slotted structure which is indexed by an integer t ∈ T = {1, 2, ..., T }, and the slot length is ∆T . At each time slot t, the computation task is assumed to be generated by ED u k as the independent and identically distributed sequence of Bernoulli random variables with common parameter ζ k (t) ∈ [0, 1], and the task generation indicator G k (t) = {0, 1} (i.e., ζ k (t) = P(G k (t) = 1) = 1 − P(G k (t) = 0)). Specifically, G k (t) = 0 indicates that there is no task generated by ED u k at time slot t. Moreover, G k (t) = 1 represents that a task with constant size A k (t) in bits is generated. Since the BS which has much higher computation ability, the intensive tasks generated at EDs can either be processed locally or offloaded to the BS through NOMA. In this work, we consider the partial offloading 1 , which means each ED makes its own decision on whether to offload the task (i.e., O

A. Communication Model
The overall bandwidth B is equally divided into N channels, each with the bandwidth as In this work, we focus on exploring the sub-channel NOMA-based IRS-aided system. To be specific, we assume that there are K users occupying the same sub-channel. The movement of EDs are assumed to follow Gauss-Markov mobility model [25]. To be specific, the velocity of u k at time slot t can be modeled as v is the velocity vector, W k (t) = [w x k (t), w y k (t)] ∼ N (0, ς 2 ) denotes the uncorrelated random Gaussian process, and α = [α x , α y ] represents the memory level, v = [v x ,v y ] denotes the asymptotic mean, andκ = [κ x ,κ y ] represents the asymptotic standard deviation of the velocity. Thus, given the velocity of v k (t), the position P k (t) = [x k (t), y k (t)] at time slot t of ED u k can be updated as The IRS in the system has M passive reflecting elements with the reflection coefficients matrix as Θ = diag{Θ 1 , Θ 2 , ..., Θ M } = diag{λ 1 e jθ 1 , λ 1 e jθ 2 , ..., λ 1 e jθ M }, where λ m ∈ [0, 1] and θ m ∈ [0, 2π] are amplitude and phase shift of the m th reflection element on IRS [26]. The equivalent baseband time-domain channel consists of three parts in the system: ED-IRS, ED-BS, and IRS-BS. To be specific, the channel gain between ED u k and the IRS is g k ∈ C 1×M , the channel gain between ED u k and the BS is h d k , and the channel gain between IRS and BS is h r ∈ C 1×M . Moreover, in practice, the IRS is installed near the EDs and may be located far away from BS.
Hence, the ED-IRS channels would be LoS dominant that can be denoted as g k = β (DI) kg k with β (DI) k as the large-scale path loss andg k as the deterministic LoS channel between ED u k and the IRS. Differently, EDs may be far away from the BS with rich scatterers. Thus, we assume that the ED-BS channels are Rayleigh fading, which can be written as h d k = β (DB) kh d k with large-scale fading coefficient β (DB) k and small-scale fading channelh d k . Furthermore, the IRS-BS channel can be modeled as Rician fading, which can be expressed as with β (IB) as the distance-dependent path loss factor between IRS and BS, δ as the Rician factor 2 which represents the ratio between LoS componenth r,(LoS) and non-LoS (NLoS) component Conventionally, for NOMA systems, SIC decoding order is very crucial which is determined by the channel gains. In IRS-aided system, the IRS reflection coefficients can affect the combined channel gains. When the signal from ED u k needs to be decoded at the receiver BS, the achievable capacity in our proposed system for ED u k , when prior decoding order were decoded, can be formulated as where p k is the power allocated to ED u k , h k = h d k + h r Θg k is the offloading channel which consists of both ED-BS and cascaded ED-IRS-BS paths, I k is the interference of the signals from EDs whose decoding order is latter than k, and the noise is assumed to be white Gaussian channel noise with variance as σ 2 . Specifically, I k = {q∈K|η(q)>η(k)} |h q | 2 p q , which represents the overall interference of the rest EDs whose decoding order is latter than ED u k in this subchannel (i.e., η(q) and η(k) are the decoding order mapping 3 for ED u q and ED u k , respectively.
Thus, η(q) > η(k) denotes that ED u q has the decoding order latter than ED u k ). Therefore, the overall system communication throughput can be formulated as R = K k=1 R k .

B. Task Execution Model
Since the main purpose of MEC is that each ED can execute its intensive tasks locally or offload to the edge server who has powerful computation ability. Hence, this subsection formulates the local computation model and edge computing model, respectively.

1) Local Computing: When O
(l) k (t) = 1, ED u k executes task locally at time slot t, the computing energy consumption at each time slot t can be formulated as E where π (l) is the effective capacitance coefficient of processor chip [27]. Moreover, the computing rate can be given as R indicates the CPU frequency on ED u k , and c k is the CPU cycles required to process a single bit of data in cycles/bit for ED u k . Therefore, the whole system energy consumption for local computing can be formulated as represents the EDs set who decide to execute the task locally. The system overall local computing rate can be denoted as k (t) = 1, ED u k decides to offload the task for edge computing at time slot t. Compared to local EDs, edge server which is BS in this paper, has much powerful computation ability. The computing energy consumption at each time slot t at BS can be formulated as E (e) (t) = π (e) f (e) (t) 3 , where π (e) is the effective capacitance coefficient of processor chip on BS. Moreover, the computing rate at BS can be given as where f (e) (t) indicates the CPU frequency allocated to process the offloaded tasks on BS, c 0 is the CPU cycles required to process a bit of data in cycles/bit on BS.

C. Task Offloading Communication Model
Since different tasks generated by u k can be assigned to the same sub-channel with different transmit power, the achievable capacity in (2) should be further presented as with time index t, Q (l) k (t) represents the task queue size of ED u k at time slot t, and p k (t) represents the transmit power of ED u k at time slot t. The system overall arrival workloads at BS within a time slot can be written as where K (e) (t) denotes the subset of EDs who decide to offload the task at time slot t. The offloading transmission energy consumption can be given as E (o) (t) = k∈K (e) (t) p k (t) + M̟, where M̟ is the overall power consumption of IRS with ̟ denoting the power consumption of each passive element on IRS [28]. During the offloading process, we assume that each ED adopts discrete transmit power control, with the transmit power values bounded by ρ max . At time slot t, each end ED u k chooses its transmit power p k (t) from the range [0, ρ max ] with different transmit power levels to determine whether to execute the task locally (i.e., p k (t) = 0) or to offload with selected power (i.e., p k (t) ∈ (0, ..., ρ max ]).

D. Task Queue Model
The task queue models at both ED u k and BS are illustrated in Fig. 2. The size of the tasks that can be processed at ED u k and BS during time slot t can be respectively formulated as and Each ED u k has the local queue Q (l) k (t) for unprocessed workloads. Similarly, the BS has a buffering offloaded queue Q (e) (t) to maintain the workloads received from EDs for edge computing. The input to Q (l) k (t + 1) for the next time slot t + 1 consists of both generated task data and the remaining data, which can be updated as with R k (t) as the amount of data transmitted to BS during a time slot interval ∆T and [·] + denotes max{·, 0}. Moreover, the task queue at the BS can be updated as which is the combination of offloading workloads and the remaining data after edge computing.

III. PROBLEM FORMULATION
In this paper, we formulate the system overall energy efficiency maximization problem in the NOMA-based IRS-aided system by jointly optimizing the task offloading decisions (O (e) , O (l) ), and transmit power of the EDs p = {p 1 , ..., p K }, the CPU frequency on EDs k }, the CPU frequency on BS f (e) , as well as the IRS reflection coefficients Θ. We , Θ} to denote the whole solution set.

A. Problem Formulation
The long-term overall energy efficiencyῩ is the sum of the system rate over the sum of the system energy consumption, which can be formulated as where is the sum of the overall system rate at time slot t, and is the overall system energy consumption at time slot t.
Remark 1. As shown in Definition 1, task queue stability is guaranteed if the length of the queue is finite. According to the Little's Theorem, the system average delay is proportional to the average queue length based on the given traffic arrival rate. Therefore, queue stability is a crucial constraint to describe the system service delay.
Hence, the optimization problem can be formulated as Constraints (10b,10c) denote the offloading decision set and local computing decision set. Constraint (10d) denotes the transmit power of offloading. Constraints (10e) and (10f) mean that the maximum CPU frequency constraints of both EDs and BS. Constraint (10g) restricts the IRS reflection coefficient. Constraint (10h) denotes that the task queues are guaranteed to be stable. It is noticed that (10a) is a long-term stochastic optimization problem in which multiple decision variables are adaptive to the dynamic system. However, due to the high stochastic and unpredictable of the dynamic system, it is impossible to acquire comprehensive knowledge for offline optimization. Moreover, (P0) is a mixed integer nonlinear programming problem , and Θ. Therefore, it is even challenging and inefficient to find the optimal solution through conventional optimization techniques which decouple the problem into sub-problems and solve them separately.

B. Problem Reformulation
It is observed that the energy efficiency optimization problem in (P0) is a stochastic nonlinear programming problem, which is also a NP-hard problem with nonlinear fractional property. For , Θ * } denotes the feasibility solutions of (P0) with the optimal solution sequence {F * (t) t∈0,1,2,... } that can maximizeῩ from a long-term perspective under all constraints, which can be written as Theorem 1. The optimal energy efficiencyῩ opt can be achieved if and only if Proof. This theorem can be proved from two aspects: the necessity and the sufficiency. Firstly, based on the necessity, for any feasible solutions F , we havē Then, the formula can be rearranged as which completes the necessity proof. Secondly, based on sufficiency proof, we firstly assume that F ′ is the optimal solution of (12). Thus, for any feasible solution, we havē Then, we can obtainR Therefore, F ′ is the optimal solution of (10a), which completes the sufficient proof.
Thus, with the support of Theorem 1, we further transform the energy efficiency maximization problem as However, the optimal value ofῩ opt cannot be known in advance. To overcome such difficulty, with t ∈ {1, 2, ..., T } and Υ (0) = 0. It can be observed that in (18), Υ (t) can be treated as the parameter which is mainly determined by the past task offloading and resource allocation decisions. Then, the optimization problem (P1) can be reformulated as Similar to the transformation in [30] on handling the stochastic optimization problems in renewal systems, it has been proved that the reformulated problem (P2) can effectively solve the original problem.

C. Lyapunov Optimization
As there is the task queue stability constraint (10h) in the optimization problem (P2), Lyapunov method can be adapted to solve the stochastic optimization problem and investigate the energy efficiency and task queue length trade-off. Therefore, with the system queue length Q(t) = Q (l) at time slot t, the Lyapunov-function and the Lyapunov drift can be written as and ∆L (Q(t)) = L (Q(t + 1)) − L (Q(t)) .
Furthermore, the Lyapunov drift-plus-penalty function at each time slot t can be formulated as where V is the positive weight control parameter to tune the trade-off between the drift and the penalty. To ensure the feasibility of the Lyapunov drift-plus-penalty optimization with finite queue capacity, the enforcement constraint has been proposed in [31]. Specifically, this paper indicates that the queue stability capacity constraint can be realized by carefully adjust the factor V value. In Lyapunov theorem, to minimize the function in eq. (22), the strategy can be refined by deriving an upper bound which can be provided as where Q The proof can be found in Appendix A. Therefore, the optimization problem in (P2) can be achieved by minimizing IV. PROPOSED ALGORITHM Deep Q-network (DQN) is an advanced RL algorithms which can handle high-dimensional discrete state spaces. Although the problem with continuous action spaces in our system can be transferred to discrete format through quantization, the dropping precision and significantly increased complexity makes it infeasible to exploit DQN to solve our problem directly. Actorcritic RL structure has been proposed to solve the problem with continuous action spaces. This is because the actor network outputs an action through the actor deep neural network (DNN) at each time step and the critic network evaluates the reward or the Q-value of a given input state.
As the critic network learns which states are better or worse, the actor uses this information to teach the agent to seek out good states and avoid bad states. There are three popular algorithms under this actor-critic structure: Proximal Policy Optimization (PPO) [32], advantage actor critic (A2C) algorithm [33] and deep deterministic policy gradient (DDPG) [34]. Among three actorcritic based algorithms, PPO is more sample expensive which is because the set of trajectories should be obtained by running the current policy over certain time slots. However, after the policy being updated, the old trajectories generated based on previous policy are no longer applicable for training the network. Thus, in this paper, considering the dynamic IRS-assisted NOMA-aided MEC work, it's more practical to consider the sample efficient algorithms. The main difference between A2C and DDPG is the way they apply critic network. In A2C, the critic is treated as a baseline from empirical trajectories, so whether to use only critic is just an option to improve the robustness of the training as the random selection prevents the back-propagation for training.
However, in DDPG, the policy is deterministic and the gradient can be acquired from the Qvalue, where Q-value is obtained from the critic network and the actions generated by the actor network. More importantly, DDPG is off-policy which means that many history trajectories can be utilized for training.
In order to adapt RL-based techniques to solve (P2), we reformulate the problem into the RL framework which consists of a few key elements, such as agent 4    r(t) = R(s(t), a(t)) and changes the state from s(t) to s(t + 1). In other words, the transition probability P(s(t + 1), r(t)|s(t), a(t)) of the MDP which is defined as the transition from state s(t) to s(t + 1) with reward r(t) when the action a(t) is taken according to the policy.
In this section, we firstly propose the centralized framework called LMIDDPG which can solve the problem (P2) with the mixed integer variables as the action space. At the same time, the negative of the upper bound Lyapunov drift-plus-penalty function has been treated as the reward function. In order to further save the communication cost at both training and execution stages, we then propose the distributed framework called HMA-LMIDDPG with heterogeneous multi-agent.

A. Centralized RL Framework: Lyapunov-Function-Based Mixed Integer Deep Deterministic Policy Gradient
In fundamental DDPG, the two different DNNs namely critic network Q(s, a|Φ Q ) which approximates the Q-function, and actor network χ(s|Φ χ ) which approximates the policy function χ, respectively. To be noticed, Φ Q and Φ χ denote the weights of the critic and actor network DNNs. The process of centralized training and execution framework is shown in Fig. 3(a). Under this framework, the BS is the agent that can observe the environment by receiving the EDs states.
With this interaction, BS will obtain the experience buffer memory at each time slot t to train the critic network Q(s, a|Φ Q ) and the actor network χ(s|Φ χ ). During the training and execution stages, EDs receive the sub-action instructions which are determined by the actor network from the BS. Before giving the details of the training process of the proposed LMIDDPG algorithm, the design of the centralized agent and mixed integer need to be explained:

1) Centralized Agent:
In the centralized RL framework, the BS acts as the agent. The details definition of the state space, action space, and the reward are defined as follows: • State space: The state s(t) at each time slot t can be defined as s(t) = {R (o) (t), Q(t), Q(t− 1), P(t)} with the information acquired from the EDs and BS itself.
• Action space: The action vector of the whole system can be formulated as a(t) = {F (t)} with mixed integer property.
• Reward: The reward function of the whole system at time slot t can be defined as r(t) = R(s(t), a(t)) = −L(t), which is the negative value of the function derivated from the upper-bound Lyapunov drift-plus-penalty function in eq. (24).    Thus, a mapping procedure is designed to recover the true mixed integer action space a(t).

2) Mixed Integer
The structure of the actor network and mapping details can be found in Fig. 3(b). After different activation function (i.e., Θ is activated through function tanh, the rest sub-actions are activated through sigmoid), the sub-actions p   Initialization: Actor (policy) network χ(s|Φ χ ) with parameters Φ χ , critic (Q-function) network Q(s, a|Φ Q ) with parameters Φ Q , target actor network χ ′ with parameters Φ χ ′ ← Φ χ , target critic network Q ′ with parameters Φ Q ′ ← Φ Q , empty replay buffer memory M. 1: for episode = 1 to E do 2: Initialize environment, i.e., initial observation state s(1).

3:
for time t = 1 to T do % Mixed integer action space mapping: 4: Observe current normalized system state s (nor) (t) and choose normalized continuous action a (con) (t) by the actor network. 5: Map the normalized continuous action a (con) (t) into true mixed integer action space a(t). Execute a(t) in the environment. % Lyapunov queue stability award training: 6: Observe the Lyapunov-function-based reward r(t) and new state s(t).

10:
Update critic by minimizing the loss:

11:
Update actor policy using the sampled policy gradient: ∇ Φχ end for 14: end for buffer memory M at BS. Then, during the training process, a mini-batch B will be sampled from the buffer memory M to compute the target value where γ ∈ [0, 1] is the discounting factor, Q ′ indicates the Q-value calculated based on the target critic network, χ ′ s (nor) (t + 1)|Φ χ ′ indicates the normalized continuous action instructed by the target actor network χ ′ given the next normalized state s (nor) (t + 1). Then, the critic network parameters can be updated based on a learning rate ̺ Q by minimizing the loss function In DDPG, the actor network is approximated by a DNN and the policy gradient of the actor can be updated by applying the chain rule to the reward with respect to the actor network parameters Φ X . The sampled policy gradient can be derived as  Then, the actor network parameter can be updated based on the learning rate ̺ χ . To improve the stability of the learning process, the weights of both target actor network and target policy network are then updated by slowly tracking the learned parameters: with ς as the update rate.

B. Distributed RL Framework: Heterogeneous Multi-agent Lyapunov Mixed Integer Deep Deterministic Policy Gradient
Apart from the privacy and security concern, there are mainly two reasons of considering the decentralized RL framework. Firstly, during the training stage, the centralized manner allows the BS to broadcast the sub-actionS to EDs for state update. However, in the distributed manner, the BS firstly receives the EDs state and action to calculate the environment award, then only broadcasts the cooperative reward for distributed training. Secondly, during the execution stage, in the centralized manner, BS still needs to broadcast sub-actions determined by the trained actor network to EDs. Differently, in the distributed manner, EDs can make action decisions based on its own distributed actor networks with their local observation.
Most existing works have focused on homogeneous multi-agent RL with agents sharing the same state and action space. However, in our system, due to the diversity of the action variables, the multi-agent, which includes EDs and BS, have a heterogeneous composition. Different categories of the agents need to leverage their unique information and rely on other agents' specializations to collaboratively find effective policies. Therefore, we propose the distributed training and execution framework called HMA-LMIDDPG as illustrated in Fig. 4(a). Under this Algorithm 2 Heterogeneous Multi-agent Lyapunov-Function-Based Mixed Integer Deep Deterministic Policy Gradient (HMA-LMIDDPG) Initialization: Agent i (i ∈ {1, ..., K + 1}. When i = K + 1, the agent is BS. Otherwise, the agent is ED u k .) Actor network χ i (s i |Φ χi ) with parameters Φ χi , Agent i critic network Q i (s i , a i |Φ Qi ) with parameters Φ Qi , target Agent i actor network χ Each Agent i has its own initial observation state s i (1). Specially, compared to homogeneous EDs as the Agent 1 to Agent K, BS as the Agent K + 1 has different dimension of state space and action space. 3: for Time t = 1 to T do 4: for Agent i = 1 to K + 1 do     for Agent i = 1 to K + 1 do 10: Observe the new state s i (t + 1), and normalize the new state as s               17: end for 18: end for 19: end for framework, both EDs and BS act as the heterogeneous agents to train their independent critic and actor networks. During the execution stage, each ED u k can take the action decision based on its own local observation without listening to the sub-action instruction from the BS. Specifically, the heterogeneous multi-agent design is specified. In the proposed distributed RL framework, there are two categories of the heterogeneous multi agents. The first agent category includes all homogeneous EDs: • State space: In this paper, when ED acts as the agent, the state s k (t) at each time slot t can be defined as s k (t) = {Q k (t), Q k (t − 1), P k (t)}. Each ED can only access its own state.
• Action space: The action vector of the ED u k can be formulated as a k (t) = {O   The second agent category includes the BS, which is heterogeneous compared to EDs: • State space: The state s K+1 (t) at each time slot t can be defined as s K+1 (t) = {R (o) (t), Q(t), Q(t − 1), P(t)}, with the information acquired from the EDs and BS itself.
• Action space: The action vector of the whole system can be formulated as a 0 (t) = {f (e) (t), Θ(t)}.
What's more, the system reward is defined as the cooperative reward that can be acquired at the BS: In our proposed distributed RL framework, during the training stage, the uplink information transmitted between EDs and BS is the EDs' state space s 1→K (t) and action space a 1→K (t).
Moreover, the structure of the actor network for agent EDs and agent BS are illustrated in Fig. 4 (b) and (c), respectively. It can be observed that the defined heterogeneous multi-agent has separate actor network structures from input, activation function, mapping and output. The pseudocode of the training process of the proposed distributed RL algorithm named as HMA-LMIDDPG can be found in Algorithm 2. Although the actor and critic network parameters training process for each agent is similar to the centralized framework in Algorithm1, there are two differences we would like to emphasize. Firstly, as the main purpose of our work is to acquire optimal energy efficiency of the whole system with queue stability, the reward is the cooperative reward which is calculated with the knowledge of the whole system. Therefore, during the training stage, the cooperative Lyapunov-function-based reward is calculated at the BS and then broadcasted to all EDs for training. Secondly, in the designed heterogeneous multiagent structure, each agent i ∈ {1, ..., K + 1} has its own replay buffer memory M i to store their individual experience {s   The distance-dependent path-loss factor equal to β (DI) , and d (IB) are the distances of ED u k to IRS, ED u k to BS, and IRS to BS links, respectively. Moreover, the Racian factor δ = 1 and the channel noise power is σ 2 = −174dBm [36]. The summary of other system simulation parameters is listed in Table I. The parameters of the critic and actor networks in proposed LMIDDPG and HMA-MIDDPG can also be found in Table I. To show the superior energy efficiency performance of our proposed algorithms, there are several benchmarks we have compared: • No IRS: there is no IRS in the MEC system, • All Offloading + Random Phase: all the users decide to offload and the phase shifts on the IRS are randomly assigned, • Random Phase: the phase shifts on the IRS are randomly assigned, • LMIA2C: the A2C RL algorithm [33] with Lyapunov-function-based reward and mixed integer mapping.

A. Benchmarks Comparison
As shown Fig. 5 and Fig. 6, the convergence performance and the energy efficiency have been compared among different benchmarks, respectively. It can be observed that in Fig. 5

System Parameters Training Parameters Parameters
Value [19] 10 −26 , 10 −27 L ac ,L cr 5, 5 c k , c 0 100, 100 cycles/bit n ac,1→(Lac−1) 256 V 10 n cr,1→(Lcr−1) 256  and algorithms show higher system energy efficiency, which proves the benefits of deploying IRS in the MEC system. For example, the random phase IRS can achieve around 25% energy efficiency enhancement. The proposed LMIDDPG has around 56% improvement. Furthermore, the proposed distributed framework HMA-HMIDDPG can achieve further 10% improvement compared to centralized framework LMIDDPG.
In Fig. 7, we further compare the average system throughput over T = 300 and the energy efficiency at T = 300 based on different IRS element sizes M. It can be observed in Fig. 7 (a), the All Offloading scenario always achieves higher system throughput due to the constant task offloading transmission. However, the energy efficiency of which is very low as shown in Fig. 7 (b). Differently, the proposed algorithms, LMIDDPG and HMA-MIDDPG, can acquire the highest energy efficiency, which proves the superior resource allocation decision. Moreover, in Fig. 7 (a), apart from all offloading benchmark, it can be observed that the lower throughput proves the different offloading and local computing decisions have been made. What's more, with lower system throughput, the energy consumption is much lower to achieve higher energy efficiency performance.

B. Impact of Different Parameters
In Fig. 8, the average system throughput over T = 300 and the energy efficiency at T = 300 based on different user number K have been compared among the proposed centralized framework LMIDDPG and distributed framework HMA-LMIDDPG. As shown in Fig. 8 (a), with the user number increasing, the overall throughput increases accordingly under both frameworks.
Moreover, from Fig. 8 (b), HMA-LMIDDPG can achieve higher system throughput with higher energy efficiency, which proves the superior decision making based on the proposed distributed framework.
In Fig. 9, the average system throughput over T = 300 and the energy efficiency at T = 300 based on different task generation probability ζ k have been compared among the proposed centralized framework LMIDDPG and distributed framework HMA-LMIDDPG. It is easily noticed in Fig. 9 (a) that with higher task generation probability, the system throughput becomes  higher. However, Fig. 9 (b) shows that the energy efficiency has similar performance among different probabilities.
As shown in Fig. 10 and Fig. 11, the convergence performance, the average system throughput over T = 300 and the energy efficiency at T = 300 based on different Lyapunov weights V have been illustrated, respectively. In Fig. 10, compared to LMIA2C, the proposed algorithms can converge at higher rewards. Moreover, with V increases the rewards decreases. This trend is because the design of the reward function is the negative eq. (22) with negative V . Moreover, from Fig. 11 (a) that with V increases, the system throughput also increases under centralized framework LMIDDPG. From Fig. 11 (b), the proposed distributed framework HMA-LMIDDPG always achieves better energy efficiency performance.  In Fig. 12 and Fig. 13, the convergence performance, the average system throughput over T = 300 and the energy efficiency at T = 300 based on different arrival data sizes A k (t) are plotted, respectively. It is easily noticed in Fig. 12 that with A k (t) increases, the rewards decreases. This is because of the design of the reward function is the negative eq. (22) with upper bounds obtained from eq. (29) with arrival data. Moreover, as shown in Fig. 13, with arrival data size increases, the distributed framework HMA-LMIDDPG can always achieve higher average system throughput and higher energy efficiency compared to LMIDDPG.  Fig. 15 (a). However, from Fig. 15 (b), compared to LMIDDPG, the distributed HMA-LMIDDPG maintains higher system throughput with time evolving.

VI. CONCLUSIONS
In this paper, we investigated the joint offloading, communication and computation resource allocation for IRS-assisted NOMA-aided MEC system. We proposed the mixed integer energy efficiency maximization problem with system queue stability constraint. The centralized RL framework algorithm called LMIDDPG has been proposed with the award function defined as the upper-bound of the Lyapunov drift-plus-penalty function. In addition, we further designed the mixed integer action space mapping which contains both continuous mapping and integer mapping. To ensure less communication interaction overheads between BS and EDs during the execution stage, HMA-LMIDDPG has been further proposed with heterogeneous multi-agent that includes both homogeneous agents EDs and heterogeneous agent BS. Under this framework, each agent can make its individual action decisions based on its local observation. Numerical results clearly proved that, our proposed algorithms can greatly increase the system energy efficiency with maintained queue stability compared to benchmark algorithms. Moreover, the distributed structure HMA-LMIDDPG can achieve more energy efficiency gain compared to centralized structure LMIDDPG.

APPENDIX A LYAPUNOV DRIFT-PLUS-PENALTY UPPER BOUND
The drift item in eq.(22) can be extended as ∆L (Q(t)) = 1 2 Q 2 (t + 1) − Q 2 (t) = 1 2 Q(t) + Q(t) (in) − Q(t) (out) 2 − Q 2 (t) Specifically, b is a constant upper bounds which can be obtained as where the first two items in the inequality denotes the maximum input and the processed (i.e. local computing and offloaded) data size at K EDs task queues, the last two items represents the maximum arrived and the processed data size in the edge task queue at the BS.