{'fraction_non_alphanumeric': 0.06130922817125261, 'fraction_numerical': 0.01648952021456243, 'mean_word_length': 4.295739084692268, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multicast Scheduling for Multi-Message over Multi-Channel: A Permutation-based Wolpertinger Deep Reinforcement Learning Method

Multicasting is an efficient technique to simultaneously transmit common messages from the base station (BS) to multiple mobile users (MUs). The multicast scheduling problem for multiple messages over multiple channels, which jointly minimizes the energy consumption of the BS and the latency of serving asynchronized requests from the MUs, is formulated as an infinite-horizon Markov decision process (MDP) with large discrete action space and multiple time-varying constraints, which has not been efficiently addressed in the literatures. By studying the intrinsic features of this MDP under stationary policies and refining the reward function, we first simplify it to an equivalent form with a much smaller state space. Then, we propose a modified deep reinforcement learning (DRL) algorithm, namely the permutation-based Wolpertinger deep deterministic policy gradient (PW-DDPG), to solve the simplified problem. Specifically, PW-DDPG utilizes a permutation-based action embedding module to address the large discrete action space issue and a feasible exploration module to deal with the time-varying constraints. Moreover, as a benchmark, an upper bound of the considered MDP is derived by solving an integer programming problem. Numerical results validate that the proposed algorithm achieves close performance to the derived benchmark.


I. INTRODUCTION
The number of mobile users (MUs) and mobile data traffic keep increasing rapidly in recent years. In 2023, the number of smartphones is projected to grow from 4.9 billion in 2018 to 6.7 billion [1] and the mobile data traffic is expected to reach 130 exabytes per month, which is 5-fold of the monthly data traffic in 2018 [2]. Remarkably, video data traffic, including live streaming, virtual reality (VR) video, and augmented reality (AR) video, which has a high similarity among the contents requested by different MUs, is expected to make up 79 percentage of the total mobile data traffic in 2023 [2]. The current cellular systems can hardly fulfill such heavy tasks with the unicast technique for the downlink transmissions, by which the base station (BS) assigns individual channel to each MU's request. By exploiting the similarity of the data requests, multicast technique can temporarily hold the MU's requests until there are sufficient amount of them, and then serve them simultaneously during one multicast transmission. Apparently, exploiting multicast in the future cellular network could significantly improve the energy efficiency (defined as the average energy cost per request) [3]. Notably, the asynchronous requests from multiple MUs would be served with different latencies in this multicast mechanism. Therefore, it is of great interest to balance the energy efficiency and the average latency of the considered multicast scheme in the cellular networks.

A. Related works
In recent years, many works have studied the multicast scheduling problem in the cellular networks, especially in the 5G era where the cellular networks are assisted by the advanced techniques such as millimeter wave (mmWave) communication [4], non-orthogonal multiple access (NOMA) [5], and intelligent reflecting surface (IRS) [6]. One major goal is to maximize the system throughput with limited resources, and approaches including deep reinforcement learning (DRL), beamforming, were validated to have promising performances in [4], [6]- [8]. Along another research avenue, due to the increasing public attention to green communications, energy efficiency maximization becomes more compelling in the recent researches [5], [9]- [11] and aims to serve more requests with as little energy as possible. To tackle this issue, the authors in [5], [9] utilized beamforming and designed the optimal precoders at the BS to maximize the energy efficiency for the NOMA and mmWave systems, respectively. The authors in [10] discussed the device-to-device assisted multicast scheduling problems for the mmWave system, and proposed a supervised learning method to maximize the energy efficiency. The authors in [11] discussed the joint unicastmulticast scheduling for the NOMA system, and proposed a rate-splitting method to optimize both the spectral and energy efficiency. Another major concern for multicast scheduling is the fairness issue, which aims to balance the latencies among different MUs [12], [13]. Obviously, the methods only optimizing the energy efficiency in [5], [9]- [11] prefer to allocate resources to the MUs with better channel condition, whereas the MUs with poor channel condition would be served with more latencies and may even be abandoned by the BS. To jointly optimize the energy efficiency and the fairness, namely the proportional fairness [14], a wise choice is to optimize the summation of the logarithm transmission rates of all MUs [15], where, due to the "strictly concave" and the "increasing" features of the logarithm function, improving the already high transmission rate for the good-channel-condition MUs is inferior to serving the MUs with poor channel condition. Another approach is to gather the MU's requests into request queues [16]. By utilizing the Lyapunov optimization, the rate stability of the request queues was established, which ensures finite latencies for the MUs' requests. However, the proposed algorithms in [15], [16] only gave a rough control on the latency to be finite, and the precise dominations, such as minimizing the latency subject to certain constraints, are impossible under the considered framework.
To jointly optimize the energy efficiency and the latency in a more flexible way, the most common approach is to minimize the weighted summation of them [17]- [20]. However, since the information about the accumulated requests from the MUs is required to compute the latency, the multicast scheduling problem was formulated as an infinite-horizon Markov decision process (MDP), which is a commonly known difficult problem due to "the curse of dimensionality" [21]. The authors in [17]- [20] proposed some low-complexity methods to solve this problem. Specifically, the authors in [17] studied the multicast scheduling problem where multiple MUs request one common message and the MUs have mixed latency penalties, and a suboptimal scheduling policy was developed based on the optimal stopping theorem. The authors in [18] discussed the multicast scheduling problem with multiple requested messages and single available channel, and the optimal scheduling policy was constructed based on the classical relative value iteration (RVI) algorithm. The authors in [19], [20] deployed DRL algorithms to study the joint multicast scheduling and message caching problems in the ultra dense networks and the heterogeneous networks, respectively, and again focused on the scenarios with multiple requested messages and single available channel. Notably, the authors in [17]- [20] assumed that the multicasting of each message on each channel starts and ends synchronously, i.e., all channels are always released from multicastings simultaneously and available at every scheduling moment.
To the best of our knowledge, these is no existing work discussed the scheduling problem for asynchronous multicasting, where the multicastings of different messages over different channels consume different time durations and thus different multicastings would start and end in an asynchronous manner.

B. Main contributions
This paper studies the most general multicast scheduling problem in the cellular networks: a) multiple MUs randomly request multiple messages, which are stored at the BS in advance with the cache technique [19], [20]; b) multiple channels are available for multicast transmissions at the BS; c) the multicastings of different messages over different channels consume different time durations, i.e., asynchronous multicasting. d) the optimal multicast scheduling policy is supposed to strike the optimal tradeoff between the energy efficiency and the latency penalty. We summarize our contributions as follows: • We formulate the multicast scheduling problem as an infinite-horizon MDP, which is challenging due to the large discrete action space and multiple time-varying constraints induced by stochastic availability of channels. By studying the intrinsic features of the considered MDP under stationary policies and refining the reward function, we first simplify this MDP to an equivalent form with a much smaller state space. Then, we propose a modified DRL algorithm based on deep deterministic policy gradient (DDPG) [22], which is named as the permutationbased Wolpertinger DDPG (PW-DDPG), to address the simplified MDP. Compared with the conventional DDPG, PW-DDPG contains a permutation-based action embedding module and a feasible action exploration module, where the former module embeds actions in the "permutation" manner during both the offline training and online validating procedures of PW-DDPG and can well address the large discrete action space issue. To address the timevarying constraints issue, the two modules separately restrict the exploitation and exploration modes of PW-DDPG to only selecting feasible actions. • We derive an upper bound of the considered MDP, which also works as the benchmark to validate the proposed PW-DDPG. Specifically, we first prove that the timevarying constraints of the considered MDP can be released as a set of time-invariant constraints, which reveal the statistical features of the considered MDP under the policies satisfying the time-varying constraints. Then, we convert the released problem to a two-step optimization problem, where the optimization problem of the first step purely minimizes the latency penalty and the second step jointly minimizes the energy consumption and the latency penalty. Finally, we successively solve the twostep optimizations and obtain an upper bound of the considered MDP by solving an integer programming problem. The remainder of this paper is organized as follows. Section II introduces the system model and formulates the multicast scheduling problem. Section III presents the proposed algorithm. Section IV discusses the upper bound benchmark of the multicast networks. Section V evaluates the performance of the proposed algorithm. Finally, Section VI concludes this paper.

