Joint Offloading and Resource Allocation Based on Lyapunov Algorithm in Delay-Sensitive SAGIN

I. INTRODUCTION

AS a new network architecture, space-air-ground integrated network (SAGIN) can break through the limitations of the ground to realize wide-area coverage, high-speed transmission, and heterogeneous interconnections [1]–[3]. According to the ITU mobile telecommunication system (IMT- 2030) promotion group, the seamless global coverage through SAGIN can meet the demand for safe and reliable humanmachine- object unlimited connectivity [4], and it is one of the main directions of the research on the 6G network architecture in the future. In recent years, more and more organizations have begun to carry out research in the field of SAGIN, such as Telesat [5], OneWeb [6], SpaceX [7], and spawned emerging business scenarios, for instance, high-precision satellite positioning, wide-area Internet of things (IoT), and full-scenario autonomous driving.

Meanwhile, with the development of various up and coming applications, the requirements for low latency are also increasing. Simultaneously, the SAGIN, characterized by its extensive scale and intricate architecture, faces challenges related to high latency and high usage costs within satellite networks, which significantly impact the application of SAGIN. mobile edge computing (MEC) emerges as a pivotal approach to augmenting device computing capabilities while mitigating latency. By relocating computational tasks to mobile edge nodes, MEC achieves several benefits, including decreased network transmission delays, reduced energy consumption and minimized data jitter [7]–[9]. Therefore, introducing MEC into SAGIN can avoid long-distance transmission and congestion by making it possible that the computational tasks of ground users can be directly offloaded to satellites for processing, without the need to offload the tasks to cloud data centers via satellite backhaul links [10], [11]. However, due to the highly dynamic and heterogeneous nature of SAGIN, the following challenges must be addressed in combining SAGIN with MEC: 1) firstly, the high-speed network nodes of SAGIN lead to frequent switching; 2) in addition, the introduction of MEC raises the complexity of the network, and adds a new dimension in terms of optimization choices [12], [13].

In recent years, researchers have made significant progress in the field of SAGIN. The satellite relay transmission scheme proposed in the [14]–[16] effectively utilizes the powerful centralized arithmetic resources of cloud computing centers by routing the computational tasks of the ground devices to be processed in remote cloud data centers, while the feasibility of task processing directly on the satellite is ignored. The core idea of introducing MEC technology in SAGIN is to sink the rich computing and caching resources of cloud data centers to the edge of the network closer to mobile users, thus effectively reducing the latency of task processing and thus improving the quality of user experience.

For such purpose, the resource management decision of SAGIN for IoT scenarios with only one LET satellite in consideration was proposed [17], [18]. In which the unmanned aerial vehicles (UAV) act as “flying gateways” to collect computing tasks from ground IoT devices and decide whether perform local edge computing or forward tasks to the satellite. In this network, the offloading policy of the UAV was optimized by a linear programming approach and a deep reinforcement learning approach, respectively. In [19], an air-space-assisted hybrid cloud edge computing system was proposed, in which UAVs were equipped with MEC servers to provide low-latency edge computing services, while a low Earth orbit (LEO) satellite was equipped with cloud servers to provide ubiquitous cloud computing services for IoT devices. The article presented a fair-aware resource scheduling problem to minimize the maximum task execution delay among IoT devices.

SAGIN with multiple satellites forming a satellite network has also been studied by some scholars. In [20], an integrated geosynchronous Earth orbit (GEO) satellites space-airsea network architecture with edge computing capability was proposed to provide flexible computing services for maritime communications. This paper considered minimizing the total task execution delay under the constraints of communication rate and computational power, and designed a deep reinforcement learning-based solution to solve this complex optimization problem. [21] presented a satellite constellation design problem in a satellite Internet of things (SIoT) scenario. The paper used space-time graph (STG) to describe the network topology and then used various Tabu Search (TS) algorithms to obtain the satellite constellation configuration with the best QoS level.

It can be seen that current research on MEC and SAGIN lacks consideration of air-based networks such as high-altitude platforms, and generally pays less attention to impact of the periodic motion of satellite constellations on system performance. To address the above problems, this paper proposes a SAGIN architecture with LEO satellites carrying MEC servers and using high altitude balloons (HAB) as fixed airborne relays. This design fully accounts for the impact of satellite constellation motion. Additionally, we introduce the Lyapunov optimization method to address optimization challenges over a long-term time horizon. The main contributions of this paper are summarized as follows:

· Aiming at the requirements of SAGIN users on service quality, this paper proposes to optimize user delay as well as usage cost in the lower layer network composed of HABs and users. To accommodate the highly dynamic nature of the satellite constellation, the problem is modeled as a long-term time-averaged optimization problem.

