\section{Introduction}
\label{sec:introduction}
Emerging ICPSs, such as smart grid, smart manufacturing, and smart transportation are spatially distributed and high-dimensional. These systems require high reliability, communication between numerous devices, low latency, power-efficient communication, and high computational load \cite{tranoris2018enabling, zerihun2020effect}. To manage these requirements, 5G and beyond 5G networks present a wide range of services that are classified as 1) enhanced mobile broadband (eMBB), 2) ultra-reliable and low-latency communication (URLLC), and 3) massive machine-type communication (mMTC).
The eMBB, URLLC, and mMTC services provide a high data rate with a moderate latency, a communication with low end-to-end delay, and connecting many devices respectively. Because the type of most communications in control sub-systems is URLLC, 5G network is a good choice for exchanging data between a controller and sensors or actuators.\\
However, there are some serious challenges to deploying 5G network in ICPSs. Specifically, ICPSs with 5G networks have limited network resources and lack the desired stability and performance guarantees \cite{eisen2019control,liu2020latency}. The performance of the control sub-system is defined as achieving the required dynamics response, which is specified by measures of performance such as a desirable steady-state tracking error.
The stability and performance of a control sub-system may be guaranteed through periodic transmissions with a high data rate. It, however, comes at the cost of a higher packet loss rate due to limited resources in wireless networks \cite{lu2015real}. Furthermore, in many applications of ICPSs, most wireless devices rely on batteries and their battery life may be significantly reduced by the increased transmission rate \cite{funk2021learning}. Consequently, the event-triggered control (ETC) method is proposed, in which the transmission times of the control sub-system are triggered based on a predefined event instead of a periodic transmission. This event is characterized according to the stability and performance requirement of the control sub-system.\\ In recent years, extensive research has concentrated on different classes of ETC strategies; see \cite{ge2019distributed, asurvey} and the references therein. Besides, in this context, there are substantial works, which prove better energy-saving and performance of ETC in comparison with the traditional periodic control \cite{antunes2014rollout, dolk2017event}. Nevertheless, these works analyze only low-dimensional or linear models of control sub-systems \cite{funk2021learning, xu2020remote}. Moreover, the analysis of the event-triggered control becomes too complicated when the volatile properties of wireless communication such as delay, limited resources, packet drops, and unreliable links are considered.\\
The design of the event-triggered control in the presence of unreliable links and packet losses has been recently drawn a lot of attention \cite{zhu2020event, wu2020event, xu2020remote}. However, in addition to packet drops, there exist many other features of wireless communication, such as the delay and limited resources, which make a direct impact on the stability and performance of control sub-systems. To deal with these interconnections between the control and communication sub-systems, the joint design method is taken in ICPSs \cite{xu2020reinforcement, liu2019controller, funk2021learning,mamduhi2020cross}. However, developing an analytical model of all control and network features is a fundamental challenge to this method. This is because the sub-systems are typically high-dimensional and the conditions of radio resources are continuously and randomly changing.\\
Therefore, researchers have used model-free reinforcement learning (RL) in the joint design of ICPS' sub-systems \cite{demirel2018deepcas, vamvoudakis2018model,funk2021learning, xu2020reinforcement,leong2020deep}. In \cite{demirel2018deepcas}, RL is used for proposing a sensors scheduler while the controller is designed beforehand. The actor-critic RL method is also used in \cite{vamvoudakis2018model} to learn the event-triggered control. In \cite{funk2021learning}, option method of Deep RL (DRL) is used for joint optimization of an event policy and a control policy. The event policy determines when the control input should be transmitted and the control policy determines what the control input value should be sent. Nonetheless, the varying characteristics of the wireless network are not considered in \cite{funk2021learning, vamvoudakis2018model}.
In \cite{xu2020reinforcement}, RL approach is used to jointly design the sampling rate of the control sub-system and the modulation type of the wireless network.\\ Although stability is an essential property for every control sub-system, RL methods could hardly guarantee the stability and reliability of a learning-based controller \cite{han2020actor}. Nonetheless, in \cite{han2020actor, han2021reinforcement,zhang2021safe}, a learning-based controller with uniformly ultimate boundedness (UUB) stability guarantee is proposed, which can be usefully employed in ICPSs with safety constraints. In general, UUB stability says that if the norm of starting state variables of a control sub-system is less than a specified value, then the state variables will eventually enter the neighborhood of the sub-system's equilibrium within a finite time and will never escape from this neighborhood set afterwards \cite{han2021reinforcement}.\\
The goal of this paper is to jointly design the event-triggered control and the energy-efficient allocation of radio resources in an ICPS. To the best of our knowledge, this joint design problem has not yet been studied. We propose to use a novel Hierarchical RL (HRL) approach with UUB stability guarantee to solve the problem.
Our contributions are as follows.
 \begin{itemize}
    \item We assume an ICPS containing multiple eMBB users and a control plant with multiple URLLC users sharing a single cell Orthogonal Frequency-Division Multiple Access (OFDMA) network. We formulate the joint design of the event-triggered control and the energy-efficient resource allocation in the ICPS as a multi-objective optimization problem. The goals of the problem are both minimizing the number of updates on the actuators' input and the energy consumption in the downlink. The constraints of this problem contain the dynamics and UUB stability of the control plant, the minimum Quality of Service (QoS) demand of eMBB and URLLC users, and the power and sub-carrier constraints of the OFDMA network.
   \item The problem is high-dimensional, complicated and associated with a hybrid action space. To handle these properties, we combine Cascade Attribute Learning Network (CAN) method and option-critic method to develop a novel model-free HRL approach with
   UUB stability guarantee. First, we use CAN method and decouple the problem into two low-dimensional sub-problems of control and resource allocation. We show that using the decoupling method leads to a Pareto solution to the optimization problem. In the second step, we use option-critic method, which is reformulated as Double Actor-Critic (DAC) architecture, to address each sub-problem with a hybrid action space.
   \item The novel model-free HRL with
   	UUB stability guarantee can simultaneously learn four policies: 1) update time policy on the actuators' input, 2) control policy, which determines the value of control input, 3) energy-efficient sub-carrier allocation policy, and 4) energy-efficient power allocation policy.
   \item We demonstrate the effectiveness and capability of the proposed approach by several simulation results. In comparison with a disjoint and model-based method, our numerical simulation results show that both the number of updates on the actuators' input and the downlink energy consumption are reduced significantly by applying the proposed approach. Moreover, we show the capability of the proposed approach compared with the soft actor-critic algorithm.
\end{itemize}
\indent This paper is outlined as follows. The system model and problem formulation is described in Section II. The proposed approach is presented in Section III. In Section IV, simulation results are discussed. In Section V, the paper's conclusion and future work are given.
\begin{figure}[t]

	\centerline{\includegraphics[width=0.5\textwidth]{Fig1.jpg}}
	\label{fig:System model}
	\caption{System model of the considered ICPS}
\end{figure}
\section{SYSTEM MODEL and PROBLEM FORMULATION}
\subsection{System Model}
 Consider a model of ICPS that consists of 1) a control plant, 2) a central event-triggered controller, and 3) a downlink model of an OFDMA cellular network (Fig. 1). The state values of the control plant, are measured by multiple sensors and sent to the event-triggered learning-based controller. Next, the controller calculates and sends the control input to actuators, whenever required, through the OFDMA single cell network.
 Following the 5G architecture explained by International Telecommunication Union (ITU), the central learning-based controller is supposed to run on a specific or shared hardware in the central office data center layer, which is placed near the network's Base Station (BS) \cite{HUAWEI}.\\
 \textbf{Control Plant}: We suppose dynamics of the control plant is unknown, that is:
 \begin{equation}
 \begin{split}
 \ x[k+1] & = f(x[k], u[k], \omega[k]),\\
   y[k] &= v(x[k]),
 \label{eqn:pl1}
 \end{split}
 \end{equation}
where $f(.)$ and $v(.)$ are unknown functions, $x[k]\in \mathbb{R}^{n}$, $u[k]\in \mathbb{R}^{m}$, and $y[k]\in \mathbb{R}^{q}$ denote the vector of the control state, control input, and sensors' output at discrete time $k\in\mathcal{K}$ ($\mathcal{K}= \{1,2,...,K\}$) respectively. Also, vector $\omega[k]\in \mathbb{R}^{m}$ is actuation disturbances at discrete time $k\in\mathcal{K}$.
We assume the control plant, described by dynamics (\ref{eqn:pl1}), is completely state observable, as it is regularly assumed in the related literature, e.g. \cite{radac2020robust}.\\
\textbf{Event-triggered Controller}: When the sensors' output vector ($y[k]$) is received by the central event-triggered controller, it decides whether actuators' input should update ($\beta[k]= 1$) or ignore the update and save wireless resource ($\beta[k]= 0$). This decision has been taken based on UUB stability guarantee of the control plant, defined in what follows.
\begin{def1} \label{D1}
	A control plant is uniformly ultimately bounded with ultimate bound $\rho$, if there are positive constants $b, \rho$ and $\forall \zeta < b:\exists T(\zeta,\rho)$, such that $||x[k_0]||<\zeta	\Rightarrow||x[k]||<\rho, \forall k \geq k_0 + T$. If $\zeta$ can be arbitrary large, then the control plant is globally uniformly ultimately bounded.
\end{def1}
\noindent In addition, if update variable $\beta[k]= 1$, then the controller calculates the control input variable ($u[k]$) considering UUB stability. We assume that Zero Order Hold (ZOH) holds actuators' input constant between two consecutive updates. This can be mathematically given by:
\begin{equation}
\ u[k]  = u[k] \beta[k] + u[k-1](1- \beta[k]): \beta[k] \in \{0,1\}.
\label{eqn:event-control}
\end{equation}
\textbf{OFDMA Network}:
 We assume the downlink model of a single cell OFDMA network with one BS. The model has $N$ downlink users denoted by $\mathcal{N}=\{1,2,...,N\}$. The downlink users have a set of $N^c$ control plant users (URLLC users) defined by $\mathcal{N}^c=\{1,2,...,N^c\}$ and a set of $N^e$ eMBB moving users (coexisted with the control plant users) defined by $\mathcal{N}^e=\{1,2,...,N^e\}$. It is noted that URLLC users and control plant users are employed interchangeability from hereon. We consider that URLLC users are fixed and eMBB users move within the range of the BS coverage area. Let dividing the total bandwidth of the network in $J$ sub-carriers forming set $\mathcal{J}=\{1,2,...,J\}$. Also, let $p_{n,j}[k]$ be the base station's transmit power for communicating with downlink user $n$ on sub-carrier $j \in \mathcal{J}$ at discrete time $k$. The variable of $p_{n,j}[k]$ is assumed continues. The overall power transmit of the BS is limited to a maximum value represented by $\overline{P}_{BS}$, which means $\sum_{n=1}^{N}\sum_{j=1}^{J}p_{n,j}[k]\leq\overline{P}_{BS}$. Moreover, the BS' total power usage in the considered ICPS is calculated as \cite{loodaricheh2014energy}:\\
 \begin{equation}

 P^\text{BS}_\text{total}[k]=P_\text{cst}+\epsilon^\text{BS}\sum\nolimits_{n \in \mathcal{N}} \sum\nolimits_{j \in \mathcal{J}}a_{n,j}[k]p_{n,j}[k],

 \label{eqn:poweD}
 \end{equation}
where $P_{\text{cst}}$ is a constant power used by BS circuit, $\epsilon^\text{{BS}}$ is the amplifier inefficiency constant, $a_{n,j}[k]$ is the sub-carrier allocation variable, which is a binary variable. $a_{n,j}[k]=1$ if sub-carrier $j$ is allocated to downlink user $n$ at discrete time $k$, or else, $a_{n,j}[k]=0$. Also, $\textbf{A}[k]$ and $\textbf{P}[k]$ are power and sub-carrier allocation matrices at discrete time $k$ respectively
($\textbf{A}[k]{\mathop:}=[{a_{n,j}}_{(n\in \mathcal{N},j\in \mathcal{J})}[k]]$ and $\textbf{P}[k]{\mathop:}=[{p_{n,j}}_{(n\in \mathcal{N},j\in \mathcal{J})}[k]]$).\\
The downlink Signal-to-Noise Ratio (SNR) for user $n$ on sub-carrier $j$ is given by \cite{tang2019service}:\\
\begin{equation}
\gamma_{n,j}[k]=\frac{p_{n,j}[k]g_{n,j}[k]}{N_0[k]},
\label{eqn:gammau}
\end{equation}
where $g_{n,j}[k]$ is the channel gain for each user $n$ on sub-carrier $j$ at discrete time $k$ and $N_0[k]$ denotes the corresponding additive white Gaussian noise power at the receiver of user $n$.
In accordance with the Shannon's formula, the achievable instantaneous transmission rate for each eMBB user $n \in \mathcal{N}^e$ is computed in bit/s as:\\
\begin{equation}
R^e_n[k]=\sum\nolimits_{j \in \mathcal{J}}w a_{n,j}[k]\log(1+\gamma_{n,j}[k]),
\label{eqn:rateu}
\end{equation}
where $w$ is the bandwidth of sub-carrier $j$. Moreover, the QoS requirement for each eMBB user $n \in \mathcal{N}^e$ is computed in terms of a minimum transmission rate \cite{tang2019service}. Therefore, the required QoS of eMBB users is represented by:
\begin{equation}
	R^e_n[k] \geq\overline{R}^e_{n}[k],\forall n\in \mathcal {N}^e,
	\label{eqn:ratecons}