II. SYSTEM MODEL AND PROBLEM FORMULATION
Consider a slotted cellular network as depicted in Fig. 1, which consists of one BS that caches N messages and a set of E MUs, denoted by E = {MU 1 , · · · , MU E }. During each time slot, a subset E n ⊆ E of the total MUs with size E n = |E n | are likely to send the request to the BS for downloading the n th message with probability α n . We consider the case that each MU only has interest in downloading one message and thus we have E n ∩ E m = ∅, ∀n = m, and E =   Figure 1: Schedule over the multi-message and multi-channel multicast system. of the new requests for N messages in N request buffers. Instead of immediately multicasting the messages requested by the buffered requests, the BS may wait for several time slots and simultaneously multicast these messages until there are sufficient number of requests in the request buffers. We also consider that M downlink channels are available for multicast transmissions and the multicastings of N messages may last for multiple consecutive time slots depending on their sizes. We formulate the multicast scheduling problem for the above multi-message and multi-channel multicast system as an infinite-horizon MDP described as follows.

Request buffers
1) State: The state contains request matrix, channel availability, and channel status.
Request matrix: Denote the number of the requests for the n th message arriving within the t th time slot as q n (t), which follows binary distribution as q n (t) ∼ Binary(E n , α n ), and the number of time slots after the previous multicasting of the n th message as l n (t) ∈ N, where N is the set of 0 and all natural numbers. Then, within the past l n (t) time slots, all the arrived requests for the n th message are stored in the n th request buffer of the BS, and we represent these requests by the request vector q n (t), which is defined as Here, M * is the size of the request buffer and if l n (t) > M * , the requests arriving before the (t − M * ) th time slot are accumulated in the last entry of q n (t). Finally, define the request matrix as Q(t) [q 1 (t), q 2 (t), · · · , q N (t)] T , which has the dimensionality of N -by-M * . Channel availability: Each multicasting of the n th message over the m th channel is considered to consume T n,m ∈ Z + consecutive time slots with Z + being the set of all positive integers, during which the m th channel is not available for new multicasting. Denote the availability of the m th channel for multicasting the n th message in the beginning of the t th time slot as c n,m (t) ∈ {0, 1, · · · , T n,m − 1}: if the BS is not multicasting the n th message over the m th channel in the beginning of the t th time slot, c n,m (t) is set to 0; otherwise, c n,m (t) is equal to the number of the remaining time slots for the release of the m th channel. Then, define the channel availability as a N -by-M -dimension matrix C(t), where the (n, m) th entry of C(t) is c n,m (t), i.e., [C(t)] (n,m) c n,m (t).
Channel status: Denote the number of MUs to be served in the beginning of the t th time slot by multicasting the n th message as K n (t) ∈ N, where K n (t) = M * τ =1 [q n (t)] τ holds with [q n (t)] τ being the τ th entry of q n (t), and the indices of these MUs as I n,1 (t), I n,2 (t), · · · , I n,Kn(t) (t) ∈ {1, 2, · · · , E n }. Denote the downlink channel coefficient of the link between the BS and the I n,k (t) th MU and over the m th channel as h n,m,I n,k (t) ∈ C, where 1 ≤ k ≤ K n (t) and C is the set of all complex numbers. Then, the smallest downlink channel gain among K n (t) BS-MU links is given as Finally, we denote the channel status as a N -by-M -dimension matrix G(t), with the (n, m) th entry of G(t) being g n,m (t), i.e., [G(t)] (n,m) g n,m (t).
To summary, the state of the considered system at the t th time slot is the triple s(t) (Q(t), C(t), G(t)) and the state space S is S Q N ×M * ×C N ×M ×G N ×M with dimensionality of N M * + 2N M , where Q, C, and G are the value spaces of [q n (t)] τ , c n,m (t), and g n,m (t), respectively.
2) Action: Denote the multicast scheduling decision for the m th channel in the beginning of the t th time slot as a m (t) ∈ {0, 1, · · · , N }. Specifically, a m (t) > 0 means to multicast the a m (t) th message over the m th channel in the beginning of the t th time slot; a m (t) = 0 means that the m th channel is not scheduled to multicast any message. Denote the multicast scheduling decision for all channels as a(t), i.e., a(t) [a 1 (t), a 2 (t), · · · , a M (t)] T , and obviously, a(t) is the action of the considered system. Notably, if the m th channel has been reserved for multicasting in previous slots, i.e., N n=1 c n,m (t) > 0, it cannot be scheduled to multicast any new message. That is, the action a(t) is constrained by N n=1 c n,m (t)a m (t) = 0, m ∈ {1, 2, · · · , M } . (3) Accordingly, we denote the action space at state s(t) as A s(t) , which contains all the actions satisfying the time-varying constraints in (3).