· Responding to the diffculty of solving long-term timeaveraged problems, Lyapunov’s theory is utilized to decouple the coupling relationship of each time slot, simplifying the problem form, and methods such as logsum- exp smoothing function as well as successive convex approximation (SCA) algorithm are used to achieve the solution of each sub-optimization problem.

The remainder of the paper is organized as follows. Section II describes the system model and the construction of the optimization problem. In Section III, we resolve the optimization problem based on the Lyapunov optimization algorithm and perform the performance analysis. Section IV presents the simulation results. Finally, the conclusion and the future work in Section V.

II. SYSTEM MODEL

A. SAGIN Architecture

The SAGIN consists of the following three parts: the satellite constellation consists of S LEO satellites carrying MEC servers with the set [TeX:] $$\mathcal{S}=\{1,2, \cdots, S\},$$ and there exist laser inter satellite links (ISL) between the satellites for data diversion and computational task continuation; the set of G HABs denoted as [TeX:] $$\mathcal{G}=\{1,2, \cdots, G\},$$ as an air relay for data forwarding; the set of U ground users is denoted as [TeX:] $$\mathcal{U}=\{1,2, \cdots, U\},$$ and the ground users forward the computation tasks to the LEO satellite for edge computation through the HABs and obtain the computation results in the reverse link.

The high speed movement of satellites in SAGIN makes the network topology change over time, so we divide the optimization process into T time slots of [TeX:] $$\mathcal{T}=\{1,2, \cdots, {T-1}\}.$$ Here, we consider the data transmission in a unit latitude and longitude grid and assume that all ground users in the grid can establish connections with any HAB overhead, but users can only establish connections with at most one HAB in a single time slot. The system model is shown in Fig. 1.

B. SAGIN Link Model

We assume that HABs and ground users communicate with each other over the C-band using OFDMA, and that HABs allocate the available bandwidth equally to the users connected. Since HABs can take advantage of the high wind speed in the stratosphere to adjust their positions, the relative motion between ground users and HABs is negligible [22]. In addition, since there is a LOS link between the HABs and the ground users, it is assumed that only free-space path loss as well as rain attenuation are affected in the transmission link. Therefore, according to Shannon’s formula, the uplink transmission capacity from the ground user u to the HAB g can be obtained as

where [TeX:] $$x_{u, g} \in\{0,1\},$$ is the indicator variable between the user u and the HAB g, with 1 denoting the connected state and 0 denoting the disconnected state; [TeX:] $$W_{bool}$$ means the transmission bandwidth between the user u and the HAB [TeX:] $$g ; p_u$$ represents the transmit power of the user [TeX:] $$u ; \mathbb{E}\left\{I_{u, g}\right\}$$ denotes the expectation of the interference that this link is subject to; [TeX:] $$n_0$$ indicates the the standard deviation of additive white Gaussiannoise; [TeX:] $$G_{u,g}$$ expresses the channel gain coefficient between the user u and the HAB g, which can be given by

where [TeX:] $$G_T \text{ and } G_R$$ are the transmit and receive antenna gains, and [TeX:] $$L_{\text{rain}}$$ means the rain attenuation. [TeX:] $$P L_{u, g}$$ represents the free space path loss between the user u and the HAB g, which can be given by

where λ indicates the wavelength and [TeX:] $$d_{u, g}$$ is the distance between the user u and the HAB g. The total transmission rate of the user u in the slot t is

In terms of interference, there is no inter-user interference within the service area of a single HAB because the HABs allocate the available bandwidth equally to the ground users connected to them. However, the number of channels for HABs is limited, and although there is no interference between users connected to a single HAB, but terrestrial users connected to different HABs may reuse the same sub-channel, thus generating inter-user interference between the service ranges of different HABs. Therefore, the average inter-user interference received by the link between the terrestrial user u and the HAB g can be given by

where [TeX:] $$\mathcal{U}(g)$$ denotes the set of ground users connected to the HAB [TeX:] $$g, P_I$$ represents the probability of other ground users interfering with the current user’s air-ground link, numerically equal to the reciprocal of the number of ground users served by the HAB g.

C. SAGIN Computational Offloading Model

Firstly, three task attributes [TeX:] $$\left\{C_u, f_u, O_u\right\}$$ contained in the user u are defined, which are task computational complexity [TeX:] $$C_u$$ i.e., the number of CPU cycles required to compute 1bit data, local computational capability [TeX:] $$f_u$$ i.e., the local CPU clock frequency of the user, and task computation result [TeX:] $$O_u$$ i.e., the amount of data returned by the computational task. In order to make full use of the local computing capabilities of ground users, we adopt a partial offloading method and define the computation task offloading ratio of ground users in slot t as [TeX:] $$\boldsymbol{\eta}(t)=\eta_1(t), \eta_2(t), \cdots, \eta_U(t), \text { and } \eta_u(t) \in[0,1] .$$ In the following, local and edge computing are modeled separately.