\end{equation}
where $\overline{R}^e_{n}[k]$ is the minimum required QoS of eMBB user $n$ at discrete time $k$.
The packet size of URLLC users are generally short so the Shannon's formula cannot exactly describe their transmission rate \cite{polyanskiy2010channel,tang2019service}. The achievable transmission rate of URLLC users with the finite blocklength channel coding method is derived in \cite{polyanskiy2010channel} as:
\begin{equation}
\begin{split}
R^c_n[k]\!=&w\!\sum\nolimits_{j \in \mathcal{J}}\!a_{n,j}[k](\log(1+\gamma_{n,j}[k])\\&-\! \sqrt{\!\dfrac{V_{n,j}[k]}{C_{n,j}}}Q^{-1} (\epsilon)\log e),
\label{eqn:rated}
\end{split}
\end{equation}
where $C_{n,j}$ is the number of symbols in each codeword block, $Q^{-1}$ is the inverse of Gaussian Q-function, $\epsilon$ is the
error probability, and $V_{n,j}$ is dispersion of sub-carrier $j$ for user $n \in \mathcal{N} ^c$ given by:
\begin{equation}
V_{n,j}[k]=1-\dfrac{1}{(1+\gamma_{n,j}[k])^2}.
\label{eqn:V}
\end{equation}
In a single time slot $k$, to satisfy the required QoS of URLLC users, it is necessary to provide the achievable instantaneous data rate condition as below:
\begin{equation}
R^c_n [k]\geq \dfrac{L_c}{T_c[k]},\forall n\in \mathcal {N}^c,
\label{eqn:ratecon_urllc}
\end{equation}
where $L_c$ is the length of actuator's packet size in bits and $T_c[k]$ is the maximum tolerable transmission delay for the packet. We calculate $T_c[k]$ according to the given maximum tolerable end-to-end (e2e) delay between the controller and actuators.
Let $T^\text{{comp}}_\text{{max}}$ be the maximum queuing and computation delay that is $T^\text{{comp}}[k]\leq T^\text{{comp}}_\text{{max}}$ and the propagation delay is negligible. Thus, we conservatively assume the e2e delay is:
 \begin{equation}

 T_\text{e2e}[k]=T_c[k]+T^\text{{comp}}_\text{{max}}.

 \label{eqn:e2edelay}
 \end{equation}
Noticeably, we assume the minimum reliability requirement for URLLC users is satisfied through some enabler techniques such as low-rate codes.
\vspace{-5mm}
\subsection{Problem formulation}
We now formally state the joint design problem of the event-triggered control and the energy-efficient resource allocation of the OFDMA network, as a multi-objective optimization problem. It aims to minimize both the number of updates on the actuators' input and the total downlink power usage, subject to the dynamics and UUB stability of the control plant, the QoS demands of eMBB and URLLC users, power and sub-carrier constraint, and the maximum practicable level of the BS' transmit power. This problem is formulated as:\\
\begin{equation}
\begin{aligned}
\begin{split}
&\underset{\{\beta[k]\},\{u[k]\}}{\text{minimize}} \sum_{i=1}^K\beta[i] \quad \\
&\underset{\{\textbf{A}[k]\},\{\textbf{P}[k]\}}{\text{minimize }} \sum_{i=1}^K P^\text{BS}_\text{total}[i]\\
&\text{subject to}:\\
&C_{1}:x[k+1] = f(x[k], u[k], \omega[k])\\
&\quad \quad y[k] = v(x[k]):\forall k \in \mathcal{K}\\
&C_{2}:u[k] = u[k] \beta[k] + u[k-1](1- \beta[k]):\forall k \in \mathcal{K}\\
&C_{3}:||x[0]||<\zeta \Rightarrow ||x[k]||<\rho: \forall k \geq T(\zeta,\rho)\\
&C_{4}: \beta[k] \in\{0,1\} :\forall k \in \mathcal{K}\\
&C_{5}:R^e_n[k]\geq\overline{R}^e_{n}[k] :\forall k \in \mathcal{K}, \forall n\in \mathcal {N}^e\\
&C_{6}:R^c_n[k] \geq \dfrac{L_c}{T_c[k]}\beta[k] :\forall k\in \mathcal{K}, \forall n\in \mathcal {N}^c\\
&C_{7}:\sum_{n=1}^{N}a_{n,j}[k]\leq 1 :\forall k\in \mathcal{K}, \forall j\in \mathcal {J}\\
&C_{8}:a_{n,j}[k]\in\{0,1\} :\forall n\in \mathcal {N},j\in \mathcal {J},k\in \mathcal{K}\\
&C_{9}:\sum_{n=1}^\mathcal{N}\sum_{j=1}^\mathcal{J}p_{n,j}[k]a_{n,j}[k]\leq\overline{P}_{BS} :\forall k\in \mathcal{K},
\label{eqn:optt}
\end{split}
\end{aligned}
\end{equation}
where constraints $C_{1}$, $C_{2}$, and $C_{3}$ illustrate the plant dynamics, the event-triggered controller function, and UUB stability requirement of the control plant respectively. Constraint $C_{4}$ shows update variable $\beta[k]$ takes binary value. $C_{5}$ and $C_{6}$ represent the required QoS of eMBB and URLLC users respectively.
Constraints $C_{7}$ and $C_{8}$ are related to the exclusive assignment of the sub-carrier in the OFDMA network. And constraint $C_{9}$ shows the maximum allowable transmit power of the BS. \\
In multi-objective optimization problem (\ref{eqn:optt}), thanks to minimizing the second objective, the transmit power of the control plant's users is reduced. Consequently, the downlink transmission rates are reduced and the transmission delay is increased. Accordingly, to guarantee UUB stability of the control plant  ($C{3}$), the number of updates on the actuators' input is increased in future time steps and the first objective function is increased. Due to the trade-off between these two objective functions, the idea of the Pareto optimality is employed as a solution for problem (\ref{eqn:optt}) \cite{boyd2004convex}. The Pareto optimal solution is defined as follows.
\begin{def2}\label{D2}
	Assuming a multi-objective optimization problem with $f_{i}(k), i \in \{1,2,...,I\}$, as its objective functions and considering all objectives are minimizing functions, a feasible solution, $k^*$, can dominate another one, $k^{**}$, (or $k^*$ is better than $k^{**}$) if:
	\begin{enumerate}
		\item $f_{i}(k^*)\leq f_{i}(k^{**})$ for all $i\in \{1,2,...,I\}$ and
		\item $f_g(k^*)< f_{g}(k^{**})$ for at least one $g\in \{1,2,...,I\}$.
	\end{enumerate}
	$k^*$ is named as a Pareto optimal solution when any other solution cannot be found to dominate $k^*$. In other words, $k^*$ is a Pareto optimal solution if and only if it is a feasible solution and there exists no better feasible solution.