3) Transitions:
The transitions are to update request matrix, channel availability, and channel status. First, we define the multicasting status of the n th message at the t th time slot as b n (t) I M m=1 I n (a m (t)) , n ∈ {0, 1, · · · , N }, (4) where I n (x) equals 1 if x is n, and otherwise, equals 0; I(x) equals 1 if x is positive, and otherwise, it equals 0. Apparently, if the n th message is selected for multicasting over some channels in the beginning of the t th time slot, b n (t) equals 1, and otherwise, it equals 0. Then, based on b n (t), the transitions are derived as follows.
Request matrix: We update the request matrix in the following two cases: if the n th message is selected for multicasting in the beginning of the t th time slot, i.e., b n (t) = 1, only q n (t) new requests will be buffered in the beginning of the next time slot and q n (t + 1) = [q n (t), 0, · · · , 0] T ; if the n th message is not selected for multicasting, i.e., b n (t) = 0. The buffered requests in the beginning of the next time slot contain both the previously buffered ones and the newly arrived ones, i.e., q n (t where [q n (t) T ] 1:M * −2 is the 1-by-(M * − 2)-dimension vector containing the first to the (M * − 2) th entries of q n (t) T . To summary, we have Channel availability: We update the channel availability in the following three cases: if the m th channel has been reserved for multicasting the n th message in previous slots, i.e., c n,m (t) > 0. The remaining time for the release of the m th channel decreases by one in the beginning of the next time slot, i.e., c n,m (t + 1) = c n,m (t) − 1; if the m th channel has not been reserved for multicasting the n th message and is going to multicast the n th message in the beginning of the t th time slot, i.e., c n,m (t) = 0 and a m (t) = n, the m th channel will be released after T n,m − 1 time slots counting from the (t + 1) th time slot, i.e., c n,m (t + 1) = T n,m − 1; and if the m th channel has not been reserved for multicasting the n th message and is not going to multicast the n th message, i.e., c n,m (t) = 0 and a m (t) = n, c n,m (t + 1) remains 0. To summary, we have Channel status: We model the channel coefficient between the BS and any MU and over any channel as a stationary ergodic process, i.e., where Pr n,m,k {h j |h i } is a constant and h i , h j ∈ H with H being a finite complex number set. Then, g n,m (t+1) is derived as 4) Reward: Both the energy consumption and the latency penalty are the concerned performance metrics.
Energy consumption: The energy consumption of multicasting the n th message over the m th channel within the t th time slot is is a constant with C, B, and τ being the information bits of the n th message to be transmitted within one time slot, the bandwidth of the m th channel, and the duration of one time slot, respectively. Note that the above energy consumption exists only if a m (t) = n or c n,m (t) > 0 holds, when the m th channel starts to multicast the n th message or has been reserved for multicasting the n th message. Thus, the total energy consumption within the t th time slot is Zn,m gn,m(t) (I n (a m (t)) + I(c n,m (t))). Latency penalty: The requests buffered at the BS will produce instant latency penalty in every time slot before they get served. Denote a function p n (τ ) ∈ R as the instant latency penalty for delaying one request for the n th message over τ time slots. Then, according to (1), the buffered requests in q n (t) will produce instant latency penalty of M * τ =1 [q n (t)] τ p n (τ ) at the t th time slot. Thus, the total instant latency penalty produced at the t th time slot is N n=1 The reward r(t) of the considered system is defined as the weighted sum of the energy consumption and the latency penalty, i.e., where V is the tradeoff parameter. 5) Problem formulation: We aim to minimize the long-term average reward and thus formulate the multicast scheduling problem as where the expectation is taken with respect to the request arrival process {q n (t)} N n=1 and the transitions on channel status, i.e., Pr{G |G}. To solve problem (P1), we need to find the optimal policy π (π 1 , · · · , π t , · · ·), where π t is the optimal scheduling rule at the t th time slot and maps the history h(t) (s(1), a(1), · · · , s(t − 1), a(t − 1), s(t)) to the optimal distribution of a(t), i.e., π t : H(t) × A s(t) → [0, 1] with H(t) being the set of all histories h(t). Apparently, the scope of the required history for π t grows linearly with respect to t and to derive such history-based optimal policies is too difficult [21]. Therefore, we only investigate the stationary policies for problem (P1), which determines the action a(t) based on the current state s(t) and utilizes the same scheduling rule, say π, at every time slot, and try to derive the stationary policy that achieves with Π S being the set of all stationary policies. This modified problem, i.e., problem (P1) with solution space containing only the stationary policies, is named as problem (P1-1), and to solve it, we first simplify it to another MDP (P2): where (11) is derived by (3) where (12) is derived by (6). • Reward:r(t), which is defined aŝ with u(n, m, g) Here, g n,m (τ )| gn,m(0)=g is a random variable representing the value of g n,m (τ ) given the value of g n,m (0) as g and its distribution can be easily derived according to Pr{G |G}. • Problem formulation: problem (P2) is formulated as (7), (8), (11), (12).
Compared with problem (P1-1), problem (P2) has a much simpler structure: the dimensionality of the state space is reduced from N M * +2N M to N M * +M +N M ; the reward function no longer involves the high-dimension matrix C(t). Moreover, we obtain the following proposition. Proposition 2.1: Problems (P1-1) is equivalent to problem (P2), i.e., their optimal values are equal.
Proof: Please see Appendix A.
Remark 2.1: Problem (P2) is still challenging due to its three features: 1) large state spaceŜ; 2) large discrete action space A s ; and 3) multiple time-varying constraints on the actions in (11), and existing approaches cannot efficiently address it.
(1) Dynamic programming [21] requires value evaluation for each state s, which is defined as J(s) given s as the initial state. Therefore, it can barely deal with the MDP with very large state space, such as problem (P2). Lyapunov optimization determines the action at each time slot by solving a onestep optimization problem, which is only related to the current state and thus has a low computational complexity. Notably, to deploy the Lyapunov optimization for MDP, a queue model is first formulated and the independency of the "departure variables" over time is required. However, the departure variables of problem (P2) are the numbers of the served MUs at each time slot, which based on (5), are not independent over time.
(2) DRL algorithms cannot efficiently address the MDPs with large discrete action space or with time-varying constraints on the actions. First, only a few literatures have discussed the MDP with large discrete action space, whereas the most valuable attempt is the Wolpertinger policy proposed in [24] and it requires the MDP satisfying that "if two policies always select actions close in L 2 distance at each state, the MDP would transit to close states and feedback close rewards, and thus the two policies have similar performances". However, the MDP in problem (P2) does not have this feature, which can be easily verified by checking the state transitions in (5), (8), (12) and the reward function in (13). Moreover, almost all existing DRL algorithms utilize the exploitation versus exploration mechanism, where both the exploitation and the exploration are unpredictable during the training phase and should randomly select action. Therefore, for our proposed MDP with multiple time-varying constraints on the actions, it is of high probability for the DRL to select infeasible actions, and consequently, terminate the training without converging to a stable policy. Thus, DRL algorithms cannot be directly applied to our problem (P2).