1) Local Computation

The local computation process does not go through data transmission, so only the computation delay is included, and the local computation time of the user u in the slot t can be expressed as

where [TeX:] $$D a_u(t)$$ is the total task data volume to be processed by the user u in the slot t.

2) Edge Computation

The edge computation duration of the user u in the slot t can be expressed as

where [TeX:] $$\mathbb{E}\left\{d_{g, s}\right\}$$ denotes the average distance from the HAB g to the multiple LEO satellites it is connected to. Considering that the computational power of the MEC server is much larger than that of the ground user’s local computation, and the amount of data in the computation results of the task is much smaller than that of the task itself, the computation delay and the transmission delay of the computation results can be ignored compared to the propagation delay.

Since the local computation and edge offloading computation are carried out simultaneously, the total latency of the ground user u in the slot t can be given by

In addition to the latency, the high usage cost is also an important factor for limiting SAGIN, for which the following usage cost of the user u in the slot t is given by

The user-side cost [TeX:] $$ta_{UE}$$ is the power consumption of the user to consume a unit of power, which consists of the transmit power and the computational power that is proportional to the third power of the clock frequency, and the effective switching capacitance [TeX:] $$\kappa=10^{-11}$$ is set in this paper [23]. In addition, the x² function is adopted to describe the user fairness during the optimization process, which prevents the users with better channel conditions from blindly increasing their transmit power, and at the same time, enables the users with poor channel conditions to obtain satisfactory throughputs at a lower cost. The cost required by a ground user to connect to a single HAB is defined as [TeX:] $$ta_{HAB}.$$ the cost required to use a single LEO satellite is [TeX:] $$t a_{\text {sat }}, \text { and } \Phi(g)$$ denotes the number of LEO satellites to which the HAB g is connected.

D. Data Cache Queue Model

Considering that the computation tasks of ground users arrive randomly, this paper constructs a data caching queue model for each ground user. Define the user-side queue vector as [TeX:] $$\boldsymbol{Q}(t)=\left\{Q_1(t), \cdots, Q_U(t)\right\},$$ and assume that the task data arrivals of the ground users obey the Poisson distribution with intensity [TeX:] $$\nu,$$ which satisfies the independent homogeneous distribution in each time slot, so the task data arrival vector of all users in slot t is [TeX:] $$\boldsymbol{a}(t)=\left\{a_1(t), \cdots, a_U(t)\right\} .$$ Similarly, the data departure vector for each ground user in the slot t is [TeX:] $$\boldsymbol{D} \boldsymbol{a}(t)=\left\{D a_1(t), \cdots, D a_U(t)\right\} .$$ In summary, the queue update formula for the ground user u in the slot t can be obtained as follows

E. Optimization Problem Modeling

In this paper, we optimize the performance in terms of both latency and usage cost from the QoS metrics that users care about. We propose the QoS utility function of the user u in the slot t with a weighting factor ω

From the definitions of (8) and (9), we can see that the network latency and usage cost are related to the user’s transmit power and local computational capability. The greater the user’s transmit power and local computational capability, the smaller the user’s local computing and edge computing latency, but at the same time the user’s usage cost is also greater. Therefore, it can be seen that the optimization of latency and usage cost has an antagonistic relationship.

In the following, we take minimizing the long-term timeaveraged QoS utility function of the ground users as the optimization objective, and builds the following optimization problem with the ground user’s offloading decision (offloading direction [TeX:] $$\boldsymbol{x}(t)$$, amount of processed data [TeX:] $$\boldsymbol{D} \boldsymbol{a}(t),$$ offloading ratio [TeX:] $$\boldsymbol{\eta}(t)$$) and the ground user’s transmit power [TeX:] $$\boldsymbol{p}(t)$$ as the decision variables

The constraint C1 is a long-term time-averaged throughput constraint which is set to ensure that each ground user can obtain sufficient edge computing resources, and [TeX:] $$R_{th}$$ represents the throughput threshold value; the constraint C2 is the matching variable constraint between HABs and ground users, ensuring that a single ground user can only connect to at most one HAB in the slot t; the constraint C3 represents the stability constraint of the user data cache queue, which ensures that the data queue backlog of each user will not increase infinitely and cause huge queuing delay; the constraint C4 denotes the task computing delay constraint, which ensures that the computing tasks dispatched in the current slot can be completed within the current slot; the constraints C5–C8 are the value range constraints of the user’s transmit power, utility function weight factor, matching variable and task offloading ratio.