\end{def2}
\vspace{-10pt}
\section{The Proposed approach}
In optimization problem (\ref{eqn:optt}), the dynamics model of the control plant and its interconnection with the network is unknown. To address this problem, we propose a novel model-free HRL approach. Specifically, a Markov Decision Process (MDP) is first constructed associated with problem (\ref{eqn:optt}). Due to the state and action spaces of the MDP are large, we first apply CAN method and decompose problem (\ref{eqn:optt}) into two sub-problems. Then, DAC architecture is used for solving each sub-problem with a hybrid action space.
\vspace{-10pt}
\subsection{RL-related Definition}
The joint design problem can be described by MDP $\mathcal{M}= (\mathcal{S}; \mathcal{A}; \mathcal{R}; \mathcal{P}_0; \mathcal{P}_{ss^{'}})$, where $\mathcal{S}$ is the set of possible states, $\mathcal{A}$ is the set of actions, $\mathcal{R}$ is a reward function ($\mathcal{S} \times \mathcal{A} \rightarrow \mathbb{R}$), $\mathcal{P}_0$ is an initial distribution ($\mathcal{S} \rightarrow [0,1]$), $\mathcal{P}_{ss^{'}}$ is the probability of states transition ($\mathcal{S} \times \mathcal{A}\times \mathcal{S} \rightarrow [0,1]$). The state at time step $k$, $\mathcal{S}[k]$, is defined as:
\begin{equation}
\mathcal{S}[k] = \{s^{c}[k], s^{N}_n[k] \},
\label{eqn:state}
\end{equation}
where $s^{c}[k] = y[k]$, $s^{N}_n[k] \in\{0,1\}: n \in\mathcal{N}$ denotes status of URLLC and eMBB users at environment time step $k$, which $s^{N}_n[k] = 1$ if user $n$ receives its minimum required rate; otherwise $s^{N}_n[k] =0$. We consider the learning agent action at time step $k$, $\mathcal{A}[k]$, as follows:
\begin{equation}
\mathcal{A}[k] = \{\beta[k], u[k], \textbf{A}[k], \textbf{P}[k]\}.
\label{eqn:action}
\end{equation}
An action is taken, at each time step $k$, on the basis of policy $\pi(\mathcal{A}[k]| \mathcal{S}[k])$, which is a likelihood function of each action for every possible state. By choosing $\mathcal{A}[k]$, the environment state is transmitted from current state $\mathcal{S}[k]$ to $\mathcal{S}[k+1]$ according to the probability of $\mathcal{P}(\mathcal{S}[k+1] | \mathcal{S}[k], \mathcal{A}[k])$ and also a reward of $\mathcal{R}[k+1]$ is gotten ($\mathbb{E}(\mathcal{R}[k+1])= \mathcal{R}(\mathcal{S}[k],\mathcal{A}[k])$). Assuming the transition trajectory as $\tau = (\mathcal{S}[0], \mathcal{A}[0], ..., \mathcal{S}[K])$, the goal of RL is to obtain a policy (${\pi^0}$), which maximize the expected receiving cumulative reward trough the trajectory, which is given by $\mathcal{R}(\tau) = \mathbb{E}(\sum_{j=k}^{K}\gamma^{j-k}\mathcal{R}[k])$
where $0 \leq \gamma  < 1$ denotes the discount factor showing the important weight of future rewards. $\mathcal{R}(\tau)$ is the cumulative reward of an episode between the step of $k$ and the terminal step of $K$.
\vspace{-10pt}
\subsection{Applying CAN Method and Decomposing the Problem}
It is obvious that the size of the state and action spaces of the joint design problem may be too large in practical cases. In such a high-dimensional and complex problem, the speed of learning is considerably reduced. Furthermore, the training process generally consumes an unreasonable amount of computation power in the high-dimensional problem. To manage these challenges, CAN method is used as explained in \cite{chang2020cascade}. In CAN method, the learning process of a complicated problem is decomposed into low-dimensional attribute modules, which are linked in cascade series. The state space of every attribute is determined as minimum as possible provided that the space can completely describe the attribute, indicated by $\mathcal{S}=\{\mathcal{S}^0, \mathcal{S}^1, \mathcal{S}^2,...\}$. Also, every attribute enjoys its own reward function ($\mathcal{R}=\{\mathcal{R}^0, \mathcal{R}^1, \mathcal{R}^2,...\}$).
Moreover, the transition probability distribution in every attribute is indicated by $\mathcal{P}=\{\mathcal{P}^0, \mathcal{P}^1, \mathcal{P}^2,...\}$.
Although it is shown that CAN method makes the training process significantly faster and more simple in \cite{chang2020cascade}, it is not mathematically proven that applying the decoupling method results in an optimal/sub-optimal solution. Here, however, we demonstrate this through the following lemmas in the case of problem (\ref{eqn:optt}).
\begin{lemma}\label{L1}
The second objective of problem (\ref{eqn:optt}) is decreasing with respect to $\beta[k]$.
\end{lemma}
\begin{IEEEproof}
By decreasing $\beta[k]$, the number of control users that require to communicate  decreases, through which the number of downlink users, $\textit{N}$, is decreasing. Consequently, the total power consumption of the BS will decrease in accordance with equation (\ref{eqn:poweD}).
\end{IEEEproof}
\begin{lemma}\label{L2}
Assuming $\mathcal{H}=\{\beta^*,u,\textbf{A},\textbf{P}\}$ be a subset of feasible solution set of problem (\ref{eqn:optt}), if $\beta^*\in \mathcal{H}$ minimize the first objective function, then $\mathcal{H}$ is a subset of the Pareto solution set of optimization problem (\ref{eqn:optt}).
\end{lemma}
\begin{figure}[t]

	\centerline{\includegraphics[width=0.5\textwidth]{Fig2}}
	\label{fig:Proposed approach}
	\caption{Proposed approach}
\end{figure}
\begin{IEEEproof}
Lemma \ref{L2} will be demonstrated by the contradiction. Assuming that there is a feasible solution $k^{**}$, which dominant $k^*\in\mathcal{H}=\{\beta^*,u,\textbf{A},\textbf{P}\}$ ($\beta^*$ minimize the first objective function). But,
in accordance with Lemma \ref{L1}, the second objective function is decreasing by decreasing $\beta$ and also the first objective is optimized in $\beta^*$. Thus, the conditions presented in Definition 2 are not fulfilled. Accordingly, $k^{**}$ does not dominate $k^*$. As a result, the initial presumption that there is a feasible point, which dominate $k^{*}$ is contradicted.
\end{IEEEproof}
Lemma \ref{L2} allows us to decouple optimization problem (\ref{eqn:optt}) into two sub-problems as:
\begin{equation}
\begin{aligned}
&\underset{\{\beta[k]\},\{u[k]\}}{\text{minimize}} \sum_{i=1}^K\beta[i] \\
&\text{subject to}:
C_{1}, C_{2}, C_{3}, C_{4},\quad \quad \quad
\label{eqn:sub1}
\end{aligned}
\end{equation}
\text{and}
\begin{equation}
\begin{aligned}
&\underset{\{\textbf{A}[k]\},\{\textbf{P}[k]\}}{\text{minimize }} \sum_{i=1}^K P^\text{BS}_\text{total}[i]\\
&\text{subject to}: C_{4}, C_{5}, C_{6}, C_{7}, C_{8}, C_{9}.
\label{eqn:sub2}
\end{aligned}
\end{equation}
The architecture of the proposed approach applying CAN method is shown in Fig. 2. The training process of the proposed approach has two parts. In the first part, DRL policy of the base attribute module is trained to address sub-problem (\ref{eqn:sub1}). The base module is fed with $\mathcal{S}^0\subset \mathcal{S}$ and output $\mathcal{A}^0 \subset \mathcal{A}$, considering reward function $\mathcal{R}^0$. Notably, $\mathcal{A}^0$ contains $\beta[k]$ and $u[k]$, which $u[k]$ is a continues variable and $\beta[k]$ is a binary variable. Having decided $\beta[k]$, DRL policy of the first attribute module is trained subsequently, which is accountable to solve sub-problem (\ref{eqn:sub2}). This module is fed with $\mathcal{S}^1 \subset \mathcal{S}$ and output power matrix $\textbf{P}[k]$ along with sub-carrier matrix $\textbf{A}[k]$, considering reward function $\mathcal{R}^1$.\\
The action space of each sub-problem is a hybrid space, and the majority of regular RL-based solutions are not appropriate to solve these hybrid problems \cite{funk2021learning}. Therefore, to address each sub-problem, we propose to use option-critic method, which is reformulated as DAC architecture in \cite{zhang2019dac}, since it is well-suited to deal with hybrid action space \cite{hou2020data, funk2021learning}.
\vspace{-3mm}
\subsection{The Base Attribute Module}
To handle sub-problem (\ref{eqn:sub1}), the state and action spaces of the base module are defined as:
\begin{equation}
\begin{split}
\mathcal{S}^0[k] &= \{s^{c}[k]\},\\
\mathcal{A}^0[k] &= \{\beta[k], u[k]\}.
\label{eqn:state_base}
\end{split}
\end{equation}
The base module is responsible for learning a policy (${\pi^0}$) over $\beta[k]$ and $u[k]$. The policy aim to maximize the expected receiving cumulative reward through transmission trajectory $\tau^0=\{\mathcal{S}^0[0],\beta[0], u[0], \mathcal{S}^0[1],\beta[1], u[1], ..., \mathcal{S}^0[K]\}$. The reward function of the base module is defined as:
	\begin{equation}

	\!\mathcal{R}^{0}(\mathcal{S}^0[k],\mathcal{A}^0[k])\!= \mathcal{R}^{ctrl}[k] -\mu_1\beta[k],
	\label{eqn:reward_base}

	\end{equation}
where the first term ($\mathcal{R}^{ctrl}[k]$) is the control reward and the second term ($-\mu_1\beta[k]$) is to minimize the number of updates on the actuators' input. The control reward is defined to encourage the control plant to reach its specified targets. Also, $\mu_1$ is a hyper-parameter denoting the penalty weight of the number of actuators updates.\\ To guarantee UUB stability of the learning controller with policy ${\pi^0}$, we use a more general definition of UUB stability presented in \cite{han2021reinforcement}. Indeed, in \cite{han2021reinforcement}, the classical definition of UUB stability (Definition 1) is extended for general cases in which the stability constraint functions are not necessarily the norm of the control state ($||x[k]||$).
Let $\mathcal{C}_{\pi^0}(\mathcal{S}^0[k])\doteq \mathbb{E}_{\mathcal{A}^0[k]\sim{\pi^0}}\mathcal{C}(\mathcal{S}^0[k],\mathcal{A}^0[k])$ be the constraint function under the policy ${\pi^0}$ and $\mathcal{C} (\mathcal{S}^0[k],\mathcal{A}^0[k])$ be a continuous nonnegative constraint function, which is defined to measure how good or bad a state$-$action pair of the base module is. The general definition of UUB stability with respect to $\mathcal{C}_{\pi^0}(.)$ is stated in what follows.
	\begin{def3}\label{D4}
		A control plant is UUB with respect to $\mathcal{C}_{\pi^0}(.)$, if there are positive constants $b, \rho$ and $\forall \zeta < b$ : $\exists T(\zeta,\rho)$, such that $\mathcal{C}_{\pi^0}(\mathcal{S}^0[k_0])<\zeta \Rightarrow\mathcal{C}_{\pi^0}(\mathcal{S}^0[K])<\rho, \forall k \geq k_0 + T$.
	\end{def3}
It is shown that Definition 3 is an inherent feature of the control plant when it is UUB stable. Thus, if the control plant is UUB with respect to $\mathcal{C}_{\pi^0}(.)$, then the closed-loop control is UUB \cite{han2021reinforcement, zhang2021safe}. It is noted that UUB points to the property defined by Definition 3 from hereon.
\begin{theorem1}\label{T1}
	Assuming that the Markov chain induced by policy $\pi^0$ is ergodic, $\Lambda()\doteq\{\mathcal{S}^0[k]\in \mathcal{S}^0 |\mathcal{C}_{\pi^0}(\mathcal{S}^0[k])\geqslant\rho\}$, and $\mathbb{I}_{\Lambda}(\mathcal{S}^0[k])=\left\{
	\begin{array}{rl}
	&1 \quad\quad \qquad \mathcal{S}^0[k]\in\Lambda\\
	&0 \quad \quad\qquad \mathcal{S}^0[k]\notin\Lambda
	\end{array}\right.$, if there are a function $\Gamma(\mathcal{S}^0[k]):\mathcal{S}^0\rightarrow \mathbb{R}^+$ and positive constants $\alpha_1, \alpha_2, \alpha_3$, and $\rho$, such that
	\begin{equation}
	\alpha_1\mathcal{C}_{\pi^0}(\mathcal{S}^0[k])\leq\Gamma(\mathcal{S}^0[k]) \leq\alpha_2\mathcal{C}_{\pi^0}(\mathcal{S}^0[k]), \forall \mathcal{S}^0[k]\in \mathcal{S}^0,
	\end{equation}
	and
	\begin{equation}
	\begin{aligned}

	&\mathbb{E}_{\mathcal{S}^0[k]\sim\Omega_{\overline{\rm K}}}(\mathbb{E}_{\mathcal{S}^0[k+1]\sim\mathcal{P}^0}\Gamma(\mathcal{S}^0[k+1])\mathbb{I}_{\Lambda}(\mathcal{S}^0[k+1])-\\
	&\Gamma(\mathcal{S}^0[k])\mathbb{I}_{\Lambda}(\mathcal{S}^0[k]))
	<-\alpha_3 \mathbb{E}_{\mathcal{S}^0[k]\sim\Omega_{\overline{\rm K}}}\mathcal{C}_{\pi^0}(\mathcal{S}^0[k])\mathbb{I}_{\Lambda}(\mathcal{S}^0[k]),

	\end{aligned}
	\label{eqn:stability}
	\end{equation}
where $\Omega_{\overline{\rm K}}$ shows the average distribution of $\mathcal{S}^0[k]$ over the finite ${\overline{\rm K}}$ time steps,
		$\Omega_{\overline{\rm K}}(\mathcal{S}^0[k])\doteq \dfrac {1}{\overline{\rm K}} \sum_{k=1}^{{\overline{\rm K}}} \mathcal{P}^0(\mathcal{S}^0[k]|\mathcal{P}_0^0,\pi^0,k)$,
and ${\overline{\rm K}} = max \{k:\mathcal{P}^0(\mathcal{S}^0[k]\in\Lambda|\mathcal{P}_0^0,\pi^0,k)>0\}$,
then $\pi^0(\mathcal{A}^0[k]| \mathcal{S}^0[k])$ guarantees UUB stability of the control plant with ultimate bound $\rho$. If for any $\epsilon$, there is a $k>\epsilon$, such that  $\mathcal{P}^0(\mathcal{S}^0[k]\in\Lambda|\mathcal{P}_0^0,\pi^0,k)>0$, then ${\overline{\rm K}}= \infty$.
\end{theorem1}
Similar to \cite{han2021reinforcement}, a fully connected deep neural network is used to construct function $\Gamma_C(\mathcal{S}^0[k],\mathcal{A}^0[k])$, which satisfies $\Gamma(\mathcal{S}^0[k])= \mathbb{E}_{\mathcal{A}^0[k]\sim\pi^0}\Gamma_C(\mathcal{S}^0[k],\mathcal{A}^0[k])$ and the function $\Gamma_C$ is parameterized by $\upsilon$. A ReLU activation function is employed in the output layer of the deep neural network to guarantee positive output. To update $\upsilon$, the
following objective function is minimized:
\begin{equation}
\mathcal{L}(\upsilon)=\mathbb{\tilde{E}}_{\Lambda} (\dfrac{1}{2} (\Gamma_C(\mathcal{S}^0[k],\mathcal{A}^0[k])-\mathcal{C} (\mathcal{S}^0[k],\mathcal{A}^0[k]))^2),\\
\label{eqn:upsilon}
\end{equation}
where $\mathbb{\tilde{E}}_{\Lambda}(.)$ is the average over a mini-batch of samples collected from the sampling distribution $\Omega_{\overline{\rm K}}(s)$.\\
In the following, an approach based on option-critic method, which is reformulated as DAC architecture, is proposed to obtain $\pi^0$. In the obtaining procedure of policy $\pi^0$, we employ Theorem 1 to guarantee UUB stability.\\
\textbf{Option-Critic Method}:
Option-critic method is an HRL that has three policies: a master policy, an intra$-$option policy, and an option termination function \cite{zhang2019dac, funk2021learning}. The master policy decides which option should be performed. On the basis of this decision, an action is taken through intra$-$option policy until the option is terminated by the termination function. Accordingly, in the context of sub-problem (\ref{eqn:sub1}), the master policy specifies the probability of choosing update variable $\beta[k]$ at each time step $k$ and then control input $u[k]$ is determined by the intra$-$option policy. Furthermore, the termination function is omitted (similar to \cite{hou2020data} and \cite{funk2021learning}) because of the binary type of update variable $\beta[k]$. Indeed, when the master policy chooses one option ($\beta[k] = 1$ or $\beta[k] = 0$), it terminates another option simultaneously. Considering this performing model, we have:
 \begin{equation}
 \begin{split}
\mathcal{P}^0(\mathcal{S}^0[k+1] | \mathcal{S}^0[k], \beta[k]) = \sum_{a} \pi (a=&u[k] | \mathcal{S}^0[k],\beta[k])\times\\
\mathcal{P}^0(\mathcal{S}^0[k&+1] | \mathcal{S}^0[k],u[k]),\\
\mathcal{P}^0(\mathcal{S}^0[k+1], \beta[k+1]| \mathcal{S}^0[k], \beta[k]) =  \\\mathcal{P}^0 ( \mathcal{S}^0[k+1] | \mathcal{S}^0[k],\beta[k])\times
 \mathcal{P}^0(\beta [k+1]&| \mathcal{S}^0[k+1], \beta[k]).
 \label{eqn:proab}
 \end{split}
 \end{equation}
In \cite{zhang2019dac}, it is demonstrated that option-critic method can be reformulated as DAC architecture, which contains two augmented MDPs. The MDPs contain the high-level MPD,  $\mathcal{M}^{H_0}$, and the low-level MPD, $\mathcal{M}^{L_0}$, which are employed for choosing the option and the action respectively. Consequently, the high-level MPD of the base module is defined as:
\begin{equation}
\begin{split}
\mathcal{M}^{H_0}\!&\doteq \{\mathcal{S}^{H_0}; \mathcal{A}^{H_0}; \mathcal{R}^{H_0}; \mathcal{P}^{H_0}_0; \mathcal{P}^{H_0}\},\\
\mathcal{S}^{H_0}[k]\!&\doteq \{\beta[k-1], s^{c}[k] \},\\
\mathcal{A}^{H_0}[k]\!&\doteq \{\beta[k] \},\\
\mathcal{R}^{H_0}\!&\doteq \mathcal{R}^{H_0}(\mathcal{S}^{H_0}[k], \mathcal{A}^{H_0}[k])\\ &\doteq \mathcal{R}^0(\mathcal{S}^0[k], \beta[k]),\\
\mathcal{P}^{H_0}_0(\mathcal{S}^{H_0}[0])\!&\doteq \mathcal{P}_0 (\beta[-1], s^{c}[0]),\\
\mathcal{P}^{H_0} (\mathcal{S}^{H_0}[k+1]|\mathcal{S}^{H_0}[k],\mathcal{A}^{H_0}[k])\!&\doteq\!\textbf{1}_{\mathcal{A}^{H_0}=\beta[k]}\mathcal{P}^0 (\mathcal{S}^0[k+1]|\\&\qquad \mathcal{S}^0[k],\!\beta[k]),
\label{eqn:highlevel_C}
\end{split}
\end{equation}
where $\textbf{1}_{(.)}$ is the indicator function. Also, the high-level policy on $\mathcal{M}^{H_0}$ is defined as:
\begin{equation}
\pi^{H_0} (\mathcal{A}^{H_0}[k] | \mathcal{S}^{H_0}[k]) \doteq \mathcal{P}^0 (\beta[k] | \beta[k-1], s^{c}[k]).
\label{eqn:policy_C}
\end{equation}
The low-level MPD and policy of the base module are respectively stated as:
\begin{equation}
\begin{split}
\mathcal{M}^{L_0} &\doteq \{\mathcal{S}^{L_0}; \mathcal{A}^{L_0}; \mathcal{R}^{L_0}; \mathcal{P}^{L_0}_0; \mathcal{P}^{L_0}\},\\
\mathcal{S}^{L_0}[k] & \doteq \{s^{c}[k]\}\times \{\beta[k]\},\\
\mathcal{A}^{L_0}[k] & \doteq \{u[k] \},\\
\mathcal{R}^{L_0} & \doteq \mathcal{R}^{L_0}(\mathcal{S}^{L_0}[k], \mathcal{A}^{L_0}[k]) \\&\doteq \mathcal{R}^0(\mathcal{S}^0[k], \mathcal{A}^0[k]),\\
\mathcal{P}^{L_0}_0(\mathcal{S}^{L_0}[0]) & \doteq \mathcal{P}_0 (\mathcal{S}^0[0])\mathcal{P}_0(\beta[0]|\mathcal{S}^0[0]),\\
\mathcal{P}^{L_0} (\mathcal{S}^{L_0}[k+1]|\mathcal{S}^{L_0}[k],\mathcal{A}^{L_0}[k])& \doteq\\
\mathcal{P}^0 ((\beta[k+1], \mathcal{S}^0[k+1])| (\beta&[k], \mathcal{S}^0[k]), \mathcal{A}^{L_0}=u[k]) =\\
 \mathcal{P}^0 (\mathcal{S}^0[k+1]|\mathcal{S}^0[k],u[k]) \times\mathcal{P}^0 &(\beta[k+1] | \mathcal{S}^0[k+1],\beta[k]),
\label{eqn:Lowlevel}
\end{split}
\end{equation}
and
\begin{equation}
\pi^{L_0} (\mathcal{A}^{L_0}[k] | \mathcal{S}^{L_0}[k]) \doteq \mathcal{P}^0 (u[k] | \mathcal{S}^0[k],\beta[k]).
\label{eqn:policy_L}
\end{equation}
Considering trajectories of $\psi^0 = \{\tau^0| \mathcal{P}^0(\tau^0 | \pi^0 ,\mathcal{M}^0)\}$, $\psi^{H_0} = \{\tau^{H_0}| \mathcal{P}^{H_0}(\tau^{H_0} | \pi^{H_0} ,\mathcal{M}^{H_0})\}$ and $\psi^{L_0} = \{\tau^{L_0}| \mathcal{P}^{L_0}(\tau^{L_0} | \pi^{L_0} ,\mathcal{M}^{L_0})\}$, two bijection functions as $\mathcal{J}^{H_0} $ and $\mathcal{J}^{L_0} $ are obtained, which map $\tau^0$ to $\tau^{H_0}$ and $\tau^0$ to $\tau^{L_0}$ respectively. Here, the following lemmas holds, which appear similar to \cite{zhang2019dac}:
\begin{lemma}\label{L3}
Assuming the bijection function $\mathcal{J}^{H_0} $, we have $\mathcal{P}^0(\tau^0 | \pi^0 ,\mathcal{M}^0) = \mathcal{P}^{H_0}(\tau^{H_0} | \pi^{H_0} ,\mathcal{M}^{H_0})$ and $\mathcal{R}^0 (\tau^0)= \mathcal{R}^{H_0} (\tau^{H_0})$.
\end{lemma}
\begin{figure}[h!]
\centerline{\includegraphics[width=0.25\textwidth]{Fig3}}
\label{fig:actor-critic architecture}
\caption{DAC architecture of the base module}
\end{figure}
\begin{lemma}\label{L4}
	Assuming the bijection function $\mathcal{J}^{L_0} $, we have $\mathcal{P}^0(\tau^0 | \pi^0 ,\mathcal{M}^0) = \mathcal{P}^{L_0}(\tau^{L_0} | \pi^{L_0} ,\mathcal{M}^{L_0})$ and $\mathcal{R}^0 (\tau^0)= \mathcal{R}^{L_0}(\tau^{L_0})$.
\end{lemma}
The proof of above lemmas are provided in Appendices A and B respectively. These lemmas specify that $\{\pi^{H_0} ,\mathcal{M}^{H_0}\}$ and $ \{\pi^{L_0} ,\mathcal{M}^{L_0}$\} can share the same samples with $\{\pi^0 ,\mathcal{M}^0\}$. In the same way of the provided proof, Theorem 2 can be simply driven as follows.
\vspace{-10pt}
\begin{theorem2}
\begin{equation}
\begin{split}
 \Phi^0&\doteq \int \mathcal{R}^0 (\tau^0) \mathcal{P}^0(\tau^0 | \pi^0 ,\mathcal{M}^0)d\tau^0 \\
&= \int \mathcal{R}^{H_0} (\tau^{H_0}) \mathcal{P}^{H_0}(\tau^{H_0} | \pi^{H_0} ,\mathcal{M}^{H_0})d\tau^{H_0}\\
&=  \int \mathcal{R}^{L_0} (\tau^{L_0}) \mathcal{P}^{L_0}(\tau^{L_0} | \pi^{L_0} ,\mathcal{M}^{L_0})d\tau^{L_0}.
\label{eqn:eq-objectives}
\end{split}
\end{equation}
\end{theorem2}
\vspace{-1.5mm}
\noindent Following Lemma \ref{L3}, Lemma \ref{L4}, and Theorem 2, to handle sub-problem (\ref{eqn:sub1}), the learning agent alternately optimize $\pi^{H_0}$ (decide on option variable $\beta[k]$) and $\pi^{L_0}$ (decide on continues variable $u[k]$). Thereby option-critic method is reformulated to DAC architecture (see Fig. 3). To optimize policies in each augmented MDP ($\mathcal{M}^{H_0}$, $\mathcal{M}^{L_0}$), Proximal Policy Optimization (PPO) method is used similar to \cite{zhang2019dac}.\\
In DAC architecture, two parameterized polices $\pi^{H_0}(\mathcal{S}^{H_0}[k], \mathcal{A}^{H_0}[k],\theta^{H_0})$, which is summarized as $\pi_{\theta^{H_0}}(\mathcal{S}^{H_0}[k], \mathcal{A}^{H_0}[k])$ and $\pi^{L_0}(\mathcal{S}^{L_0}[k], \mathcal{A}^{L_0}[k],\theta^{L_0})$, which is summarized as $\pi_{\theta^{L_0}}(\mathcal{S}^{L_0}[k], \mathcal{A}^{L_0}[k])$ are estimated in two actor neural networks. Additionally, the parameterized value function ($V^{0}(\mathcal{S}^{0}[k],\omega^{0})$) and parameterized function $\Gamma_C(\mathcal{S}^0[k],\mathcal{A}^0[k],\upsilon)$ are estimated in two critic neural networks.
 Therefore, to update two parameters $\theta^{H_0}$ and $\theta^{L_0}$, considering UUB stability constraint, the following objective functions are minimized respectively:
\begin{equation}
\begin{split}
\mathcal{L}_{H_0}(\theta^{H_0}) &= \mathbb{\tilde{E}} [min (\frac{\pi_{\theta^{H_0}}(\mathcal{A}^{H_0}[k] | \mathcal{S}^{H_0}[k])} {\pi_{\theta^{H_0}_{old}}(\mathcal{A}^{H_0}[k] | \mathcal{S}^{H_0}[k])} ,\\& clip (\frac{\pi_{\theta^{H_0}}(\mathcal{A}^{H_0}[k] | \mathcal{S}^{H_0}[k])} {\pi_{\theta^{H_0}_{old}}(\mathcal{A}^{H_0}[k] | \mathcal{S}^{H_0}[k])}, 1-\epsilon, 1+\epsilon)) \widehat{A}_0[k]]\\& + \lambda \mathbb{\tilde{E}}_{\Lambda}[\Gamma_C(\mathcal{S}^0[k+1],\mathcal{A}^0[k+1])\mathbb{I}_{\Lambda}(\mathcal{S}^0[k+1])\\&-(\Gamma_C(\mathcal{S}^0[k],\mathcal{A}^0[k])-\alpha_3\mathcal{C}_{\pi^0}(\mathcal{S}^0[k]))\mathbb{I}_{\Lambda}(\mathcal{S}^0[k])],
\label{eqn:eq-objective func1}
\end{split}
\end{equation}
and
\vspace{-8pt}
\begin{equation}
\begin{split}
\mathcal{L}_{L_0}(\theta^{L_0}) =& \mathbb{\tilde{E}} [min (\frac{\pi_{\theta^{L_0}}(\mathcal{A}^{L_0}[k] | \mathcal{S}^{L_0}[k])} {\pi_{\theta^{L_0}_{old}}(\mathcal{A}^{L_0}[k] | \mathcal{S}^{L_0}[k])} ,\\& clip (\frac{\pi_{\theta^{L_0}}(\mathcal{A}^{L_0}[k] | \mathcal{S}^{L_0}[k])} {\pi_{\theta^{L_0}_{old}}(\mathcal{A}^{L_0}[k] | \mathcal{S}^{L_0}[k])}, 1-\epsilon, 1+\epsilon)) \widehat{A}_0[k]]\\&+\lambda \mathbb{\tilde{E}}_{\Lambda}[\Gamma_C(\mathcal{S}^0[k+1],\mathcal{A}^0[k+1])\mathbb{I}_{\Lambda}(\mathcal{S}^0[k+1])\\&-(\Gamma_C(\mathcal{S}^0[k],\mathcal{A}^0[k])-\alpha_3\mathcal{C}_{\pi^0}(\mathcal{S}^0[k]))\mathbb{I}_{\Lambda}(\mathcal{S}^0[k])],
\label{eqn:eq-objective func2}
\end{split}
\end{equation}
where $\mathbb{\tilde{E}}(.)$ is the average over a mini-batch of samples (the size of the mini-batch is $\mathcal{Y}^0$), function $clip(.)$ constrains the ratio of $\frac{\pi_{\theta^{L_0}}(\mathcal{A}^{L_0}[k] | \mathcal{S}^{L_0}[k])} {\pi_{\theta^{L_0}_{old}}(\mathcal{A}^{L_0}[k] | \mathcal{S}^{L_0}[k])}$ between the interval of $(1-\epsilon, 1+\epsilon)$, and $\epsilon$ is a hyper-parameter. $\lambda$ is a positive Lagrangian multiplier, which is adjusted via gradient ascent to maximize the following objective function \cite{han2021reinforcement}:
\begin{equation}
\begin{split}
\mathcal{L}(\lambda) = &\lambda \mathbb{\tilde{E}}_{\Lambda}[\Gamma_C(\mathcal{S}^0[k+1],\mathcal{A}^0[k+1])\mathbb{I}_{\Lambda}(\mathcal{S}^0[k+1])\\&-(\Gamma_C(\mathcal{S}^0[k],\mathcal{A}^0[k])-\alpha_3\mathcal{C}_{\pi^0}(\mathcal{S}^0[k]))\mathbb{I}_{\Lambda}(\mathcal{S}^0[k])].
\label{eqn:lambda}
\end{split}
\end{equation}In equations (\ref{eqn:eq-objective func1}) and (\ref{eqn:eq-objective func2}), $\widehat{A}_0[k]$ is the advantage function at time step $k$ and is estimated via the Generalized Advantage Estimation (GAE) as:
\begin{equation}
\begin{split}
\widehat{A}_0[k] = & \delta^0 [k] + (\gamma \xi)\delta^0 [k+1] (\gamma \xi)\delta^0 [k+1]+...\\&+(\gamma \xi)^{(\mathcal{Y}^0-k+1)}\delta^0 [\mathcal{Y}^0-1],
\label{eqn:advantag_0}
\end{split}
\end{equation}
where $\xi \in [0,1]$ is the GAE parameter and $\delta^0[k]$ is the temporal difference(TD) error, given by $\delta^0[k] = \mathcal{R}^0[k+1] +\gamma V^0_{\omega^0} (\mathcal{S}^0[k+1]) - V^0_{\omega^0} (\mathcal{S}^0[k])$. $\omega^0$ is updated by an Stochastic Gradient Descent (SGD) algorithm as:
\begin{equation}
\omega^0 = \omega^0 - \zeta_{\omega^0} \nabla L^V(\omega^0),
\label{eqn:updat_omaga}
\end{equation}
where $\zeta_{\omega^0}$ is the learning rate and $L^V(\omega^0)$ is the objective function calculated as:\\
\begin{equation}
L^V(\omega^0) = \mathbb{\tilde{E}}[ \vert \widehat{V}^{\text{target}}_{\omega^0}(\mathcal{S}^0[k] )- V^0_{\omega^0}(\mathcal{S}^0[k] )\vert ],
\label{eqn:objective_omaga}
\end{equation}
where $\widehat{V}^{\text{target}}_{\omega^0}(\mathcal{S}^0[k] ) = \mathcal{R}^0[k+1] + \gamma V^0_{\omega^0}(\mathcal{S}^0[k+1])$ is the target value of time-difference error.\\
In summary, at time step $k$, each actor network selects its action according to its current state using its policy. This leads to the state transition to $\mathcal{S}^{0}[k+1]$ and a new reward value, which is estimated by the critic network via value function $V^{0}(\mathcal{S}^{0}[k],\omega^{0})$. Afterward, the TD error is calculated, which is the critic network feedback to optimize $\pi^{H_0}$ and $\pi^{L_0}$ using optimizing problems (\ref{eqn:eq-objective func1}) and (\ref{eqn:eq-objective func2}). In addition, the selected actions and the state transition ($\mathcal{S}^{0}[k+1]$ and $\mathcal{A}^0[k+1]$) lead to updating $\upsilon$ and $\lambda$ through (\ref{eqn:upsilon}) and (\ref{eqn:lambda}). Then, the updated $\Gamma_C(.)$ is sent to actor networks as feedback of $\Gamma_C(.)-$Critic Network to optimize $\pi^{H_0}$ and $\pi^{L_0}$.
\vspace{-10pt}
\subsection{The First Attribute Module}
Having decided update variable $\beta[k]$, DRL policy of the first attribute module is trained to address sub-problem (\ref{eqn:sub2}). The state and action spaces of this module are given by:
\begin{equation}
\begin{split}
\mathcal{S}^1[k] &= \{s^{N}_n[k] \},\\
\mathcal{A}^1[k] &= \{\textbf{A}[k], \textbf{P}[k]\}.
\label{eqn:state_first}
\end{split}
\end{equation}
Considering the objective and the constraints of sub-problem (\ref{eqn:sub2}), the reward function of the first module is calculated as:
\begin{equation}
\begin{split}
\mathcal{R}^{1}[k]  = \left\{
\begin{array}{rl}
&-\mu_2 \quad\qquad \overline{P}_\text{BS}\leq\sum_{n=1}^\mathcal{N}\sum_{j=1}^\mathcal{J}p_{n,j}[k]a_{n,j}[k] \\
&-\mathcal{P}^\text{BS}_\text{total}+ \mu_3 \sum_{n=1}^{N}s^{N}_n[k] \quad \quad\qquad o/w,
\end{array}\right.
\label{eqn:reward_first}
\end{split}
\end{equation}
where $\mu_2$ is a hyper-parameter denoted the penalty weight on
crossing the limitation of the BS' power consumption. Also, $\mu_3$ is a hyper-parameter indicated the weight on the number of users received their required rate.\\
\textbf{Option-Critic Method}: The action space of this module is also hybrid. Accordingly, to address sub-problem (\ref{eqn:sub2}), option-critic method reformulated as DAC architecture is employed too. Assuming sub-carrier allocation matrix $\textbf{A}[k]$ is the option variable, we have:
\begin{equation}
\begin{split}
\mathcal{P}^1(\mathcal{S}^1[k+1] | \mathcal{S}^1[k], \textbf{A}[k]) = \sum_{a} \pi (a=&\textbf{P}[k] | \mathcal{S}^1[k],\textbf{A}[k])\times\\
\mathcal{P}^1(\mathcal{S}^1[k&+1] | \mathcal{S}^1[k],\textbf{P}[k]),\\
\mathcal{P}^1(\mathcal{S}^1[k+1], \textbf{A}[k+1]| \mathcal{S}^1[k], \textbf{A}[k]) =  \\\mathcal{P}^1 ( \mathcal{S}^1[k+1] | \mathcal{S}^1[k],\textbf{A}[k])\times
\mathcal{P}^1(\textbf{A}[k+1] &| \mathcal{S}^1[k+1], \textbf{A}[k]).
\label{eqn:proab_first}
\end{split}
\end{equation}
Option-critic method can be reformulated as two augmented MDPs. The high-level MPD ($\mathcal{M}^{H_1}$) is used for the sub-carrier assignment ($\textbf{A}[k]$) and the low-level MPD ($\mathcal{M}^{L_1}$) is used for the power allocation ($\textbf{P}[k]$). The high-level MPD of the first module is given as:
\vspace{-6pt}
\begin{equation}
\begin{split}
\mathcal{M}^{H_1}\doteq &\{\mathcal{S}^{H_1};\mathcal{A}^{H_1}; \mathcal{R}^{H_1};\\ &\quad \mathcal{P}_0^{H_1};\mathcal{P}^{H_1}\},\\
\mathcal{S}^{H_1}[k]\doteq &\{\textbf{A}[k-1], s^{N}_n[k] \},\\
\mathcal{A}^{H_1}[k]\doteq &\{\textbf{A}[k] \},\\
\mathcal{R}^{H_1}\doteq &\mathcal{R}^{H_1}(\mathcal{S}^{H_1}, \mathcal{A}^{H_1})\\\doteq &\mathcal{R}^1(\mathcal{S}^1[k],  \textbf{A}[k]),\\
\mathcal{P}_0^{H_1}(\mathcal{S}^{H_1}[0])\doteq &\mathcal{P}_0^1 ((\textbf{A}[-1], s^{N}_n[0])),\\
\mathcal{P}^{H_1} (\mathcal{S}^{H_1}[k+1]|\mathcal{S}^{H_1}[k],\mathcal{A}^{H_1}[k] )\doteq&\textbf{1}_{\mathcal{A}^{H_1}[k]=\textbf{A}[k]}\mathcal{P}^1 ( \mathcal{S}^1[k+1]|\\& \mathcal{S}^1[k],\textbf{A}[k]).
\label{eqn:highlevel_Cfirst}
\end{split}
\end{equation}
\vspace{-2pt}
And the high-level policy of the first module on $\mathcal{M}^{H_1}$ is defined as:
\begin{equation}
\pi^{H_1} (\mathcal{A}^{H_1}[k] | \mathcal{S}^{H_1}[k]) \doteq \mathcal{P}^1 (\textbf{A}[k] | \textbf{A}[k-1], s^{N}_n[k]).
\label{eqn:policy_C_first}
\end{equation}
The low-level MPD and policy of the first module are respectively defined as:
\begin{equation}
\begin{split}
\mathcal{M}^{L_1} \doteq& \{\mathcal{S}^{L_1}; \mathcal{A}^{L_1}; \mathcal{R}^{L_1}; \mathcal{P}^{L_1}_0;\\ & \quad \mathcal{P}^{L_1}\},\\
\mathcal{S}^{L_1}[k] \doteq& \{s^{N}_n[k]\}\times \{\textbf{A}[k]\},\\
\mathcal{A}^{L_1}[k]  \doteq& \{\textbf{P}[k] \},\\
\mathcal{R}^{L_1} \doteq& \mathcal{R}^{L_1}(\mathcal{S}^{L_1}, \mathcal{A}^{L_1})\\=& \mathcal{R}^{L_1}((\mathcal{S}^{L_1},\textbf{A}[k]),\textbf{P}[k]) \\\doteq& \mathcal{R}^1(\mathcal{S}^1[k], \mathcal{A}^1[k]),\\
\mathcal{P}^{L_1}_0(\mathcal{S}^{L_1}[0]) \doteq& \mathcal{P}^1_0 (\mathcal{S}^1[0])\mathcal{P}^1_0(\textbf{A}[0]|\mathcal{S}^1[0]),\\
\mathcal{P}^{L_1} (\mathcal{S}^{L_1}[k+1]|\mathcal{S}^{L_1}[k],\mathcal{A}^{L_1}[k]) \doteq&\\
\mathcal{P}^1 ((\textbf{A}[k+1], \mathcal{S}^1[k+1])| (\textbf{A}[k]&, \mathcal{S}^1[k]), \mathcal{A}^{L_1}=\textbf{P}[k]) \\
= \mathcal{P}^1 (\mathcal{S}^1[k+1]|\mathcal{S}^1[k],\textbf{P}[k]) \times\mathcal{P}^1 &(\textbf{A}[k+1] | \mathcal{S}^1[k+1],\textbf{P}[k]),
\label{eqn:Lowlevel_first}
\end{split}
\end{equation}
and
\begin{equation}
\pi^{L_1} (\mathcal{A}^{L_1}[k] | \mathcal{S}^{L_1}[k]) \doteq \mathcal{P}^1 (\textbf{P}[k] | \mathcal{S}^1[k], \textbf{A}[k]).
\label{eqn:policy_L_first}
\end{equation}
Considering trajectories of $\psi^1 = \{\tau^1| \mathcal{P}^1(\tau^1 | \pi^1 ,\mathcal{M}^1)\}$, $\psi^{H_1} = \{\tau^{H_1}| \mathcal{P}^{H_1}(\tau^{H_1} | \pi^{H_1} ,\mathcal{M}^{H_1})\}$ and $\psi^{L_1} = \{\tau^{L_1}| \mathcal{P}^{L_1}(\tau^{L_1} | \pi^{L_1} ,\mathcal{M}^{L_1})\}$, two bijection functions as $\mathcal{J}^{H_1} $ and $\mathcal{J}^{L_1} $ can be found, which map $\tau^1$ to $\tau^{H_1}$ and $\tau^1$ to $\tau^{L_1}$ respectively. Similar lemmas to Lemma \ref{L3} and Lemma \ref{L4} can be stated, which indicate that $\{\pi^{H_1} ,\mathcal{M}^{H_1}\}$ and $ \{\pi^{L_1} ,\mathcal{M}^{L_1}$\} can share the same samples with $\{\pi^1 ,\mathcal{M}^1\}$. Similarly, Theorem 3 can be driven as:
\begin{theorem3}
\begin{equation}
\begin{split}
\Phi^1&\doteq \int \mathcal{R}^1 (\tau^1) \mathcal{P}^1(\tau^1 | \pi^1 ,\mathcal{M}^1)d\tau^1 \\
&= \int \mathcal{R}^{H_1} (\tau^{H_1}) \mathcal{P}^{H_1}(\tau^{H_1} | \pi^{H_1} ,\mathcal{M}^{H_1})d\tau^{H_1}\\
&=  \int \mathcal{R}^{L_1} (\tau^{L_1}) \mathcal{P}^{L_1}(\tau^{L_1}| \pi^{L_1} ,\mathcal{M}^{L_1})d\tau^{L_1}.
\label{eqn:eq-objectives_first}
\end{split}
\end{equation}
\end{theorem3}
Therefore, to solve sub-problem (\ref{eqn:sub2}), the learning agent alternatively optimize $\pi^{H_1}$ (decide on sub-carrier assignment $\textbf{A}[k]$ in $\mathcal{M}^{H_1}$) and $\pi^{L_1}$ (decide on power allocation $\textbf{P}[k]$ in $\mathcal{M}^{L_1}$). To optimize the policies in each augmented MDP of the first module ($\pi^{H_1}$ and $\pi^{L_1}$), PPO method is employed analogous to the previous sub-section.\\
\indent The proposed approach to address problem (\ref{eqn:optt}) is called Cascade Stable Double Actor-Critic (CSDAC). In summary, Algorithm 1 provides the pseudo-code of CSDAC approach.
\begin{algorithm}\label{nfv-alg}
	\BlankLine
	\SetKwFunction{Range}{range
	\SetKw{KwTo}{in}\SetKwFor{For}{for}{\string:}{}%
	\SetKwIF{If}{ElseIf}{Else}{if}{:}{elif}{else:}{}%

	\SetAlgoNoEnd

	\SetAlgoNoLine%

	\textbf{Input:} The mini-batch sizes of $\mathcal{Y}^{\Lambda}$, $\mathcal{Y}^0$ and $\mathcal{Y}^1$.\\
	Set iteration index $i \leftarrow  0 $.\\
	\textbf{Repeat}\\
	\Indp
	for each time step do:\\
	\Indp
	Observe state $\mathcal{S}^0[k]$.\\
	Execute action $\mathcal{A}^{H_0}[k]$ and $\mathcal{A}^{L_0}[k]$ and get reward $ \mathcal{R}^0[k+1]$ (according to (\ref{eqn:reward_base})) and transmit to $\mathcal{S}^0[k+1]$. \\
	Store transition $(\mathcal{S}^0[k], \mathcal{A}^{0}[k], \mathcal{R}^0[k+1], \mathcal{S}^0[k+1])$ in $\mathcal{D}^0$.\\
	Record the largest instant $\overline{\rm K}$ at which $s\in\Lambda$ and Store transition $(\mathcal{S}^0[k], \mathcal{A}^{H_0}[k], \mathcal{A}^{L_0}[k], \mathcal{C}[k], \mathcal{R}^0[k+1], \mathcal{S}^0[k+1]): k<\overline{\rm K}$ in $\mathcal{D}^{\Lambda} $.\\
	Sample mini-batches of transitions from $\mathcal{D}^0$ and $\mathcal{D}^{\Lambda}$ and for each update step do:\\
	\Indp
	\Indp
	Optimize $\pi_{\theta^{H_0}}$ (update parameter $\theta^{H_0}$ through the minimization of (\ref{eqn:eq-objective func1})).\\
	Optimize $\pi_{\theta^{L_0}}$ (update parameter $\theta^{L_0}$ through the minimization of (\ref{eqn:eq-objective func2})).\\
	Optimize $\Gamma_C(.)$ is (update parameter $\upsilon$ through the minimization of (\ref{eqn:upsilon})).\\
	Update $\omega^0$ (according to (\ref{eqn:updat_omaga})) and Calculate the TD error $\delta^0[k]$.\\
	Optimize $\mathcal{L}(\lambda)$ (update parameter $\lambda$ through the minimization of (\ref{eqn:lambda})).\\
	\Indm
	\Indm
	Observe $\beta[k]$ and $\mathcal{S}^1[k]$.\\
	Execute action $\mathcal{A}^{H_1}[k]$ and $\mathcal{A}^{L_1}[k]$ and get reward  $\mathcal{R}^1[k+1]$ (according to (\ref{eqn:reward_first})) and transmit to $\mathcal{S}^1[k+1]$.\\
	Store transition $(\mathcal{S}^1[k], \mathcal{A}^{H_1}[k], \mathcal{A}^{L_1}[k], \mathcal{R}^1[k+1], \mathcal{S}^1[k+1])$ in $\mathcal{D}^1$.\\
	Sample mini-batches of transitions from $\mathcal{D}^1$ and for each update step do:\\
	\Indp
	\Indp
	Optimize $\pi^{H_1}$ (similar to $\pi_{\theta^{H_0}}$). \\
	Optimize  $\pi^{L_1}$(similar to $\pi_{\theta^{L_0}}$).\\
	Update $\omega^1$ (similar to (\ref{eqn:updat_omaga})) and Calculate TD error $\delta^1[k]$.\\
	\Indm
	\Indm
	\Indm
	end for.\\
	$i \leftarrow i + 1$.\\
	\Indm
	\textbf{until}(\ref{eqn:stability}) is satisfied and $i$ exceeds a designed threshold.
	\caption{Pseudocode of Cascade Stable Double Actor-Critic (CSDAC) }
\end{algorithm}
\section{Validation Results}
In this section, we evaluate the general applicability and efficacy of our proposed approach (CSDAC) through various simulation results. First, we show the general applicability of CSDAC through a rather simple simulated ICPS including OpenAI Gym Cart-Pole environment \cite{gym} communicating over a single cell OFDMA. This low-dimensional environment allows us to compare the performance of CSDAC with a disjoint design of control and communication methods based on a classical ETC method. Second, the simulated ICPS consists of PyBullet-Gym Ant environment \cite{pybullet} communicating over a single cell OFDMA. We show the capability of CSDAC in this challenging high-dimensional environment. For both simulated ICPSs, we consider a single cell OFDMA as their network sub-systems with the same specifications. In the downlink model of each cellular network, eMBB and URLLC users are considered randomly dispersed in a square cell. The channel gain for each user $n$ on sub-carrier $j$ at time step $k$ is calculated as $g_{n,j}[k] = hd_n^{-3}[k]$, where $d_n[k]$ is the distance between user $n$ and the BS as well as $h = 0.09$ is the loss factor. The distance is fixed for URLLC users; however, it varies for eMBB users. Moreover, we assume that eMBB users move only within the BS coverage area during the simulation time. The other parameters of the simulated networks are given in Table II.
\begin{table}\begin{center}
		\caption{Values of Simulation Parameters}
		\begin{tabular}{||c  c ||}
			\hline
			Parameter & Value \\ [0.7ex]
			\hline\hline
			Number of eMBB users & 8\\

			\hline
			Bandwidth of sub-carrier ($\omega$) & 180 KHz \\
			\hline
			Maximum transmission power of the BS ($\overline{P}_{BS}$) & 44 dBm\\
			\hline
			Constant power of the BS ($P_{cst}$) & 20 dBm  \\
			\hline
			Noise power ($N_0[k]$)  & -62 dBm \\
			\hline
			Distance between a user and the BS ($d_n[k]$) & 10-50 m \\
			\hline
			Minimum data-rate for an RC user ($\overline{R}^e_{n}[k]$) & 100 bit/s \\
			\hline
			Packet length of control users ($L_c$)  & 70 bit \\
			\hline
			Number of sub-carriers ($J$)  & 8\\ [1ex]
			\hline
		\end{tabular}
		\label{tab: Table11}
\end{center}\end{table}
\vspace{-15pt}
\subsection{Cart-Pole Environment}
\vspace{-3pt}
Cart-Pole is a classical control problem, which includes a pole placed on a cart. This system is a two-degree-of-freedom system containing the linear movement of the cart through the X axis ($x_d$) and the rotational movement of the pole on the X-Y axes ($x_r$). The goal of the control problem is to keep the pole upright via moving the cart left/right. Each episode of the Gym simulation ends when the rotational movement of the pole is more than 0.261 radians or the linear movement of the cart is more than 2.4 units from the origin.\\
As a proof of concept, we apply CSDAC, described in Sec.III, to jointly design the event-triggered control and the resource allocation in the simulated ICPS including Cart-Pole environment communicating over the single cell OFDMA. Moreover, the efficacy of CSDAC is compared to a disjoint and model-based method, which is a combination of the suggested methods in  \cite{luo2019periodic} and \cite{kai2019joint}. We consider the control reward ($\mathcal{R}^{ctrl}$) is identical to the default reward provided by the Gym environment. The constraint function is calculated by $\mathcal{C}= max(||x_d[k]||-1.1,0)$. To estimate each hyper-parameter in equations (\ref{eqn:reward_base}) and (\ref{eqn:reward_first}), 25 values between 0.01 and 100 are examined. The efficiency of each hyper-parameter is analyzed through 5 randomized training processes using different random inputs. We evaluate the performance of CSDAC by carrying out 50 randomized test episodes and each episode is terminated after 300 discrete time steps. Fig. 4 and Fig. 5 depict the results from one test episode. Fig. 4 illustrates the control state responses obtained with the learning event-triggered controller. Having 46 updates on control input for 25s, the state variables of Cart-Pole system remain stable within the range.\\
To compare CSDAC to a disjoint and model-based method, we use Matlab/Simulink applying a linear model of Cart-Pole as \cite{shi2009output}:\\
\begin{equation}
\begin{split}
\begin{bmatrix}
\dot{x_d}(t)\\
\ddot {x}_d(t)\\
\dot {x_r}(t)\\
\ddot {x}_r(t)
\end{bmatrix}
\!=& \!\begin{bmatrix}
1 & 0.1 & -0.0166 & -0.0005 \\
0 & 1 & -0.3374 & -0.0166\\
0 & 0 & 1.0996 & 0.1033 \\
0 & 0 & 2.0247 & 1.0996
\end{bmatrix}\!\begin{bmatrix}
x_d(t)\\
\dot {x}_d(t)\\
x_r(t)\\
\dot{x}_r(t)
\end{bmatrix}\!+\!\begin{bmatrix}
0.0045\\
0.0896\\
-0.0068\\
-0.1377\\
\end{bmatrix}\!u,\\
\begin{bmatrix}
x_d(t)\\
x_r(t)
\end{bmatrix} =& \begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 \\
\end{bmatrix}
(\begin{bmatrix}
x_d(t),
\dot{x}_d(t),
x_r(t),
\dot{x}_r(t)
\end{bmatrix})^T,
\label{eqn:pl2}
\end{split}
\end{equation}
where $x_d(t)$ is the horizontal displacement of the cart from origin, $x_r(t)$ is the rotational movement of the pole, $u(t)$ denotes the enforced control input to the cart. We employ the event-triggered method proposed in \cite{luo2019periodic} with the linear quadratic regulation method to guarantee the closed-loop stability and an upper bound of performance, which is defined in \cite{luo2019periodic}. Moreover, the optimal resource allocation is carried out using MATLAB/CVX similar to \cite{kai2019joint}. Indeed, we start with the sub-carrier assignment for a given power matrix. Next, the power matrix is allocated while using the previous step sub-carrier assignment matrix. These two steps are done iteratively up to reach a convergence criteria.\\
The simulation results, illustrated in Fig. 4 and Fig. 6, show both the model-based method (in accordance with \cite{luo2019periodic, kai2019joint}) and CSDAC are properly controlled the ICPS. However, the total number of updates on the actuator's input, under CSDAC algorithm and the model-based method, are respectively 46 and 102.
In other words, the number of updates on the actuator's input is reduced by around 55\% when CSDAC is employed. Also, the downlink power consumption in CSDAC algorithm and the model-based method are compared in Fig. 8, where a roughly 64\% downlink power saving is observed.
\begin{figure}[h!]
	\centerline{\includegraphics[width=0.5\textwidth]{fig4}}
	\label{fig:learnig_cartpol}
	\caption{State variables' response obtained with CSDAC}
\end{figure}
\begin{figure}[h!]

	\centerline{\includegraphics[width=0.5\textwidth]{fig5}}
	\label{fig:learnig_cartpolg}
	\caption{Variation of control input and option variables in CSDAC}
\end{figure}
\begin{figure}[h!]

	\centerline{\includegraphics[width=0.5\textwidth]{fig6}}
	\label{fig:resp model-based event-triggered}
	\caption{State variables' response obtained with the model-based method \cite{luo2019periodic, kai2019joint}}
\end{figure}
\begin{figure}[h!]

	\centerline{\includegraphics[width=0.5\textwidth]{fig7}}
	\label{fig:control model-based event-triggered}
	\caption{Variation of the control input in the model-based method \cite{luo2019periodic, kai2019joint}}
\end{figure}
\begin{figure}[h!]

 	\centerline{\includegraphics[width=0.5\textwidth]{fig8}}
 	\label{fig:power_compare}
 	\caption{Compare downlink power consumption in CSDAC with the model-base method (in accordance with \cite{luo2019periodic, kai2019joint})}
 \end{figure}
 \begin{figure}[h!]
	\centering

	\centerline{\includegraphics[width=0.25\textwidth]{fig9}}
	\label{fig:ant}
	\caption{Ant environment screen-shot}
\end{figure}
\subsection{Ant Environment}
We now focus on PyBullet-Gym Ant environment. In such nonlinear and high-dimensional environments, employing known ETC methods is usually unsuccessful as the settings are too complicated \cite{funk2021learning}.
Ant is a three-dimensional quadruped robot (see Fig. 9), which is rewarded to run forward as fast as possible while a safety constraint on its speed ($v<2.3$) needs to be ensured. We apply CSDAC algorithm in the simulated ICPS including Ant environment communicating over the single cell OFDMA. Moreover, the capability of CSDAC is compared to the soft actor$-$critic (SAC) algorithm \cite{haarnoja2018soft}. We assume the control reward ($\mathcal{R}^{ctrl}$) is the same with
the default reward in the PyBullet-Gym. Also, the safety
constraint on forwarding speed is $\mathcal{C}= max(||v[k]||-2.3,0)$. To estimate each hyper-parameter in equations (\ref{eqn:reward_base}) and (\ref{eqn:reward_first}), we follow the same procedure described in the previous sub-section. Also, each episode is terminated after 300 discrete time steps.\\
As illustrated in Fig. 10, CSDAC performs stably with respect to the safety constraint. On the other hand, SAC fails to find UUB stable policy concerning the safety constraint. The simulation results shown in Fig. 10, Fig. 11, and Fig. 12, demonstrate that the number of constraint violations is approximately zero in CSDAC policy while maintaining reasonable rewards of the base and first modules. In terms of downlink power consumption saving and users' QoS satisfying, as shown in Fig. 12, CSDAC performs much better than SAC and returns significantly more reward in the first module.
\begin{figure}[h!]
	\centering

	\centerline{\includegraphics[width=0.5\textwidth]{fig10}}
	\label{fig:const_break}
	\caption{Number of constraint violations in CSDAC compared to SAC}
\end{figure}
\begin{figure}[h!]
	\centering

	\centerline{\includegraphics[width=0.5\textwidth]{fig11}}
	\label{fig:base_reward}
	\caption{The base module's reward trained by CSDAC and SAC}
\end{figure}
\begin{figure}[h!]
	\centering

	\centerline{\includegraphics[width=0.5\textwidth]{fig12}}
	\label{fig:first_module_reward}
	\caption{The first module's reward trained by CSDAC and SAC}
\end{figure}
\section{Conclusion and Future Work}
This paper has concerned with the joint design problem of the event-triggered control and the energy-efficient resource allocation in the 5G-based network. So, a multi-objective optimization problem was formulated to minimize both the number of updates on the actuators' input and the total downlink power usage. The problem constraints contained the dynamics and UUB stability of the control plant, the QoS demands of eMBB and URLLC users, power and sub-carrier constraints, and the BS transmit power limitation. We proposed a novel model-free HRL approach with UUB stability guarantee using CAN method to address the problem. In CAN, the problem was decoupled into two sub-problems. In this regard, we proved that using the decoupling method resulted in a Pareto solution. As the action space of each sub-problem was hybrid, we used DAC architecture to deal with each sub-problem. We demonstrated the effectiveness of the proposed approach by simulation results. Considering joint uplink and downlink transmission will be concentrated on in future work.
\begin{appendices}
\section{Proof of Lemma\ref{L3}}
\vspace{-1mm}
In the higher-level MDP of the base module, we define:\\
\begin{equation*}
 \quad \quad \quad \mathcal{A}^{H_0}[k]\doteq \{\beta[k] \}\\
\end{equation*}
\begin{equation*}
\mathcal{P}^{H_0} (\mathcal{S}^{H_0}[k+1]|\mathcal{S}^{H_0}[k],\mathcal{A}^{H_0}[k]) \doteq\!\textbf{1}_{\mathcal{A}^{H_0}=\beta[k]}\mathcal{P}^0 (\mathcal{S}^0[k+1]|\mathcal{S}^0[k],\!\beta[k]),\\
\end{equation*}
and
\begin{equation*}
\pi^{H_0} (\mathcal{A}^{H_0}[k] | \mathcal{S}^{H_0}[k]) \doteq \mathcal{P}^0 (\beta[k] | \beta[k-1], s^{c}[k]),
\end{equation*}
therefore:
\begin{equation*}
\mathcal{P}^0(\tau^0 | \pi^0 ,\mathcal{M}^0) = \mathcal{P}^0 (\mathcal{S}^0[0]) \prod_{k= 0}^{K-1} (\mathcal{P}^0(\beta[k] | \mathcal{S}^0[k],\beta[k-1]) \mathcal{P}^0(\mathcal{S}^0[k+1]| \\
\end{equation*}
\begin{equation*}
\mathcal{S}^0[k],\beta[k])) = \mathcal{P}^0 (\mathcal{S}^0[0]) \prod_{k= 0}^{K-1} (\mathcal{P}^0(\beta[k] | S^{0}[k],\beta[k-1])\textbf{1}_{\mathcal{A}^{H_0}=\beta[k]} \\
\end{equation*}
\begin{equation*}
\mathcal{P}^0(S^{0}[k+1]|S^{0}[k],\beta[k])) = \mathcal{P}^{H_0} (S^{H_0}[0]) \prod_{k= 0}^{K-1}(\pi^{H_0} (\mathcal{A}^{H_0}[k] | \mathcal{S}^{H_0}[k])\\
\end{equation*}
\begin{equation*}
\mathcal{P}^{H_0}(S^{H_0}[k+1]|S^{H_0}[k],\mathcal{A}^{H_0}[k])) = \mathcal{P}^{H_0}(\tau^{H_0} | \pi^{H_0} ,\mathcal{M}^{H_0}). \qquad \\
\end{equation*}
$\mathcal{R}^0 (\tau^0)= \mathcal{R}^{H_0} (\tau^{H_0})$ follows directly from the definition of $\mathcal{R}^{H_0}$.
\section{Proof of Lemma\ref{L4}}
\vspace{-1mm}
In the lower-level MDP of the base module, we define:\\
\begin{equation*}
\quad \mathcal{A}^{L_0}[k]\doteq\{u[k] \}\\
\end{equation*}
\begin{equation*}
\mathcal{P}^{L_0} (\mathcal{S}^{L_0}[k+1]|\mathcal{S}^{L_0}[k],\mathcal{A}^{L_0}[k]) \doteq\mathcal{P}^0 ( \beta[k+1],\mathcal{S}^0[k+1] |(\beta[k],\\
\end{equation*}
\begin{equation*}
\mathcal{S}^0[k]),u[k])= \mathcal{P}^0 (\mathcal{S}^0[k+1]|\mathcal{S}^0[k],u[k]) \times\mathcal{P}^0 (\beta[k+1] | \mathcal{S}^0[k+1]\\
\end{equation*}
\begin{equation*}
,\beta[k]),
\end{equation*}
and
\begin{equation*}
\pi^{L_0} (\mathcal{A}^{L_0}[k] | \mathcal{S}^{L_0}[k]) \doteq \mathcal{P}^0 (u[k] | \mathcal{S}^0[k], \beta[k]),
\end{equation*}
therefore:
\begin{equation*}
\mathcal{P}^0(\tau^0 | \pi^0 ,\mathcal{M}^0) = \mathcal{P}^0 (\mathcal{S}^0[0])\mathcal{P}^0 (\beta[0]|\mathcal{S}^0[0]) \prod_{k= 0}^{K-1} (\mathcal{P}^0(u[k] | \mathcal{S}^0[k] \\
\end{equation*}
\begin{equation*}
 ,\beta[k])\mathcal{P}^0(\mathcal{S}^0[k+1]|\mathcal{S}^0[k],u[k]) \mathcal{P}^0(\beta[k+1]|\mathcal{S}^0[k+1],\beta[k]))
\end{equation*}
\begin{equation*}
= \mathcal{P}^{L_0} (S^{L_0}[0]) \prod_{k= 0}^{K-1} (\pi^{L_0} (\mathcal{A}^{L_0}[k] | \mathcal{S}^{L_0}[k]) \mathcal{P}^{L_0}(S^{L_0}[k+1]|S^{L_0}[k], \\
\end{equation*}
\begin{equation*}
\mathcal{A}^{L_0}[k]))= \mathcal{P}^{L_0}(\tau^{L_0} | \pi^{L_0} ,\mathcal{M}^{L_0}). \qquad \qquad \qquad \qquad \qquad \qquad \qquad
\end{equation*}
$\mathcal{R}^0 (\tau^0)= \mathcal{R}^{L_0} (\tau^{L_0})$ follows directly from the definition of $\mathcal{R}^{L_0}$.
\end{appendices}
\bibliographystyle{IEEEtran}


\section{Introduction}
\label{sec:intro}

Since 2008, the measurements of cosmic-ray (CR) positrons by PAMELA~\cite{PAMELA:2008gwm}, Fermi-LAT~\cite{Fermi-LAT:2011baq}, and AMS-02~\cite{AMS:2013fma, AMS:2019iwo} have shown an unexpected excess at energies $\gtrsim 10~\si{GeV}$.
Possible interpretations for this excess include annihilating/decaying dark matter~\cite{Bergstrom:2008gr, Cholis:2008hb, Yin:2008bs, Yuan:2013eja} and astrophysical sources like nearby pulsars within kpc~\cite{Hooper:2008kg, Yuksel:2008rf}.
In particular, the middle-aged pulsar Geminga with a distance of $\sim 250~\si{pc}$ is widely assumed to produce high energy positrons that could propagate to the Earth~\cite{Hooper:2008kg, Yuksel:2008rf, Yin:2013vaa, Feng:2015uta, Hooper:2017tkg, Cholis:2017ccs, Fang:2018qco, Profumo:2018fmz, Cholis:2018izy, Tang:2018wyr, Xi:2018mii, Johannesson:2019jlk, DiMauro:2019yvh, Fang:2019ayz, Manconi:2020ipm, Wang:2021xph, Fang:2022mdg}.

In 2017, the HAWC experiment observed $\sim 10~\si{TeV}$ $\gamma$ rays spatially extended about 2 degrees around Geminga, which would be produced by positrons and electrons of energies $\sim 100~\si{TeV}$ via inverse Compton scattering (ICS) off low energy photons~\cite{HAWC:2017kbo}.
Therefore, this observation confirms that Geminga is a source of high energy positrons and electrons.
But the surface brightness profile (SBP) measured by HAWC implies a diffusion coefficient smaller than the conventional value by at least two orders of magnitude.
The recent observation of another extended halo around the middle-aged pulsar J0621+3749 by
LHAASO further established the general conclusion of slow diffusion around pulsars
\cite{LHAASO:2021crt}.
Such slow diffusion results in much less CR positrons arriving at the Earth, unlikely to explain the positron excess.
Nonetheless, by assuming a two-zone diffusion model with slow diffusion in a small zone around the source but normal diffusion outside the zone, positrons originated from Geminga can still sufficiently contribute to the positron excess~\cite{Fang:2018qco,Profumo:2018fmz,Bao:2021hey}.

In addition, positrons and electrons from Geminga are also expected to induce extended ICS $\gamma$ rays in the energy range of Fermi-LAT.
Based on two-zone diffusion templates, an analysis of 10-yr Fermi-LAT $\gamma$-ray data by Xi et al.~\cite{Xi:2018mii} (denoted as X19 below) did not found such extended emission and derive a stringent constraint on the $\gamma$-ray flux in the $\sim 5\text{--}100$~\si{GeV} energy range.
According to this constraint and the HAWC data,  $e^\pm$ from Geminga with a single power-law injection spectrum can only contribute a small faction to the CR positron spectrum observed by AMS-02.

On the other hand, taking into account both a larger region of interest and the proper motion of the Geminga pulsar, another analysis of Fermi-LAT data by Di Mauro et al.~\cite{DiMauro:2019yvh} (denoted as D19 hereafter) claimed a discovery of extended $\gamma$-ray emissions around Geminga in the energy range of $\sim 10$--$100$~\si{GeV}.
However, considering both the corresponding $\gamma$-ray flux and the HAWC data, the Geminga contribution to the position flux they obtained is not enough for the AMS-02 excess.

Both the X19 and D19 analyses assumed a single power-law Geminga $e^\pm$ injection spectrum with a high energy cutoff.
The inconsistency with the AMS-02 data may indicate that there are less low energy positrons and electrons producing GeV $\gamma$ rays.
Therefore, we will attempt to modify the injection spectrum by adding a low energy cutoff, in order to simultaneously explain the HAWC, Fermi-LAT, and AMS-02 data.
The results of the $\gamma$-ray flux from the X19 and D19 analyses will be considered separately.

This paper is organized as follows.
In Section~\ref{sec:Geminga}, we describe the propagation of positrons and electrons produced by Geminga and the $\gamma$-ray flux induced by ICS.
In Section~\ref{sec:D19}, we simultaneously interpret the HAWC data, the Fermi-LAT $\gamma$-ray observation given by D19, and the AMS-02 positron spectrum assuming an $e^\pm$ injection spectrum with a low energy cutoff.
In Section~\ref{sec:X19}, we use the Fermi-LAT $\gamma$-ray constraint given by X19 to explore how much contribution Geminga can supply to the AMS-02 positron excess.
Section~\ref{sec:sum} gives the summary and discussion.

\section{Positrons and electrons from Geminga}
\label{sec:Geminga}

The Geminga pulsar is a $\gamma$-ray source discovered by SAS-2~\cite{Fichtel:1975}.
Its age is about $342~\si{kyr}$~\cite{Manchester:2004bp}, and the distance for the Earth is $250^{+120}_{-62}~\si{pc}$~\cite{Faherty:2007}.
Geminga is expected to emit lots of positrons and electrons, which diffuse away from Geminga and lose energies by upscatter low energy photons in the cosmic microwave background (CMB) and interstellar radiation backgrounds through ICS processes.

The propagation of CR $e^\pm$ is described by the diffusion-cooling
equation
\begin{equation}
\frac{\partial N}{\partial t} - \nabla\cdot(D\nabla N) -
\frac{\partial}{\partial
	E}(bN) = Q \,,
\label{eq:prop}
\end{equation}
where $N$ is the $e^\pm$ differential density, $E$ is the $e^\pm$ energy, $D$ is the diffusion coefficient, and $Q$ is the source term.
The energy loss rate $b$ includes both contributions from synchrotron radiation and ICS.
The synchrotron energy loss rate in a magnetic field $B$ is given by~\cite{Crusius:1988}
\begin{equation}
b_\mathrm{syn} = \frac{4 \sigma_\mathrm{T} \gamma_e^2 U_B}{3 m_e c},
\end{equation}
where $\sigma_\mathrm{T}$ is the Thomson cross section, $\gamma_e = E/(m_e c^2)$ is the $e^\pm$ Lorentz factor, and $U_B = B^2/(8\pi)$ is the energy density of the magnetic field.
The energy loss rate due to ICS is estimated following Ref.~\cite{Fang:2020dmi}.
We convert the propagation equation to a difference equation, which is solved using the numerical method described in Ref.~\cite{Fang:2018qco}.

We assume a spherically symmetrical two-zone diffusion scenario with the diffusion coefficient given by
\begin{equation}
D(E, r)=\left\{
\begin{aligned}
D_1(E), & & r< r_\star, \\
D_2(E), & & r\geq r_\star.
\end{aligned}
\right.
\label{eq:diff}
\end{equation}
Here $r$ is the distance from Geminga, and $r_\star$ denotes the boundary of the two diffusion zones.
Both $D_1(E)$ and $D_2(E)$ are assumed to have a form of $D_{100} (E/100~\si{TeV})^\delta$, where $D_{100}$ is the diffusion coefficient at $E = 100~\si{TeV}$, and $\delta = 0.33$ is adopted for a Kolmogorov-type diffusion~\cite{Kolmogorov:1941}.

The morphological SBP study of the extended TeV $\gamma$-ray emissions around Geminga by HAWC gives a diffusion coefficient $D_{100} = (3.2^{+1.4}_{-1.0})\times10^{27}~\si{cm^{2}~s^{-1}}$ for $100~\si{TeV}$ $e^\pm$ around Geminga, while a similar study of another nearby pulsar Monogen (PSR B0656+14) leads to $D_{100} = (15^{+49}_{-9})\times10^{27}~\si{cm^{2}~s^{-1}}$~\cite{HAWC:2017kbo}.
The joint fit of both results in $D_{100} = (4.5\pm 1.2)\times10^{27}~\si{cm^{2}~s^{-1}}$.
Thus, $D_{100}$ for the inner zone with $r< r_\star$ is at the order of $10^{27}~\si{cm^{2}~s^{-1}}$.
For the outer zone with $r \geq r_\star$, positrons and electrons propagate through the ordinary interstellar medium (ISM), and we take the GALPROP~\cite{Moskalenko:1997gh} default value $D_{100} = \num{1.7e30}~\si{cm^{2}~s^{-1}}$,
which is consistent with the measurements of CR secondary-to-primary ratios, particularly the $\mathrm{B}/\mathrm{C}$ ratio.

The source term for high energy $e^\pm$ injected by Geminga is assumed as
\begin{equation}
Q(t,E,r) = q(t,E)\delta(r)\,,
\label{eq:source}
\end{equation}
where
\begin{equation}
q(t,E)=q_0\left(1+\frac{t}{\tau}\right)^{-2}E^{-\gamma} \exp\left(-\frac{E}{E_\mathrm{hc}}\right) \exp\left(-\frac{E_\mathrm{lc}}{E}\right).
\label{eq:time_profile}
\end{equation}
$\tau$ is the characteristic initial spin-down time scale of the Geminga pulsar, taken to be $12~\si{kyr}$ following Ref.~\cite{HAWC:2017kbo}.
$\gamma$ is the injection spectral index for $e^\pm$.
$E_\mathrm{hc}$ and $E_\mathrm{lc}$ are the high and low energy cutoffs, respectively.
$q_0$ is a constant determined by the normalization relation
\begin{equation}
\int_{E_{\rm min}}^{E_{\rm max}}q(t_\mathrm{s},E)EdE=\eta\dot{E}_s\,,
\label{eq:norm}
\end{equation}
with the Geminga age $t_\mathrm{s}=342~\si{kyr}$ and the spin-down luminosity $\dot{E}_\mathrm{s}=\num{3.2e34}~\si{erg~s^{-1}}$~\cite{Manchester:2004bp}.
Here $\eta$ is the conversion efficiency for the spin-down energy converted to $e^\pm$ energies.
We will not consider the difference between the positrons and electrons when calculating the $\gamma$-ray flux, and the positron flux $\Phi_{e^+}$ is just a half of the total $e^\pm$ flux.

The photon emissivity due to $e^\pm$ ICS based on the Klein-Nishina cross section is given by~\cite{Fang:2007sc}
\begin{equation}
Q_\mathrm{ICS} (t, E_\gamma, r) = 4\pi \sum_j \int^{\infty}_{0}d\epsilon
\, n_j(\epsilon)\int^{E_{\mathrm{max}}}_{E_{\mathrm{min}}}dE \,
J(t, E, r)F(\epsilon, E_{\gamma}, E),
\end{equation}
$n_j(\epsilon)$ is the number density of a background photon component $j$ with energy $\epsilon$, temperature $T_j$, and energy density $U_j$, expressed as
\begin{equation}
n_j(\epsilon)=\frac {15U_j}{(\pi k
	T_j)^4}\frac{\epsilon^2}{\exp(\epsilon/kT_j)-1},
\end{equation}
where $k$ is the Boltzmann constant.
The $e^\pm$ energy threshold for upscattering a target photon with energy $\epsilon$ to a photon with energy $E_\gamma$ is
\begin{equation}
E_\mathrm{min} = \frac{1}{2} \left( E_\gamma + \sqrt{E_\gamma^2 + \frac{E_\gamma m_e^2 c^4 }{\epsilon}} \right).
\end{equation}
$J(t, E, r) = v_e N(t, E, r)/(4\pi)$ is the $e^\pm$ intensity, with $v_e$ denoting the $e^\pm$ speed, which approaches the light speed $c$ for high energy $e^\pm$.
The function $F$ is given by
\begin{equation}
F(\epsilon, E_{\gamma}, E) = \frac{3\sigma_\mathrm{T}}{4\gamma_e^2 \epsilon}
\left[ 2q\ln q + (1+2q)(1-q) + \frac{\Gamma^2 q^2(1-q)}{2(1+\Gamma q)}
\right],
\end{equation}
with
\begin{equation}
\Gamma = \frac{4\epsilon\gamma_e}{m_e c^2},\quad
q = \frac{E_\gamma}{\Gamma (E_e - E_\gamma)}.
\end{equation}

Following Ref.~\cite{HAWC:2017kbo}, we consider three background photon components, including the CMB, the IR background, and the optical background, for the ICS processes.
The temperatures and energy densities are presented in Table~\ref{tab:ISRF}.
Integrating $Q_\mathrm{ICS} (t_s, E_\gamma, r)$ along the light of sight~\cite{Liu:2019sfl}, we obtain the $\gamma$-ray flux for specific energy $E_\gamma$ and angular separation $\theta$,
\begin{equation}
I(E_\gamma,\theta)=\frac{1}{4\pi}\int_{l_\mathrm{min}}^{l_\mathrm{max}}dl \,Q_\mathrm{ICS}(t_s, E_\gamma, r).
\end{equation}
Then we integrate out $\theta$ to get the energy spectrum of the $\gamma$-ray flux $\Phi_\gamma$, or  integrate out $E_\gamma$ to derive the SBP as a function of $\theta$.

\begin{table}[!t]
\begin{center}
\setlength{\tabcolsep}{1em}
\renewcommand\arraystretch{1.3}
\caption{Temperature $T_j$ and energy density $U_j$ of three background photon components~\cite{HAWC:2017kbo}.}
\label{tab:ISRF}
\vspace{1em}
\begin{tabular}{ccc}
\hline\hline
Component $j$ & $T_j$ (K) & $U_j$ ($\si{eV/cm^3}$)\\
\hline
CMB& 2.7 & 0.26\\
IR & 20 & 0.3\\
Optical & 5000 & 0.3\\
\hline\hline
\end{tabular}
\end{center}
\end{table}

\section{Result according to the D19 gamma-ray observation}
\label{sec:D19}

In this section, we try to interpret the HAWC and AMS-02 data according to the Fermi-LAT $\gamma$-ray observation from the D19 analysis~\cite{DiMauro:2019yvh}.
Both the results without and with the low energy cutoff $E_\mathrm{lc}$ in the $e^\pm$ injection spectrum are calculated for comparison.

\begin{figure}[!t]
\centering
\subfigure[~$\gamma$-ray spectrum \label{fig:DM:gamma}]{\includegraphics[width=0.48\textwidth]{DM_flux.pdf}}
\subfigure[~Surface brightness profile~\label{fig:DM:SPB}]{\includegraphics[width=0.465\textwidth]{DM_sbp.pdf}}
\subfigure[~CR positron spectrum\label{fig:DM:posi}]{\includegraphics[width=0.48\textwidth]{DM_positron.pdf}}
\caption{The $\gamma$-ray spectrum around Geminga (a), the Geminga SBP (b), and the CR positron spectrum (c) assuming $e^\pm$ injection spectra without a low energy cutoff for $\eta = 0.8$ (blue solid lines) and with $E_\mathrm{lc} = 20~\si{GeV}$ for $\eta=0.4$ (red dashed lines).
In the upper left panel, the green region denotes the $\sim 10~\si{TeV}$ spectral data measured by HAWC~\cite{HAWC:2017kbo}, and the data points and upper limits in $10~\si{GeV} \lesssim E_\gamma \lesssim \si{TeV}$ are given by the D19 analysis of Fermi-LAT data~\cite{DiMauro:2019yvh}.
The data points in the upper right panel shows the HAWC observation of the Geminga SBP~\cite{HAWC:2017kbo}.
The lower panel displays the positron spectrum measured by AMS-02~\cite{AMS:2019iwo}.}
\label{fig:DM}
\end{figure}

Firstly, we consider an $e^\pm$ injection spectrum without $E_\mathrm{lc}$, and adjusting the energy conversion efficiency $\eta$ to meet the data.
Setting the boundary radius $r_\star = 50~\si{pc}$, the diffusion coefficient at $E = 100~\si{GeV}$ in the inner diffusion zone $D_{100} = 4.5\times10^{27}~\si{cm^{2}~s^{-1}}$, the ISM magnetic field $B = 3~\si{\mu G}$, the $e^\pm$ injection spectral index $\gamma = 2.25$, and the high energy cutoff $E_\mathrm{hc} = 511~\si{TeV}$, we derive the $\gamma$-ray spectrum around Geminga, the Geminga SBP, and the CR positron spectrum at the Earth for $\eta = 0.8$, shown as the blue solid lines in Fig.~\ref{fig:DM}.

In order to compare the predictions and the observations, we show the $\sim 10~\si{TeV}$ spectral data measured by HAWC~\cite{HAWC:2017kbo} and the Fermi-LAT data points and upper limits from $\sim 10~\si{GeV}$ to $\sim \si{TeV}$ given by the D19 analysis~\cite{DiMauro:2019yvh} in Fig.~\ref{fig:DM:gamma}.
The HAWC observation of the Geminga SBP~\cite{HAWC:2017kbo} is demonstrated in Fig.~\ref{fig:DM:SPB}, while the positron spectrum measured by AMS-02~\cite{AMS:2019iwo} is displayed in Fig.~\ref{fig:DM:posi}.
For the above setup with $\eta = 0.8$, we find that the $\gamma$-ray prediction can well interpret the $\gamma$-ray spectrum and the SBP, and the predicted $e^+$ spectrum can explain the AMS-02 data at $E \gtrsim 100~\si{GeV}$.
However, a $80\%$ efficiency of the spin-down energy converted to $e^\pm$ energies looks unrealistic.

Secondly, we introduce a low energy cutoff $E_\mathrm{lc} = 20~\si{GeV}$ in the $e^\pm$ injection spectrum with other parameters unchanged, and find that the observational data can be explained for $\eta = 0.4$, as illustrated as the red dashed lines in Fig.~\ref{fig:DM}.
Such a $40\%$ conversion efficiency is much more reasonable than the previous one.
Now the predicted positron flux at $E \lesssim 100~\si{GeV}$ seems slightly lower than the blue solid line, but we can still interpret the AMS-02 data at $E \gtrsim 100~\si{GeV}$ very well.

\section{Result according to the X19 gamma-ray constraint}
\label{sec:X19}

In contrast to the D19 analysis~\cite{DiMauro:2019yvh}, the X19 analysis of the Fermi-LAT data have not found any extended $\gamma$-ray emission around Geminga, deriving a rather stringent constraint on the $\gamma$-ray flux at $\sim 5\text{--}100~\si{GeV}$~\cite{Xi:2018mii}.
In this section, we consider this constraint to see how it affects the Geminga contribution to the CR positron spectrum, assuming a low energy cutoff in the $e^\pm$ injection spectrum.
However,  we find it impossible to simultaneously explain the HAWC, Fermi-LAT, and AMS-02 data, because the X19 constraint is too strict.
Instead, we would like to know how much contribution Geminga can provide to the AMS-02 positron excess.

For this purpose, we treat $\gamma$, $E_\mathrm{hc}$, $E_\mathrm{lc}$, $\eta$, $B$, and $D_{100}$ in the inner diffusion zone as free parameters and perform a scan in the parameter space with fixed $r_\star$, utilizing the \texttt{MultiNest} algorithm~\cite{Feroz:2008xx} to improve the fitting efficiency.
The ranges for the free parameters in the scan are chosen to be
\begin{eqnarray}
&& 1.8 < \gamma < 2.2,\quad
200~\si{TeV} < E_\mathrm{hc} < 600~\si{TeV},\quad
100~\si{GeV} < E_\mathrm{lc} < 900~\si{GeV},
\nonumber\\
&& 0.1 < \eta < 0.4,\quad
3~\si{\mu G} < B < 8~\si{\mu G},\quad
10^{26}~\si{cm^2~s^{-1}} < D_{100} < 10^{27}~\si{cm^2~s^{-1}}.
\end{eqnarray}

\begin{table}[!t]
\begin{center}
\setlength{\tabcolsep}{1em}
\renewcommand\arraystretch{1.3}
\caption{Parameters in the best results for fixed $r_\star$.}
\label{tab:bestpara}
\vspace{1em}
\begin{tabular}{cccc}
\hline\hline
$r_\star$ (pc) & $50$ & $70$ & $100$ \\
\hline
$\gamma$ & $1.94$ & $1.93$ & $1.87$\\
$E_\mathrm{hc}$ (TeV) & $520$ & $537$ & $463$\\
$E_\mathrm{lc}$ (GeV) & $870$ & $302$ & $547$ \\
$\eta$ & $0.12$ & $0.21$ & $0.20$ \\
$B$ ($\si{\mu G}$) & $5.0$ & $6.9$ & $7.2$\\
$D_{100}$ ($10^{27}~\si{cm^2~s^{-1}}$) & $4.8$ & $7.8$ & $8.5$\\
\hline\hline
\end{tabular}
\end{center}
\end{table}

In order to get optimistic results,
we adopt the most loose upper limits on the $\gamma$-ray flux in $10\text{--}500~\si{GeV}$ derived by the X19 analysis, i.e., the upper limits in the upper panel of Fig.~6 in the X19 paper~\cite{Xi:2018mii}.
The parameters of the best results we obtain for $r_\star = 50,~70,~100~\si{pc}$ are listed in Table~\ref{tab:bestpara}.
The corresponding predictions for the $\gamma$-ray spectrum, the SBP, and the positron spectrum are demonstrated in Fig.~\ref{fig:Xi}.
While the HAWC data are properly fitted and the $\gamma$-ray flux in $5~\si{GeV} \lesssim E_\gamma \lesssim 100~\si{GeV}$ lies below the X19 upper limits, we find that Geminga can only supply less than $36\%$ of the AMS-02 positron flux at $E \sim 300\text{--}600~\si{GeV}$.

\begin{figure}[!t]
\centering
\subfigure[~$\gamma$-ray spectrum \label{fig:Xi:gamma}]{\includegraphics[width=0.48\textwidth]{Xi_flux.pdf}}
\subfigure[~Surface brightness profile~\label{fig:Xi:SPB}]{\includegraphics[width=0.465\textwidth]{Xi_sbp.pdf}}
\subfigure[~CR positron spectrum\label{fig:Xi:posi}]{\includegraphics[width=0.48\textwidth]{Xi_positron.pdf}}
\caption{The best results for the $\gamma$-ray spectrum around Geminga (a), the Geminga SBP (b), and the CR positron spectrum (c) assuming $r_\star = 50~\si{pc}$ (blue solid lines), $r_\star = 70~\si{pc}$ (red dashed lines), and $r_\star = 100~\si{pc}$ (purple dotted lines).
In the upper left panel, the upper limits in $5~\si{GeV} \lesssim E_\gamma \lesssim 100~\si{GeV}$ are given by the X19 analysis of Fermi-LAT data~\cite{Xi:2018mii}.
The other experimental data are the same as in Fig.~\ref{fig:DM}.}
\label{fig:Xi}
\end{figure}

These results show that the X19 constraint favor $\gamma < 2$, $\eta \lesssim 0.21$, and $E_\mathrm{lc}$ of several hundred GeV, which suppress the $\gamma$-ray flux at $\sim \mathcal{O}(10)~\si{GeV}$.
According to an approximate relation~\cite{Xi:2018mii}
\begin{equation}
E_\gamma = 20 \left(\frac{E}{100~\si{TeV}}\right)^2 \si{TeV}
\end{equation}
for $e^\pm$ ICS off CMB photons, $\mathcal{O}(10)~\si{GeV}$ $\gamma$ rays are induced by $\mathcal{O}(\si{TeV})$ $e^\pm$.
Thus, the X19 constraint implies less $\mathcal{O}(\si{TeV})$ positrons and electrons from Geminga, resulting in a lower CR positron flux for $E \sim \mathcal{O}(100)~\si{GeV}$ at the Earth, which is insufficient to explain the AMS-02 excess.

\section{Summary and discussion}
\label{sec:sum}

In this work, we attempt to explain the AMS-02 positron excess by the nearby pulsar Geminga assuming a two-zone diffusion scenario and an $e^\pm$ injection spectrum with a low energy cutoff, taking into account the $\gamma$-ray data from HAWC and Fermi-LAT.
The analyses of Fermi-LAT data for extended $\gamma$-ray emissions around Geminga by two groups have obtained different results.
While the X19 analysis found no such emission and derive upper limits on the $\gamma$-ray flux, the D19 analysis claimed an observation of the extended $\gamma$ rays.
We have considered both results separately.

By fitting the D19 observation and the HAWC data assuming no low energy cutoff in the injection spectrum, we find that the conversion efficiency $\eta$ should be as large as $80\%$ to account for the AMS-02 positron excess.
Nonetheless, if a low energy cutoff $E_\mathrm{lc} = 20~\si{GeV}$ is introduced, we would only need a $40\%$ conversion efficiency, which is much more realistic.
Therefore, it is plausible to interpret the positron excess by Geminga, according to the D19 analysis.

On the other hand, if the stringent constraint from the X19 analysis is considered, we find that Geminga could not accounting for the total positron excess.
We carry out a scan in the parameter space for the boundary radius $r_\star = 50,~70,~100~\si{pc}$ and require to fit the HAWC data and satisfy the X19 constraint.
The best results we obtain can only explain a faction of the AMS-02 positron flux lower than $36\%$ at $E\sim \mathcal{O}(100)~\si{GeV}$.
This may imply that more nearby pulsars or other sources are needed to interpret the positron excess.

Since the different conclusions obtained above come from the contradictory results of the two Fermi-LAT analyses, it is crucial to know whether result is true.
This may require a more careful data analysis and more Fermi-LAT data.

\begin{acknowledgments}

We thank Kun Fang for providing the code to solve the two-zone propagation equation.
This work is supported in part by the National Natural Science Foundation of China under Grants No.~11875327 and No.~11805288,
the Fundamental Research Funds for the Central Universities,
and the Sun Yat-Sen University Science Foundation.
Q.Y. is supported by the Program for Innovative Talents and Entrepreneur in Jiangsu.

\end{acknowledgments}

\bibliographystyle{utphys}