III. PERMUTATION-BASED WOLPERTINGER DEEP DETERMINISTIC POLICY GRADIENT
To solve problem (P2), this section proposes PW-DDPG, which mainly comprises of DDPG, the permutation-based action embedding (PAE) module, and the feasible action exploration (FAE) module. Particularly, the PAE module utilizes    the action embedding technique [24] to address the large discrete action space issue. To deal with the aforementioned drawback of the conventional action embedding in remark 2.1, the proposed action embedding method in PAE module has a permutation-based structure and is perfectly compatible with the MDP problem (P2). To address the time-varying constraints issue, the PAE module restricts the exploitation mode to only selecting the feasible actions and the FAE module utilizes a modified -greedy method to guarantee that the exploration mode also only explores the feasible actions.

Generation of experiences
In the following, we first introduce the structure of the proposed algorithm in details. Then, we present the offline training algorithm, which is based on the interaction with a simulated offline environment and aims to train the NNs embedded in the PW-DDPG. Finally, we illustrate the online validating algorithm to solve problem (P2).

A. Structure of proposed algorithm
As illustrated in the left part of Fig. 2, the proposed algorithm consists of four modules: data refinement (DR) module, DDPG, PAE module, and FAE module.
2) DDPG: DDPG consists of the actor network, the critic network and the experience replay module, and it is deployed to generate the proto-actionâ(t) based on the refined statê s(t), where the proto-action will be specified soon.
• Actor network: The actor network contains an evaluation NN and a target NN, which are parameterized by θ 1 , θ 1 and denoted as π θ1 , π θ 1 , respectively. Particularly, two NNs take the refined stateŝ(t) as inputs and thus have N M * +M +N M nodes at their input layers. The outputs of them are set to be with dimensionality of M , which equals the dimensionality of action a(t), and thus their output layers have M nodes. Finally, we define the output of the actor evaluation NN as the proto-action and denote it asâ(t) ∈ R M×1 , i.e., where π θ1 (ŝ(t)) is the output of the actor evaluation NN with input beingŝ(t). • Critic network: The critic network also contains an evaluation NN and a target NN, which are parameterized by θ 2 , θ 2 and denoted as Q θ2 , Q θ 2 , respectively. Particularly, two NNs take the state-action pairs (ŝ(t), a(t)) as inputs and thus have N M * + 2M + N M nodes at the input layers. Their output layers have only one node, which represents the evaluated Q values of the state-action pairs. • Experience replay module: The experience replay module contains the experience buffer and the mini-batch. Specifically, the experience buffer stores the generated experiences during the offline training phase, which are the four-component tuples (ŝ(t), a(t),r(t),ŝ(t+1)). The mini-batch randomly samples N 0 experiences from the experience buffer to update NNs, and these sampled experiences are denoted as 3) PAE module: This module takes the proto-actionâ(t) as the input and outputs the Wolpertinger-actionā(t), where the Wolpertinger-action will be specified soon. As illustrated in Fig. 3(b), it contains the following three procedures: • Rounding and feasibility fitting: Given the proto-action a(t), we first round the entries inâ(t) to their closest integers and store the results inâ (1) (t), i.e.,â , (t Q θ s 2ˆ) , (t Q θ s 2 L Figure 3: Conventional action embedding [24] vs. the proposed permutation-based action embedding. whereâ (2) m (t) is the m th entry ofâ (2) (t). Notably, althoughâ (2) it may still violate the time-varying constraints in (11) and thus be infeasible. Thus, we finally check the feasibility of a (2) (t) by validating the constraints in (11) and figure out a feasible actionâ (3) (t) accordingly. Particularly, denote the m th entry ofâ (3) (t) asâ It is easy to check thatĉ m (t)â m (t) = 0 holds, and thuŝ a (3) (t) is feasible. We also callâ (3) (t) as the candidateaction. • Action embedding based on permutation: Based on the candidate-actionâ (3) (t), we design action embedding to derive more feasible actions, which are called as the permutation-actions and denoted asâ (4) (t). Specifically, we first denote the index set of the currently available channels as B(t) {m|ĉ m (t) = 0, m ∈ {1, 2, · · · , M }} and gather the entries ofâ (3) (t) at these indices in another set C(t) {â Obviously, the generated permutation-actionâ (4) (t) satisfies the constraints in (11) and thus is feasible. Denote the set of all possible permutation-actions as D(t) and apparently, D(t) contains up to |B(t)|! permutationactions. Then, we randomly pick out k permutationactions from D(t), where the i th picked out permutationaction is denoted asâ (4) i (t), i ∈ {1, 2, · · · , k}, and send them to the next procedure. • Action selection based on Q θ2 andŝ(t): Send the refined stateŝ(t) together with k picked out permutationactions {â (4) i (t)} k i=1 to the critic evaluation NN Q θ2 , and then k outputs are obtained, which are denoted as {Q θ2 (ŝ(t),â (4) i (t))} k i=1 . The permutation-action with the largest output will be selected as the Wolpertinger-action a(t), i.e.,ā (t) arg max i∈{1,2,···,k} Q θ2 (ŝ(t),â i (t)).
Remark3.1: The conventional action embedding is proposed in [24] and illustrated in Fig. 3(a). Specifically, it first obtains k nearest discrete actions of the proto-action in the L 2 distance manner, and then send these actions together with the refined stateŝ(t) to the critic evaluation NN. Finally, the action with the largest output will be selected as the Wolpertingeraction. However, if the MDP is with multiple constraints, the k nearest discrete actions of the proto-action are very likely to be infeasible, and so is the obtained Wolpertingeraction, which disables the offline training procedure. In the proposed permutation-based action embedding, however, we replace the k nearest actions found in the L 2 distance manner by k permutation-actions, which are near to the proto-action in the "permutation" manner. As we validated before, these permutation-actions are always feasible. Moreover, it can be easily checked that "if two policies always select actions close in the "permutation" manner at each state, the MDP problem (P2) would transit to close states and feedback close rewards, and thus the two policies have similar performances". The proposed permutation-based action embedding takes advantage of this fact and thus is perfectly compatible with the MDP problem (P2).

4) FAE module:
This module takes the Wolpertinger-action a(t) as the input and outputs the action of the considered MDP, i.e., a(t). Specifically, it sets a(t) asā(t) with probability of 1 − and as a random feasible action with probability of , i.e., with probability 1− a random action satisfying the constraints in (11) with probability .
Compared with the action exploration methods of conventional DRL algorithms, such as -greedy exploration introduced in deep Q network (DQN) [25], Gaussian noise exploration introduced in DDPG [22], and entropy-based exploration introduced in proximal policy optimization (PPO) [26], the proposed FAE module does not explore actions in a purely random manner and can avoid selecting infeasible action based on (21).

B. Offline training
We first simulate an offline environment based on the historical observations on the MUs' requests and the channel status. Then, we alternatingly generate new experiences by utilizing PW-DDPG with the latest NN parameters to interact with the offline environment and update the NN parameters in PW-DDPG based on the latest generated experiences. The sketches of the experiences generation and the PW-DDPG update are illustrated in the right part of Fig. 2. 1) Offline environment simulation: The simulated offline environment is supposed to simulate the value of the new state s(t + 1) based on the values of s(t) and a(t). Specifically, we fist estimate the number of MUs {E n } N n=1 , the probability of MUs' requesting messages {α n } N n=1 , and the channel coefficient transition Pr n,m,k (h j |h i ) based on the historical observations, and simulate the values of {q n (t)} N n=1 and Pr{G |G} accordingly. Then, given the values of s(t) and a(t), the new state s(t + 1) can be derived by (5), (6), (7), and (8).
2) Generation of experiences: First, obtain the value of s(t) from the simulated offline environment and send it to the DR module, by which the refined stateŝ(t) is derived. Next, sendŝ(t) to the actor evaluation NN and obtain the protoactionâ(t). Then, by deploying the PAE module, derive the Wolpertinger-actionā(t). Finally, derive the action a(t) by utilizing the FAE module. With the derived a(t), the refined rewardr(t) can be derived by sending s(t) and a(t) to the data refinement module, the new state s(t+1) can be obtained from the simulated offline environment, and the refined new stateŝ(t + 1) can be derived by sending s(t + 1) to the DR module. The experience packs the above information into the tuple (ŝ(t), a(t),r(t),ŝ(t + 1)).
Algorithm 1 Offline training of PW-DDPG for multicast scheduling 1: Randomly initialize θ 1 and θ 2 for the actor and critic evaluation NNs; 2: Initialize actor target NN and critic target NN by θ 1 ← θ 1 and θ 2 ← θ 2 ; 3: Initialize replay buffer R with size of N R ; 4: Simulate the offline environment; 5: for episode= 1, 2, · · · for t = 1, 2, · · · , T 10: Sendŝ(t) to the actor evaluation NN π θ1 and derivê a(t); 11: Sendâ(t) to the PAE module and deriveā(t); 12: Sendā(t) to the FAE module and derive a(t); 13: Deriver(t) from DR module withŝ(t) and a(t) as the input; 14: Obtain s(t + 1) from the offline environment with a(t) as the input; 15: Send s(t+1) to the DR module and deriveŝ(t+1); 16: Store (ŝ(t), a(t),r(t),ŝ(t + 1)) in the experience buffer; 17: if t > N R − 1 18: Use mini-batch to sample N 0 experiences from the experience buffer; 19: Update actor evaluation NN based on (22); 20: Update critic evaluation NN based on (23); 21: Update the target NNs based on soft update rule; 22: end if 23: end for 24: end for 3) Update of the NN parameters in PW-DDPG: After each generation of one new experience, the experience buffer is refreshed, i.e., one old experience is replaced by the new one. Then, we use mini-batch to sample N 0 experiences from the experience buffer as {(ŝ i , a i ,r i ,ŝ i )} N0 i=1 and utilize them to update the NN parameters in PW-DDPG. Specifically, the actor evaluation NN is updated based on the policy gradient theorem, i.e., the critic evaluation NN is updated by backpropogating the following loss function where γ ∈ (0, 1) is a discount factor, and two target NNs utilize soft update rule [22]. The offline training algorithm is summarized in Algo-rithm 1.