The optimization problem [TeX:] $$\mathbf{P1},$$ as a long-term time averaging problem, is difficult to solve due to the existence of cache queues, the system state and the network decision between each time slot coupled with each other. Therefore, we propose to use the Lyapunov optimization method to solve the stochastic network optimization problem such as [TeX:] $$\mathbf{P1}.$$

III. RESOURCE ALLOCATION BASED ON LYAPUNOV OPTIMIZATION

In this section, we solve the above proposed problem of minimizing the long-term time-averaged QoS utility function based on Lyapunov optimization theory. In the following, the original optimization problem [TeX:] $$\mathbf{P1}$$ is first transformed to eliminate the long-term time-averaging constraint in it, and then the problem is decoupled and the decomposed problem is solved in each time slot in turn.

A. Format Transformation of The Optimization Problem

The principle of Lyapunov optimization is to make separate decisions in each time slice to achieve optimal long-term timeaveraging performance and to ensure the stability of the system queue state, which requires that the random events and the system state observed in each slot t must be independently and identically distributed so that the decoupling of each slot can be performed. However, due to the existence of the longterm time averaging constraint C1, the system states of each time slot are coupled with each other, and in order to satisfy the long-term time averaging constraint, the decision of each time slot must consider the impact brought by the decision of the previous time slot, and also facilitate the decision of the subsequent time slot. Therefore, [TeX:] $$\mathbf{P1}$$ cannot directly apply the stochastic network optimization algorithm based on Lyapunov theory, but needs to use virtual queues to transform the long-term time-averaging constraint into queue stability constraints, thus spreading the limits of the long-term timeaveraging constraint to each slot and facilitating subsequent decoupling. Therefore, a set of virtual queues is defined for constraint [TeX:] $$\mathrm{C} 1: \boldsymbol{S}(t)=\left\{S_1(t), \cdots, S_U(t)\right\},$$ and define the virtual queue update formula for the ground user u in the slot t as follows

where the initial value of the queue backlog is 0, i.e., [TeX:] $$S_u(0)=0, \forall u \in \mathcal{U} .$$ It can be seen that in the virtual queue [TeX:] $$S_u(t),$$ the queue arrivals per time slot are the throughput threshold [TeX:] $$R_{\mathrm{th}}$$ and the amount of data leaving the queue in the slot t is the throughput of the user u in that time slot.

According to Neely’s stochastic network optimisation theory [24], the original long term time averaging constraint C1 can be ensured to hold by guaranteeing the stability of the virtual queue defined by (13), i.e.,

At this point, problem [TeX:] $$\mathbf{P1}$$ is transformed into [TeX:] $$\mathbf{P2}$$, which can be decoupled and computed using the Lyapunov optimization algorithm.

B. Problem Decoupling Based on Lyapunov Theory

According to the above analysis, the system queue backlog vector for slot t is defined as [TeX:] $$\boldsymbol{\Theta}(t)=\{\boldsymbol{Q}(t), \boldsymbol{S}(t)\},$$ which consists of two parts: the real data cache queue as well as the virtual queue. According to [24], the following quadratic Lyapunov function is given by

Further define the conditional Lyapunov drift function in a single time slot as follows

and the conditional Lyapunov drift-plus-penalty function in a single time slot as follows

where the penalty factor V is used to adjust the preference of the algorithm for the original optimization objective and the stability of the queues. Since the core principle of the Lyapunov algorithm is minimizing the upper bound of the driftplus- penalty function to ensure queue stability and achieving minimization of the original optimization objective, so we give the following expression of the upper bound of the Lyapunov drift-plus-penalty function.

Theorem 1. The Lyapunov drift-plus-penalty function in (18) satisfies the following upper bound constraint

where B is a constant with the following specific expression

Proof. Please see Appendix A.

where [TeX:] $$a_{\max }(t), D a_{\max }(t), R_{\max }$$ represent the upper bound of data arrivals, data departures and user transmission rate per time slot, respectively. The idea of the solution of the Lyapunov algorithm is to ensure the stability of the queue and at the same time achieve the minimization of the original optimization objective by minimizing the upper bound of the drift plus penalty function. Therefore, the long-term time-averaging optimization objective can be minimized by observing the current system queue backlog, data arrivals, and channel state parameters in each slot and using them to select an appropriate transmission strategy. Thus, the stochastic optimization problem [TeX:] $$\mathbf{P2}$$ is transformed into a series of individual per-slot conventional optimization problems as [TeX:] $$\mathbf{P3}$$.

It can be seen that problem [TeX:] $$\mathbf{P3}$$ is a mixed-integer nonconvex min-max problem, and the optimization variables are coupled with each other, so we consider dividing problem [TeX:] $$\mathbf{P3}$$ into multiple subproblems according to the optimization variables and then adopting block coodinate descent (BCD) for alternate optimization.

C. Problem Solving within A Single Slot

1) Optimization of association variables between HABs and ground users.

The problem [TeX:] $$\mathbf{P3}$$ can be obtained by fixing the [TeX:] $$p_u, D a_u \text{ and } \eta_u$$ as follow

[TeX:] $$\mathbf{SP1}$$ is an integer planning problem due to the presence of the binary matching variable x(t). For this reason, it is considered to relax x(t) into the interval from 0 to 1, i.e., the constraint C7 is transformed into [TeX:] $$x_{u, g}(t) \in[0,1],$$ after which an inequality constraint [TeX:] $$x_{u, g}(t)\left[x_{u, g}(t)-1\right] \geq 0$$ is introduced to control [TeX:] $$x_{u, g}(t)$$ can only take the value of 0 or 1, thus transforming the integer planning problem into a continuous optimization problem. Also from (9), the delay expression within a single time slot is in the form of max[x, y], which makes the problem difficult to solve as a min-max nested problem, so here a log-sum-exp smoothing function is introduced to transform the delay expression in the form of max[x, y] [24]

By observing the expressions of [TeX:] $$D e_u^L(t) \text { and } D e_u^E(t)$$ given by (6) and (7), respectively, we can see that at slot t, when [TeX:] $$p(t), D a(t), \text { and } \eta(t)$$ are fixed, the local computation latency [TeX:] $$D e_u^L(t)$$ is a constant, while the expression for the edge computation delay [TeX:] $$D e_u^L(t)$$ can be expanded as follows

The first term on the right-hand side of the above equation is the reciprocal of the matching variable x(t), and the second term is a linear function of x(t). Both terms are convex functions of x(t). Therefore, [TeX:] $$D e_u^E(t)$$ is a convex function in the subproblem SP1. Furthermore, according to [24], the smoothed function [TeX:] $$f_u^{D e}(\mathrm{x}, \mu, t)$$ is also a convex function of x(t).

Substituting the smooth function defined by (23) into [TeX:] $$\mathbf{SP1}$$ yields the following minimization problem.

In [TeX:] $$\mathbf{SP1.1}$$ [TeX:] $$f_u^{D e}(\mathrm{x}, \mu, t)$$ is a convex function with respect to [TeX:] $$x(t), \text { and } R_u(t)$$ is a linear function with respect to x(t). Therefore, the optimization objective (25a) is a convex function, and other related constraints are also either convex or linear functions. Hence, [TeX:] $$\mathbf{SP1.1}$$ is a standard convex optimization problem that can be solved using CVX.

2) Optimization of transmit power for each user.

Next, we integrate and simplify the problem [TeX:] $$\mathbf{P3}$$ by fixing [TeX:] $$x_{u, g}, D a_u \text { and } \eta_u,$$ then obtain the following sub-optimization problem

Similarly, since the expression of [TeX:] $$D e_u(t)$$ for a single slot is in the form of max[x, y], the subproblem [TeX:] $$\mathbf{SP2}$$ is also a min-max problem. Therefore, the log-sum-exp smoothing function is first used to transform [TeX:] $$D e_u(t)$$ as follows

From the expressions of [TeX:] $$D e_u^L(t) \text{ and } D e_u^E(t)$$ given in (6) and (7), it can be seen that the local delay [TeX:] $$D e_u^L(t)$$ in a single time slot is a constant, and there exists a fractional expression about the user’s transmit power in the edge computation delay [TeX:] $$D e_u^E(t)$$. Therefore [TeX:] $$D e_u^E(t)$$ is not a convex function with respect to the transmit power vector p(t), so the smoothing function (27) is also non-convex.

The following minimization problem is obtained by substituting the smoothing function defined in (27) into [TeX:] $$\mathbf{SP2}$$.

Since both [TeX:] $$f_u^{D e}(\mathbf{p}, \mu, t) \text { and } R_u(t)$$ in [TeX:] $$\mathbf{SP2.1}$$ are nonconvex functions with respect to the user’s transmit power p(t), [TeX:] $$\mathbf{SP2.1}$$ is a non-convex optimization problem. In this regard, we use the SCA to simplify the problem by constructing a more tractable convex approximation to the nonconvex objective function at the iteration point [TeX:] $$x^i.$$ For both the nonconvex terms [TeX:] $$f_u^{D e}(\mathbf{p}, \mu, t) \text { and }-R_u(t)$$ in [TeX:] $$\mathbf{SP2.1}$$, the second-order Taylor expansion at the [TeX:] $$i^{th}$$ iteration point [TeX:] $$\mathbf{p}^i$$ is performed as follows

Substituting the approximation functions given by (29) and (30) into the problem [TeX:] $$\mathbf{SP2.1}$$ yields the following approximation problem.