C. Online validating
After the offline training, we design the online validating algorithm to solve problem (P1), which is very similar to the offline training algorithm. Specifically, we first derive the refined stateŝ(t) based on the DR module and the observed state s(t) from the practical multicast environment. Then, we sequentially utilize the PAE module and the FAE module to derive the action a(t). Finally, we execute the action a(t).

IV. DISCUSSION
In this section, we derive the upper bound of the objective function of problem (P2), and it works as the benchmark to validate the proposed PW-DDPG in the simulation part. To derive the upper bound, we first release the time-varying constraints of problem (P2) in (11) to a set of time-invariant constraints, which reveal the statistical features of problem (P2) under stationary policies satisfying the constraints in (11). Then, we convert the released problem to a two-step optimization problem, where the optimization of the first step purely minimizes the latency penalty and the second step jointly minimizes the energy consumption and the latency penalty. Finally, we successively analyze the optimizations of these two steps, and obtain an upper bound for problem (P2).
1) Release on constraints: We release the time-varying constraints in (11) with the following proposition. where γ n,m (T ) ∈ N counts how many times the n th message has been multicasting over the m th channel within the 1 st to T th time slots and is defined as with n ∈ {1, 2, · · · , N }, m ∈ {1, 2, · · · , M }, and T ∈ Z + .
Proof: Please see Appendix B. Proposition 4.1 indicates that the constraints in (24) are looser than those in (11). Thus, by replacing the constraints in (11) by (24), we derive a released version of problem (P2) as where the transitions in (12) from problem (P2) are dropped in problem (P3) since without the constraints in (11), (12) is no longer needed in problem (P3), and we also release the solution space such that the actions {a(t)} may not necessarily be picked following stationary polices.
Denote the upper bound of problem (P4) as f (γ), whereγ is a N -by-M -dimension matrix and defined by [γ] (n,m) γ n,m . Then, we jointly optimize the energy consumption and the latency penalty terms of problem (P3) by further solving the following problem: n,m T n,m ≤ 1, m ∈ {1, 2, · · · , M }, (26) where the constraints in (26) is derived by (24). Since the optimal solutions to problems (P4) and (P5) may assign different values to the design variable {a(t)}, problem (P5) is a released version of problem (P3) and thus the upper bound of problem (P5) is also an upper bound of problem (P3). In the following, we successively analyze problems (P4) and (P5) and derive their upper bounds, respectively.
3) Upper bound of problem (P4): By splitting both the objective function and the constraints of problem (P4), we divide problem (P4) into N sub-problems, where the n th subproblem is formulated as [q n (t)] τ p n (τ ) s.t. (5) γ n,m (T ) =γ n,m T, m ∈ {1, · · · , M } . (27) Note that the optimal solutions to N sub-problems may assign different values to {a(t)}. Thus, the summation of the optimal values of N sub-problems is an upper bound of problem (P4).
Problem (P4-1) can be interpreted as: if the number of the multicastings for the n th message is fixed as M m=1γ n,m T within the 1 st to T th time slots, when should we start these multicastings to minimize the average latency penalty?
To answer this question, we first denote t i with i ∈ {1, 2, · · · , M m=1γ n,m T } as the beginning time of the i th multicastings of the n th message, i.e, b n (t i ) = 1 holds, and without loss of generality, we let t M m=1γ n,m T = T . Then, it [q n (t)] τ p n (τ ), (28) where t 0 0 holds. Based on the above equation, the objective function of problem (P4-1) can be divided into M m=1γ n,m T terms, where the i th term is [q n (t)] τ p n (τ ) and represents the generated latency penalties after the (i − 1) th multicasting and before the i th multicasting of the n th message. Moreover, according to (5), the values of q n (t) in any two of M m=1γ n,m T terms are independent. Thus, it can be checked that the optimal policy for problem (P4-1) determines the time slot for new multicasting only based on the values of q n (t) since the last multicasting. Therefore, problem (P4-1) can be realigned as where (29) is derived by (27); t * is the time to start new multicasting given that the last multicasting happened at the 0 th time slot, and determined by policy π.
Exhausted search algorithm can be deployed to solve problem (P4-2). Particularly, in the case with p n (τ ) = 1, i.e., the case where the latency penalty explicitly represents the latency, we defineq n (t) M * τ =1 [q n (t)] τ and then problem (P4-2) can be rewritten as where its optimal policy is of a threshold type and starts new multicasting of the n th message at the t th time slot only if the value ofq n (t) exceeds the thresholdq n ; t * min{t ∈ Z + |q n (t) >q n } is the exact time slot to start new multicasting; (30) is derived by (5). Obviously, we can derive the optimal value ofq n , denoted asq * n , by directly solving the equality (29), and then derive the optimal value of problem (P4-3) as −E {qn(t)} min{t∈Z + |qn(t)>q * n } t=1q n (t). Therefore, in the case with p n (τ ) = 1, the upper bound of problem (P4) is derived as Therefore, problem (P5) can be released as where γ(T ) is defined as [γ(T )] (n,m) γ n,m (T ). Apparently, problem (P5-1) can be solved by integer programming [27], and the derived optimal value of problem (P5-1) is an upper bound of the objective function of problem (P2).