As the optimization objective in [TeX:] $$\mathbf{SP2.2}$$ is the linear sum of a series of convex functions and the constraints are also either convex orlinear, [TeX:] $$\mathbf{SP2.2}$$ is a convex optimization problem that can be solved using convex optimization tools such as CVX. By solving the convex approximation problem [TeX:] $$\mathbf{SP2.2}$$ in each iteration of the algorithm, the optimal solution of the original problem [TeX:] $$\mathbf{SP2}$$ can be obtained.

3) Optimization of amount of data processing for each user.

In order to avoid excessive backlog, ensure the stability of the data cache queuing system, and at the same time alleviate the increase in transmit power cost and the burden on the system and users due to the demand of data transmission rate, it is necessary to rationally arrange the amount of user data leaving [TeX:] $$Da(t)$$ per time slot. The following sub-optimization problem is obtained by fixing [TeX:] $$x_{u, g}, p_u \text { and } \eta_u$$ as follows

Similar to sub-problems [TeX:] $$\mathbf{SP1} \text{ and } \mathbf{SP2}$$, we first transform [TeX:] $$De_u(t)$$ using the log-sum-exp smoothing function.

By substituting (33) into [TeX:] $$\mathbf{SP3}$$, the following optimization problem can be obtained

From expressions (6) and (7), it can be seen that [TeX:] $$D e_u^L(t) \text{ and } D e_u^E(t)$$ are linear functions of the data processing amount Da(t). Therefore, [TeX:] $$f_u^{D e}(\mathbf{D a}, \mu, t)$$ is a convex function of Da(t). Based on this, (34a) is a combination of convex functions and linear functions, and thus [TeX:] $$\mathbf{SP3.1}$$ is a typical convex optimization problem, which can be solved using CVX tool.

4) Optimization of task splitting ratio for each user.

In the MEC service network constructed in this paper, ground users adopt partial offloading to perform task computation, so the following sub-optimization problem is obtained by further integrating and simplifying the problem [TeX:] $$\mathbf{P3}$$ with fixed [TeX:] $$x_{u, g}, p_u \text { and } D a_u$$ for each user.

Observing the problem formulation, it is found that [TeX:] $$\mathbf{SP4}$$ is different from the previous three sub-problems, as its optimization objective does not contain any content related to queue backlog. This is because queue backlog is only related to the arrival and departure of data in the queues, and the data offloading ratio will not affect the data arrival and departure in both the real and virtual queues. Therefore, there are no queue-related terms in [TeX:] $$\mathbf{SP4}$$.

Similar to the above subproblems, we first transform [TeX:] $$D e_u(t)$$ using the log-sum-exp smoothing function as follows

The next step is to substitute (36) into [TeX:] $$\mathbf{SP4}$$ to obtain the following optimization problem.

As shown in (6) and (7), both [TeX:] $$D e_u^L(t) \text { and } D e_u^E(t)$$ are linear functions with respect to the task splitting ratio [TeX:] $$\eta(t)$$. Therefore, [TeX:] $$f_u^{D e}(\boldsymbol{\eta}, \mu, t)$$ is a convex function with respect to [TeX:] $$\eta(t)$$. In addition, since all constraints are either convex or linear functions, sub-problem [TeX:] $$\mathbf{SP4.1}$$ is a convex optimization problem that can be solved using the CVX.

Based on the analysis of the three sub-problems [TeX:] $$\mathbf{SP1, SP2, SP3,} \text{ and } \mathbf{SP4}$$ above, the algorithm for minimizing the long-term time average QoS utility function of ground users based on Lyapunov optimization theory can be summarized as Algorithm 1.

Note that the parameter μ of the log-sum-exp smoothing function, i.e., [TeX:] $$\mu^i=\beta \mu^{i-1},$$ needs to be dynamically reduced during the iteration of individual time slots, as a way to make the solution results of individual subproblems as close to the optimal solution as possible.

D. Performance Analysis of Lyapunov Algorithm

The following theorem exists for the optimal QoS utility function value obtained from Algorithm 1:

Theorem 2. Assuming that in the SAGIN edge computing network established in this paper, the task data arrivals of each user and the air-to-ground channel state information are independently and identically distributed among each time slot, then the relationship between the long-term timeaveraged QoS utility function value obtained by Algorithm 1 and the actual optimal utility value [TeX:] $$I n_{\min }$$ of the original problem [TeX:] $$\mathbf{P1}$$ is as follows.

Proof. For any time slot t, we first expand (18) to obtain

Taking a telescoping summation operation over all time slots for the above equation yields

Afterwards, dividing both sides of (40) by V T at the same time and rearranging gives

Since [TeX:] $$I n_{\min }$$ is the actual optimal utility value of the original problem, making T converge to infinity simultaneously for both sides of (41) yields