V. SIMULATION AND NUMERICAL RESULTS
This section evaluates the performance of the proposed PW-DDPG and compares it with the existing multicast scheduling algorithms. Specifically, we consider that the latency penalty explicitly represents the latency, i.e., p n (τ ) = 1 holds. The channel coefficient between the BS and any MU, i.e., h n,m,k (t), is modeled as an i.i.d. random variable over time, which equals 0.5 with probability of 0.8 and 1 with probability of 0.2. In the proposed PW-DDPG, the actor and critic NNs have two hidden layers and each layer has 60 nodes, and the sizes of the experience buffer and the minibatch are set as 5000 and 400, respectively. We investigate the following benchmark algorithms: a) optimal stopping method [17], which is suboptimal for the scenario with single cached message and single available channel; b) RVI [18], which is optimal for the scenario with multiple messages and single channel and with only a few MUs; c) benchmark by solving problem (P5-1) (benchmark (P5-1) for short); d) random policy, which randomly picks one feasible action at each time slot; e) greedy policy, which picks the feasible action with the largest instant reward; f) proportion pseudo greedy policy, which first randomly selects a certain proportion of actions and then picks the feasible one with the largest instant reward. Remarkably, the RVI in b) cannot address the time-varying issue and thus can only be applied in the cases with T n,1 = 1 for all n ∈ {1, 2, · · · , N }. Fig. 4 investigates the performances of various algorithms under the scenario where 20 MUs send requests for downloading one common message with probability of 0.5 and the BS serves these MUs over one multicasting channel, i.e., E 1 = 20, N = 1, M = 1, and α 1 = 0.5. Moreover, the multicasting of the message consumes 2 time slots, i.e., T 1,1 = 2. From Fig. 4(a), it is observed that the optimal stopping method proposed in [17] achieves the greatest average reward of -168.92 and is close to -168.88 given by the benchmark (P5-1). The proposed PW-DDPG converges after around 20000 time slots,  where the first 5000 time slots are utilized to purely fill up the experience buffer and thus show no performance gain, and achieves the average reward of -190.09 with the exploration rate as 0.1, and -183.48 with as 0.05, which are very close to the benchmark (P5-1). Fig. 4(b) shows that PW-DDPG achieves a similar tradeoff between the energy consumption and the latency penalty compared with the benchmark (P5-1). Fig. 5 again investigates the single-message and singlechannel scenario while with varied tradeoff parameter V , and illustrates the tradeoff curves between the average latency and the average energy consumption of different algorithms. We simulate the cases with T 1,1 = 2 and T 1,1 = 4, where each multicasting consumes 2 and 4 time slots, respectively. And in both cases, it is observed that the tradeoff curve of the optimal stopping method coincides with the benchmark (P5-1), which indicates that the optimal stopping method achieves the optimal tradeoff. And the tradeoff curve of PW-DDPG is very close to that of the benchmark (P5-1) and the gap between them can be further shrunk by decreasing the exploration rate from 0.1 to 0.05. Fig. 6 studies the scenario with multiple messages and single channel and compares the performances of PW-DDPG with RVI. Due to the high computational complexity and multiple model-level requirements on applying RVI [18], we further specify the scenario as: a) 16 MUs send requests for downloading two common messages; b) the maximum accumulated numbers of two messages' requests are limited to be N 1 = 10 and N 2 = 10 and the exceeded requests will be abandoned; c) E 1 = 8, E 2 = 8, α 1 = 0.2, α 2 = 0.6, and T 1,1 = T 2,1 = 1. It is observed from Fig. 6(a) that RVI achieves the average latency penalty 2.80, which exactly equals the benchmark (P5-1). PW-DDPG achieves a very close value 2.90, while the greedy policy achieves 3.03. Fig. 6(b) shows that RVI has a threshold structure, i.e., the regions to multicast two messages are separated by a switch curve. By executing the online validating for 600000 time slots, PW-  DDPG is verified to have an approximate threshold structure, i.e., it has a very high frequencies (above 97.5%) to multicast message 1 in the green region and message 2 in the blue region, and these two regions are separated by almost the same switch curve with that of RVI. Notably, RVI has a very high computational complexity, i.e., it is derived by first evaluating