Algorithm 1 based on Lyapunov optimization theory minimizes the QoS utility function while also guaranteeing the stability of the cache queue [TeX:] $$\boldsymbol{Q}(t)$$ and the virtual queue [TeX:] $$\boldsymbol{S}(t)$$, while according to [26], if a queue is rate stable, its arrival rate a(t) and departure rate b(t) satisfy

Intuitively, only when the average arrival rate of the queue is less than the average departure rate, the queue backlog not increase infinitely in the long-term time course, so the second and third terms on the far right of (42) are less than or equal to 0. Based on this, inequality (42) can be reduced to

So far, Theorem 2 is proved.

Theorem 2 shows that the Lyapunov optimization algorithm has a “performance-backlog” equilibrium relationship [O(1/V),O(V)]. When the penalty factor V increases, the optimal value obtained by Lyapunov algorithm is closer to the theoretical optimal value, while when V decreases, the system takes more consideration of queue stability and the queue backlog decreases.

IV. NUMERICAL RESULTS

In this section, we perform extensive simulation for the MEC service scenario in the SAGIN to evaluate the optimization performance of the proposed Algorithm 1 in terms of ensuring system stability and user QoS performance.

A. Parameters Setting

Considering that the connection of air-to-ground is approximately the same in each region on the ground, we perform the algorithm simulation in a unit latitude and longitude grid. We set the simulation area within China, and its longitude range is [TeX:] $$116^{\circ} \mathrm{E} \sim 117^{\circ} \mathrm{E}$$ and dimension range is [TeX:] $$39^{\circ} \mathrm{N} \sim 40^{\circ} \mathrm{N}.$$ Meanwhile, the satellite constellation configuration of the SAGIN is given by the improved indicator-based evolutionary algorithm (IIBEA) [27], and the constellation configuration obtained after 200 iterations of the IIBEA is taken here as the basis for simulation, which has a satellite number of 81, a deployment altitude of 1497.8 km and an operation period of about 115.78 minutes. In addition, the distribution of HABs and ground users in the simulation area obeys the uniform Poisson point distribution. Other simulation parameters are shown in TABLE I.

B. Algorithm Performance Validation

In order to analyze the performance variation of the proposed algorithm with respect to the penalty factor V , we select the long-term time-averaging delay, user-cost, and their weighted sum as performance metrics for simulation, and make the values of the penalty factor V fall within the interval [TeX:] $$\left[0.7 \times 10^{10}, 1.6 \times 10^{10}\right] .$$ As shown in Fig. 2, the values of all three metrics decrease as V increases, with the variation in user-cost being less pronounced. This is due to the fact that the user-cost largely comprises the satellite usage cost, which is dependent on the number of satellites connected to the HABs. Furthermore, the number of connections between HABs and satellites within a single time slot is solely determined by the satellite constellation. Consequently, it can be inferred that the long-term time-averaging LEO satellite usage cost for each user remains relatively consistent.And the user’s usage cost is mainly related to the local calculation frequency and transmitting power, so the reduction of the use-cost with V is smaller. In conclusion, with the increase of V , the optimization objective of the original problem can be better optimized, leading to improved user QoS performance.

As shown in Fig. 3, the accumulation of both real data queue backlog and virtual queue backlog escalates in tandem with the penalty factor V . This phenomenon stems from the inherent “performance-backlog” trade-off characteristic of the Lyapunov optimization algorithm. As V increases, the system’s queues tend to enter an unstable state, characterized by a substantial backlog. Moreover, from a task computation perspective, the increase of V prompts ground users to curtail task computations per time slot while diminishing transmit power, aiming to optimize latency and use-cost performance. Consequently, this action results in a reduction of both data departure and transmission throughput per time slot at the user end, thereby amplifying the backlog of the data queue and virtual queue.

Fig. 4 compares the partial offloading scheme proposed in this paper with the full offloading scheme (i.e., all computation tasks are offloaded to the MEC side for computation) as well as the full local computation scheme. It can be observed that the average queue backlog brought about by either the full offloading computing scheme or the full local computing scheme is larger than that of the partial offloading scheme adopted in this paper. In terms of user QoS performance, the partially offloaded scheme also outperforms the other two types of computation schemes in terms of delay performance.

C. System Queueing Stability Analysis

Fig. 5 shows the average throughput and the average delay of each user for a penalty factor [TeX:] $$V=0.7 \times 10^{10} \text { and } V=1.6 \times 10^{10}$$ over a range of 500 time slots. It can be seen that the user throughput and user delay in each time slot satisfy the C1 and C3 constraints in the original optimization problem [TeX:] $$\mathbf{P1}$$ for both [TeX:] $$V=0.7 \times 10^{10} \text { and } V=1.6 \times 10^{10}$$, verifying that the resource management algorithm proposed in this paper can achieve optimal user QoS performance while strictly satisfying the constraints. It is also shown that the virtual queue [TeX:] $$\boldsymbol{S}(t)$$ defined by (13) is able to achieve the guarantee of the longterm time-averaging constraint C1 through queue stability.

Fig. 6 represents the variation of the average real data cache queue backlog and the average virtual queue backlog for each user in the system over 500 time slots when the penalty factor parameter is taken as [TeX:] $$V=1.1 \times 10^{10}.$$ It can be seen that both the virtual queue and the real queue in the system have a stable upper bound, and the queue backlog does not keep increasing with time, where the average backlog of the real queue in 500 time slots is 35106 bits, and the average backlog of the virtual queue in 500 time slots is 5601 bps. This proves that the optimization algorithm based on Lyapunov optimization theory proposed in this paper can guarantee the stability of the system queue.

D. Queue Backlog Change Analysis

In order to have a deeper understanding of the algorithm optimization process and to explore some potential transmission phenomena in the network, the following data within the 0∼100 time slots are intercepted to further analyze the queue variation in SAGIN edge service system. By observing Fig. 7, we can find that the real data queue backlog of users is mainly controlled by the data leaving amount Da(t). Since the queue backlog in the simulation is recorded at the end of the time slot, if the current time slot has a large queue backlog, then the next time slot algorithm will increase the data leaving amount to reduce the queue backlog so as to ensure the stability of the queue.

However, in Fig. 8, during the change of the average backlog of the virtual queue, if the queue backlog of the current time slot is larger, the average queue departure (i.e., the average transmission rate of each user) of the next time slot will be smaller instead. Unlike in the real data queue where there is no direct connection between the data leaving amount of each user, in the virtual queue, there is a mutual constraint relationship between the virtual queue leaving amount of each user due to the existence of inter-user interference. When the queue backlog is large for a few users, these users increase the transmit power in the next time slot as a way to increase the throughput of the next time slot to relieve the queue backlog, however, the high transmit power also brings large interference, which affects the throughput of other users. Therefore, in the virtual queue, a large average queue backlog in the current time slot tends to bring about an imbalance in the transmit power among the users in the next time slot, thus making the average data rate of the users in the next time slot lower. In other words, in the virtual queue, when the queue backlog is larger for certain users, the algorithm focuses on taking care of these users by allocating a larger transmission rate for them, while making sacrifices for other virtual queues with smaller queue backlogs. It can also be seen from Fig. 8 that for a single user, the change in the virtual queue backlog still satisfies the criterion that if the current time slot backlog is large, the next time slot departure is large, which also shows that the increase in transmission rate for a single user does not represent an increase in the average rate for each user. This also reflects the common "overall optimal-individual optimal" compromise relationship in communication network optimization.

V. CONCLUSION

In order to further expand the computing capability of the equipment and reduce the delay of ground users, this paper proposes a network architecture combining MEC and SAGIN, and on this basis, taking into account the high mobility of LEO satellites. We carry out the refined decision-making of offloading computational tasks and the dynamic allocation of wireless resources in order to achieve the optimization of the QoS including the delay as well as the usage cost over a long term time horizon based on the Lyapunov optimization theory. Simulation results show that the method proposed can effectively improve users’ QoS performance while ensuring the system stability, which is superior to the all-local-computing or all-offloading-computing schemes.

APPENDIX A

For any [TeX:] $$Q \geq 0, b \geq 0, a \geq 0$$

The upper bound expression for the drift plus penalty function is obtained by expanding the Lyapunov drift using this inequality.

1) For the real data cache queue of ground users, an expansion based on the inequality yields

The last of these inequalities can be derived from the finite nature of the user data cache and the finite nature of the queue backlog. Firstly, in this paper, we assume that ground users are equipped with data cache units to store the arriving computation tasks, but the data cache capacity on the user side is limited, so the data arrivals per time slot a(t) cannot exceed the cache length, and thus there exists an upper bound [TeX:] $$a_{\max }(t) .$$ Secondly, for the leaving data amount Da(t) per time slot, due to the finite local computing capacity and MEC-side computing capacity, coupled with the existence of an upper bound on the user-side data cache capacity, it can be deduced that there is also an upper bound [TeX:] $$D a_{\max }(t)$$ on the leaving data amount Da(t) per time slot.

2) For the rate-constrained virtual queue of ground users, an expansion based on the inequality yields

The last of these inequalities is also upper bounded by the user transmission rate [TeX:] $$R_u(t)$$ due to the limited transmit power of the user and the limited transmission bandwidth, which is bounded as follows

The upper bound of the drift-plus-penalty function shown in (19) is obtained by substituting the inequality results of (46) and (47) into the drift-plus-penalty function of (18).