N1+N2+2 N1+1
policies and then selecting the optimal one, where the evaluation of each policy requires one individual long-term simulation. Fig. 7 again studies the scenario with two messages and single channel, while considers E 1 = 20, E 2 = 20, N 1 = 22, and N 1 = 20. Apparently, it is computationally difficult to find the optimal policy by RVI, which requires to evaluate 44 21 = 2.0126 × 10 12 policies in advance. Fig. 7(a) shows that PW-DDPG achieves average latency penalty 7.60 and is close to the benchmark (P5-1) 6.92. Fig. 7(b) shows that PW-DDPG remains an approximate threshold structure. We also consider that the multicastings for four messages consume different time durations and the latency penalty functions of four messages are different, neither. To be specific, (T 1,m , T 2,m , T 3,m , T 4,m ) = (1, 2, 3, 4) for all m ∈ {1, · · · , 6} and (p 1 (τ ), p 2 (τ ), p 3 (τ ), p 4 (τ )) = (1, 2, 2, 2). Apparently, the MDP under this scenario has a large discrete action space of magnitude |A|= 5 6 = 15625 and multiple time-varying constraints, where existing algorithms cannot solve it efficiently. The proposed PW-DDPG again deploys two hidden layers in the actor and critic NNs, while the two hidden layers now contain 350 and 300 nodes, respectively. It is observed that PW-DDPG converges after around 75000 Thus, PW-DDPG has a much lower computational complexity. In Fig. 8(b), PW-DDPG again achieves a better tradeoff than the other two algorithms.

VI. CONCLUSIONS
We consider the multicast scheduling problem for multiple messages over multiple channels, which jointly optimizes the energy efficiency and the average latency penalty. This problem is formulated as an infinite-horizon MDP, which is challenging due to the large discrete action space and multiple time-varying constraints. Specifically, the MDPs with the former feature cannot be efficiently addressed by existing algorithms and the latter feature induces dynamic action space, which further complicates the problem. To simplify the formulated MDP, we first analyze the intrinsic features of this problem under stationary policies and succeed to reduce the state space by refining the reward function. To address the discrete action space issue, we propose the permutation-based action embedding method, which embeds action in the "permutation" manner and overcomes the drawbacks of the conventional action embedding. To address the time-varying constraints issue, we design PAE and FAE modules to separately restrict the exploitation and exploration modes to selecting feasible actions. Remarkably, the proposed permutation-based action embedding can also be combined with other DRL algorithms, such as asynchronous advantage actor-critic (A3C) [28] and PPO, to better fit the structure of the considered MDPs, and we will explore these combinations in future works.

APPENDIX A SKETCH PROOF OF PROPOSITION 2.1
To show that problems (P1-1) and (P2) are equivalent, we first prove that for any stationary policy of problem (P1-1), there exists another stationary policy of problem (P2) such that their objective functions are equal. Then, we prove that for any stationary policy of problem (P2), there also exists another stationary policy of problem (P1-1) such that their objective functions are equal. Based on the above two results, the optimal values of problems (P1-1) and (P2) are obviously equal.

T t=1
Is(s(t)) T , G and Q are components of s = (G, C, Q), and a m is the m th entry of a; equality (ii) holds since π is a stationary policy; in (33),ŝ = (G,ĉ, Q) and N (ŝ) is defined as {s ∈ S| µ π s π(s, a), whereπ(ŝ, a) and π(s, a) are the probabilities of selecting action a at stateŝ and s, respectively. It can be checked that the occurrence probability of the stateŝ under policȳ π in problem (P2) equals the summation of the occurrence probabilities of the states in N (ŝ) under policy π in problem (P1-1), i.e., µπ s = s∈N (ŝ) µ π s . Therefore, we derive that We use similar techniques to proof the inverse, i.e., for each stationary policyπ for problem (P2), there exists another policy π for problem (P1-1) such that (15)|π= (10)| π holds, and the proof is thus omitted.

APPENDIX B PROOF OF PROPOSITION 4.1
Based on the definitions of I n (x) and I(x) in Section II-3, we obtain I n (a m (t)) = T t=1 I (a m (t)), which indicates that there are in total N n=1 γ n,m (T ) times of multicastings over the m th channel during T time slots. We start these multicatings at the t 1 th , t 2 th , · · ·, t N n=1 γn,m(T ) th time slot and transmit the n 1 th , n 2 th , · · ·, n N n=1 γn,m(T ) th message, respectively, i.e., a m (t i ) = n i . Obviously, the number of the total consumed time slots by these multicastings is N n=1 γn,m(T ) i=1 T ni,m and it equals N n=1 γ n,m (T )T n,m . Without loss of generality, consider 1 ≤ t 1 < t 2 < · · · < t N n=1 γn,m(T ) ≤ t N n=1 γn,m(T ) + T n N n=1 γn,m (T ) ,m − 1 ≤ T , where the last inequality is to cover the last multicasting within T time slots. Now, we prove (24). First, it followŝ c m (t i +j i ) > 0, i ∈ {1,· · ·, N n=1 γ n,m (T )}, j i ∈ {1,· · ·,T ni,m −1}, which is straightforward from (12) and indicates that the m th channel is occupied from the t i th to (t i + T ni,m − 1) th time slot by the multicasting of the n i th message. Then, based on (12) and (35), it follows t i+1 ≥ T ni,m + t i , i ∈ {1, 2, · · · , N n=1 γ n,m (T ) − 1}, (36) which indicates that the (i + 1) th multicasting over the m th channel should start from or after the (T ni,m + t i ) th time slot. Finally, we sum up both sides of (36) over i ∈ {1, 2, · · · , Then, the proof is complete.