Towards Symbiotic SAGIN Through Inter-operator Resource and Service Sharing: Joint Orchestration of User Association and Radio Resources
Shizhao He Jungang Ge Ying-Chang Liang Dusit Niyato















 
theoremTheorem
corollaryCorollary[section]
definitionDefinition
propProposition
lemmaLemma
remarkRemark

rcl



[ [3]
gobble











Towards Symbiotic SAGIN Through Inter-operator Resource and Service Sharing: Joint Orchestration of User Association and Radio Resources
    Shizhao He, Jungang Ge, Ying-Chang Liang, Fellow, IEEE, and Dusit Niyato, Fellow, IEEE





This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which  this version may no longer be accessible.
S. He, J. Ge, and Y.-C. Liang are with the National Key Laboratory of Wireless Communications, and the Center for Intelligent Networking and Communications (CINC), University of Electronic Science and Technology of China (UESTC), Chengdu 611731, China (e-mail: heshizhao@std.uestc.edu.cn; gejungang@std.uestc.edu.cn; liangyc@ieee.org).

Dusit Niyato is with the School of Computer Science and Engineering, Nanyang Technological University (NTU), Singapore (e-mail:dniyato@ntu.edu.sg)
    April 29, 2024
==========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





The space-air-ground integrated network (SAGIN) is a pivotal architecture to support ubiquitous connectivity in the upcoming 6G era.
Inter-operator resource and service sharing is a promising way to realize such a huge network, utilizing resources efficiently and reducing construction costs.
Given the rationality of operators, the configuration of resources and services in SAGIN should focus on both the overall system performance and individual benefits of operators.
Motivated by emerging symbiotic communication facilitating mutual benefits across different radio systems, we investigate the resource and service sharing in SAGIN from a symbiotic communication perspective in this paper.
In particular, we consider a SAGIN consisting of a ground network operator (GNO) and a satellite network operator (SNO). Specifically, we aim to maximize the weighted sum rate (WSR) of the whole SAGIN by jointly optimizing the user association, resource allocation, and beamforming. Besides, we introduce a sharing coefficient to characterize the revenue of operators. Operators may suffer revenue loss when only focusing on maximizing the WSR. 
In pursuit of mutual benefits, we propose a mutual benefit constraint (MBC) to ensure that each operator obtains revenue gains.
Then, we develop a centralized algorithm based on the successive convex approximation (SCA) method. Considering that the centralized algorithm is difficult to implement, we propose a distributed algorithm based on Lagrangian dual decomposition and the consensus alternating direction method of multipliers (ADMM). Finally, we provide extensive numerical simulations to demonstrate the effectiveness of the two proposed algorithms, and the distributed optimization algorithm can approach the performance of the centralized one. The results also reveal that the proposed MBCs can enable operators to achieve mutual benefits and realize a symbiotic resource and service sharing paradigm.









  Symbiotic communication, SAGIN, inter-operator resource and service sharing, resource optimization.
  












§ INTRODUCTION
 
Wireless communication has made outstanding achievements in the past decades and stepped into the 5G era.
Although 5G can realize much better performance than former communication systems, including higher peak data rate, and ultra reliable and low latency, it still cannot satisfy some requirements in 6G, such as the ubiquitous connectivity <cit.>.
In fact, there are still 2.7 billion people who cannot access the Internet, and about 6% of the rural area worldwide lacks mobile network coverage <cit.>. The principal reasons are the limited coverage of ground networks and the high expenses of deploying networks in rural areas.
Satellite networks have emerged as an important complementary technology to enhance coverage of ground networks <cit.>, which can tackle the above issues. Compared with ground networks, satellite networks can provide seamless coverage to users without deploying costly fiber optic backhaul. 

Therefore, the space-air-ground integrated network (SAGIN) is a vital technology to realize ubiquitous connectivity in the 6G era.

Recently, satellite networks have rapidly developed and attracted considerable attention from academia to industry. For instance, SpaceX has launched over 4,000 low earth orbit (LEO) satellites to construct Starlink <cit.>. 
To realize ubiquitous connectivity, addressing the severe path loss caused by the long propagation distance is necessary.
In Starlink, a satellite terminal (ST) comprising a satellite dish and an access point is exploited to deal with this problem. The satellite dish provides high antenna gain to compensate for the path loss, and the access point offers access service to users. In <cit.>, a satellite access network for 5G and beyond is proposed based on this kind of ST.
Satellite networks are also expected to realize high throughput satellite communication. To this end, high frequency is adopted in satellite communication, such as Ka band <cit.>. Besides, multi-beam satellites also play an essential role in realizing high throughput satellite communication. Compared with conventional satellites, multi-beam satellites can generate multiple beams with high gain and thus provide better satellite communication <cit.>. 
Moreover, the shadowed-Rician (SR) fading is widely adopted to model the satellite-ground channels, based on which many works analyze the performance of satellite networks <cit.>.
Resource allocation in satellite networks is also widely studied <cit.>. Additionally, diverse multiple access techniques are adopted in satellite networks to serve users <cit.>, such as the time division multiple access (TDMA) manner.













The rapid growth of satellite networks paves the way for realizing SAGIN, which has been widely studied in recent years. Since SAGIN is a large-scale network, it is almost impossible for an individual network operator to construct it. 

Therefore, inter-operator sharing[For ease of notation, “inter-operator sharing” is synonymous with “inter-operator resource and service sharing” in this paper.] is a promising technology for constructing SAGIN.

In fact, the famous ground network operator (GNO) T-mobile announces that it aims to provide seamless coverage to users by cooperating with satellite network operator (SNO) SpaceX <cit.>. 

Optimizing the configuration of resources and services is an essential problem to fully exploit the advantages of SAGIN with inter-operator sharing.

Considering the limited spectrum, many works study the spectrum sharing problem in SAGIN <cit.>.

In <cit.>, catering to the heterogeneity and high dynamics of SAGIN, the authors propose an intelligent spectrum management framework empowered by artificial intelligence and software defined network (SDN). 
In <cit.>, the spectrum allocation is optimized to maximize the sum rate while reducing the interference received at satellite users.


On the other hand, service sharing is another critical topic in studying SAGIN. Many recent scientific literatures demonstrate that satellite networks are a promising supplement to ground networks in providing diversity services to remote areas, such as backhaul connectivity <cit.> and task offloading <cit.>. In <cit.>, a satellite assists ground networks in providing backhaul links to users in remote areas, and both user association and beamforming are optimized to maximize the sum rate.
In <cit.>, satellite networks support task offloading of remote Internet-of-Things (IoT) devices. An optimization problem about user association and spectrum allocation is formulated to minimize the  task offloading delay. The above works only focus on the overall performance of SAGIN. However, motivating different operators to share resources and services is also a vital problem. As the operators are inherently competitors with conflicting interests, they will be reluctant to construct SAGIN if resource and service sharing cannot guarantee their revenue, i.e., they cannot achieve mutual benefits.
In <cit.>, a pricing mechanism is proposed to ensure each 
operator's revenue. However, it cannot fully leverage the resources in SAGIN to achieve a high overall performance.
Therefore, it is vital to simultaneously achieve mutual benefits among operators and maximize the overall performance of SAGIN. 















The recently emerging symbiotic communication provides a promising paradigm to tackle the above problem, which aims to optimize the collective objectives of diverse radio systems and realize mutual benefits among them via sharing resources and services <cit.>.

Symbiotic communication has been widely studied in backscatter communications <cit.>. In <cit.>, the authors investigate the symbiotic communication in a passive IoT system, consisting of an active transmission and a backscatter transmission. Specifically, the active transmission shares spectrum and energy resources with the backscatter transmission, and in return the backscatter transmission provides a beneficial multi-path to the active transmission via reflecting the incident signal. As a result, the active and backscatter transmission achieve mutual benefits in this system, realizing higher sum rate.

Besides, the symbiotic communication is realized in the coexistence of a cellular network and WiFi system in <cit.>. Particularly, the WiFi system shares spectrum with the cellular network, and the cellular network adjusts transmission parameters to guarantee the performance of WiFi system.

Resource allocation in symbiotic communication is also studied.
In <cit.>, both user association and beamforming are optimized to realize the symbiotic communication in multi-operator cellular networks. Symbiotic communication can also be a prominent solution to construct SAGIN.

In <cit.>, the authors propose a blockchain and deep learning based framework to optimize resource allocation in SAGIN with the guidance of symbiotic communication.
However, it lacks an explicit formulation to strictly guarantee individual benefits after sharing resources and services.


Motivated by the above works, in this paper, we investigate joint user association, resource allocation, and beamforming design in SAGIN from a symbiotic communication perspective. Specifically, the considered SAGIN system consists of a GNO and SNO engaging in inter-operator sharing, i.e., sharing spectrum and allowing users to access arbitrary networks. The GNO has base stations (BSs) connected to the core network by fiber. On the other hand, the SNO has STs connected to the core network via a multi-beam LEO satellite, and this LEO satellite serves STs in the TDMA manner. Given the high cost of deploying fiber in remote areas, the number of BSs is fewer than that of STs.
Our objective is to maximize the weighted sum rate (WSR) of all users while achieving mutual benefits among operators in terms of revenue. To characterize revenue of operators, we introduce a sharing coefficient that reflects the impact of inter-operator sharing.

The revenue may experience losses when only pursuing the maximization of the WSR.

To achieve mutual benefits, we propose a mutual benefit constraint (MBC) to guarantee that each operator can obtain higher revenue with inter-operator sharing.





The main contributions of this work can be summarized as follows:

  
  * We investigate a symbiotic resource and service sharing paradigm in the SAGIN, which aims to maximize the overall system performance and realize mutual benefits between the GNO and SNO. We formulate an optimization problem regarding user association, resource allocation, and beamforming design to maximize the WSR, subject to the MBCs and STs' backhaul capacity constraints.
  
  * We introduce a sharing coefficient to characterize the revenue of each network operator, based on which the MBC is formulated to ensure each operator can obtain revenue gains after sharing spectrum and services.
  
  * We propose a centralized algorithm based on successive convex approximation (SCA) for the formulated WSR maximization problem. Concerning the large scale of the SAGIN, it is challenging to implement the centralized algorithm. Then, we develop a distributed algorithm, which consists of the Lagrangian dual decomposition based algorithm for user association and consensus alternating direction method of multipliers (ADMM) based algorithm for beamforming design and resource allocation.
  
  * Simulation results show that both the proposed algorithms can outperform the considered benchmarks, and the distributed algorithm approaches a performance close to that of the centralized algorithm. Moreover, in comparison to operators in the non-symbiotic case, each operator in the symbiotic case attains higher revenue after sharing spectrum and services. This implies that operators achieve mutual benefits and realize a symbiotic resource and service sharing paradigm in the SAGIN. 

  


The remainder of this paper is organized as follows. Section <ref> illustrates the system model. In Section <ref>, we formulate the WSR maximization problem for the considered SAGIN system. In Section <ref>, an SCA-based centralized algorithm and finding initial point algorithm are proposed for the WSR maximization problem. In Section <ref>, a distributed algorithm based on consensus ADMM and Lagrangian dual decomposition methods is developed. In Section <ref>, the extensive simulation results are presented.

Finally, Section <ref> concludes this paper.

Notations: The notations in this paper are listed as follows. The scalars, column vectors and matrices are represented by lowercase, bold lowercase and uppercase symbols (e.g., x, 𝐱 and 𝐗), respectively. |𝒜| denotes the cardinality of set 𝒜. ℂ denotes the set of complex numbers. Main notations are summarized in Table <ref>.







§ SYSTEM MODEL






 §.§ SAGIN Network

As shown in Fig. <ref>, we consider the downlink of a SAGIN system, which consists of a GNO and SNO, and we define 𝒪={G,S} as the set of operators. Specifically, the GNO has N_U users and deploys N_B BSs, directly connected to the core network via fiber that can provide large backhaul capacity. On the other hand, the SNO has N_U^' subscribed users and deploys N_B^' STs to serve users in remote areas, and it exploits a multi-beam LEO satellite to provide backhaul links to STs over the Ka band. Since the payload of the satellite is constrained, the backhaul capacity it provides to STs is limited. Besides, both operators have a dedicated C band for serving users.
We denote ℬ_G={1,2,…,N_B} and ℬ_S={N_B+1,N_B+2,…,N_B+N_B^'} as the set of all BSs and STs, respectively. Denote the set of subscribed users of the GNO as 𝒰_G={1,2,…,N_U} and that of the SNO as 𝒰_S={N_U+1,N_U+2,…,N_U+N_U^'}. Then, we denote ℬ≜∪_n∈𝒪ℬ_n and 𝒰≜∪_n∈𝒪𝒰_n. We consider that all BSs and STs are equipped with N_t transmit antennas while all users are equipped with a single antenna.
In the considered SAGIN, inter-operator sharing is adopted, where all BSs and STs share their C bands to serve users, and the users can access arbitrary BSs or STs for service. Furthermore, we assume that each operator has a control center to
collect necessary information (e.g., channel state information (CSI), and user locations.).











 §.§ Space-Ground Communication Model

We consider that the multi-beam LEO satellite is equipped with N_L feeds, and each can generate a spot beam to provide backhaul links to STs. Specifically, each beam serves STs within its coverage in the TDMA manner. Besides, universal frequency reuse is adopted in space-ground communication, which means all beams operate over the Ka band. Then, we assume each ST is pointed towards the satellite, and the received signal of ST i within the coverage of beam ℓ is

    y_ℓ, i   = √(p_ℓG_T(ζ_ℓ,i)G_R) f_i s_ℓ
        + ∑_ℓ^'=1,
    ℓ^'≠ℓ^N_L
        √(p_ℓ^'G_T(ζ_ℓ^',i)G_R) f_i s_ℓ^' + z_i,

where p_ℓ is the transmit power of ℓ-th spot beam; 

s_ℓ is the symbol transmitted by ℓ-th spot beam, which is assumed to be an independent random variable with zero mean and unit variance; and z_i∼𝒞𝒩(0, σ_s^2) is the additive white Gaussian noise (AWGN).
G_R is the receive antenna gain of STs, G_T(ζ_ℓ,i) is the transmit antenna gain of the ℓ-th spot beam toward ST i, and ζ_ℓ,i is the off-boresight angle between the center of ℓ-th spot beam and ST i. Specifically, G_T(ζ_ℓ,i) is given by

    G_T(ζ_ℓ,i) = G_T^0 · 4 |J_1(κ asinζ_ℓ,i)/κ asinζ_ℓ,i|^2,

where G_T^0 is the maximum transmit gain, κ = 2π/λ_Ka is the wave number, a is the radius of the dish antenna, and λ_Ka is the carrier length of Ka band.
f_i is the downlink channel between the LEO satellite and ST i, which is given by

    f_i = √(d_i)f̃_i,

where d_i accounts for the path-loss between the satellite and ST i, while f̃_i is the SR fading. The probability distribution function of SR fading gain |f̃_i |^2 is given by

    f_|f̃_i |^2(s;b,m,Ω)    = 1/2b(2bm/2bm+Ω)^mexp(-s/2b)×
        _1F_1(m,1,Ω s/2b(2bm+Ω)),

where Ω is the average power of a LoS component, 2b is the average power of the scatter components, m is the Nakagami parameter, and _1F_1(·,·,·) is the confluent hypergeometric function. Based on (<ref>), the backhaul capacity for ST i served by ℓ-th spot beam is

    C_ℓ,i   = B_Kat_ℓ,i
       ×log_2(1+p_ℓ G_T(ζ_ℓ,i)G_R| f_i |^2/∑_ℓ^'≠ℓp_ℓ^'G_T(θ_ℓ^',i)G_R| f_i |^2 + σ_s^2),

where B_Ka denotes the bandwidth of Ka band, t_ℓ,i is the fraction of time resource ℓ-th spot beam allocates to ST i.



 §.§ Ground-Ground Communication Model


Universal frequency reuse is also adopted by BSs/STs to serve users, and the downlink channel between BS/ST i and user k over the spectrum band of operator n∈𝒪 is given by

    𝐡_i,k^n = √(d_i,k)𝐡̃_i,k^n,

where d_i,k is the path-loss between BS/ST i and user k, and 𝐡̃_i,k^n∈ℂ^1× N_t is the small-scale fading coefficient. 
With inter-operator sharing, the users of one operator can access other operators' networks, and each BS/ST can operate over all spectrum bands. In this case, the received signal of user k associated with BS/ST i∈ℬ_m over the band n is 

    y_i,k^n    = 𝐡_i,k^n 𝐰_i,k^ns_k
       + ∑_j∈ℬ_m,k^'∈𝒰
    (j,k^')≠(i,k)
        𝐡_j,k^n 𝐰_j,k^'^n s_k^'
        + ∑_j∈ℬ∖ℬ_m,k^'∈𝒰
        𝐡_j,k^n 𝐰_j,k^'^n s_k^' + z^n,

where 𝐰_i,k^n∈ℂ^N_t × 1 is the transmit beamforming that BS/ST i assigns to user k over the band n; s_k denotes the symbol transmitted to user k; and z^n ∼𝒞𝒩(0, σ_t^2) is the AWGN over band n.
The instantaneous signal-to-interference-plus-noise ratio (SINR) of user k in the band n can be written as

    γ_i,k^n = |𝐡_i,k^n 𝐰_i,k^n | ^2/Ĩ_i,k^n + I_i,k^n + σ_t^2,

where Ĩ_i,k^n denotes the intra-operator interference while I_i,k^n denotes the inter-operator interference. Specifically, Ĩ_i,k^n and I_i,k^n are as follows:

  
    Ĩ_i,k^n = ∑_j∈ℬ_m,k^'∈𝒰
    (j,k^')≠(i,k)|𝐡_j,k^n 𝐰_j,k^'^n| ^2, 
       I_i,k^n = ∑_j∈ℬ∖ℬ_m,k^'∈𝒰|𝐡_j,k^n𝐰_j,k^'^n| ^2.


Thus, the achievable rate for user k with inter-operator sharing is shown as

    R_i,k = ∑_n∈𝒪log_2(1 + γ_i,k^n).









§ PROBLEM FORMULATION


As the network performance heavily depends on the configuration of resources and services, we formulate a WSR maximization (WSRM) problem through jointly optimizing user association, resource allocation, and beamforming design. Specifically, since the backhaul capacity offered by the LEO satellite is limited, the investigated WSRM problem is subject to the backhaul capacity constraints of STs. Besides, we introduce a sharing coefficient δ_z for operator z to characterize its revenue. Based on the revenue, MBCs are proposed to realize mutual benefits, realizing a symbiotic resource and service sharing paradigm.





 §.§ Revenue of Operators

While inter-operator sharing can improve the overall performance, it concurrently incurs additional costs due to inter-operator interference and serving other operators' users.
Hence, to compensate for the additional costs and incentivize operators to engage in inter-operator sharing,
a sharing coefficient δ_z ∈ [0, 1] is introduced for each operator z to characterize its revenue.
As the revenue of an operator is usually related to the transmission rate that it can offer, we define the revenue of operator z as

    U_z =∑_i∈ℬ,k∈𝒰x_i,kα_i,k^zR_i,k,

where 

x_i,k∈{0,1} is the user association variable, x_i,k=1 if user k is associated with BS/ST i, and x_i,k=0 otherwise. Let BS/ST i belong to operator m, and user k subscribe to operator n, then α_i,k^z is defined as

    α_i,k^z = 
                          1,     z=m=n,
    δ_m,     z=m,m ≠ n,
    
                          1-δ_m,     z=n,m ≠ n,
    
                          0,

From (<ref>) and (<ref>), it can be observed that the sharing coefficient δ_m can affect 
the revenue (or compensation) obtained by operator m from serving other operators' users. Thus, this sharing coefficient has a significant impact on service sharing. 
For example, when δ_m=0, other operators need not compensate operator m for offloading their users to it, and operator m 
cannot obtain any revenue. 
As a result, for higher revenue, operator m prefers to serve its users, while other operators prefer to offload their users to operator m.



 §.§ WSRM Formulation

In the studied SAGIN system, each user can be associated with any BSs/STs, and each BS/ST can operate in all bands. To efficiently maximize performance of the whole system while ensuring that operators achieve mutual benefits, we formulate the WSRM problem from the symbiotic communication perspective, which is written as

  
    𝐏1:   max_𝐩,𝐭,𝐖,𝐱∑_i∈ℬ∑_k∈𝒰 x_i,kb_i,kR_i,k
       s.t.  U_z ≥ U_z^0, ∀ z ∈𝒪, 
          ∑_i ∈ℬ_S t_ℓ,i≤ 1, ∀ℓ=1,⋯,N_L, 
          ∑_k∈𝒰 B_CR_i,k≤ C_ℓ,i, ∀ i ∈ℬ_S, 
          ∑_n∈𝒪∑_k∈𝒰∥𝐰_i,k^n∥^2 ≤ P_i, ∀ i ∈ℬ, 
          ∑_ℓ=1^N_L p_ℓ≤ P_Sat, 
          ∑_i∈ℬ x_i,k≤ 1, ∀ k ∈𝒰, 
              x_i,k∈{0, 1}, ∀ i∈ℬ,∀ k ∈𝒰,


where 𝐩≜{p_ℓ}_ℓ=1,2,⋯,N_L, 𝐭≜{t_ℓ,i}_ℓ=1,2,⋯,N_L, i∈ℬ_S, 𝐰≜{𝐰_i,k^n}_n∈𝒪, i∈ℬ,k∈𝒰, 𝐱≜{x_i,k}_∈ℬ,k∈𝒰, b_i,k denotes the weight for the rate of user k associated with BS i, and P_i and P_Sat are the maximum transmit power of BS/ST i and the LEO satellite, respectively. Besides, (<ref>) is the time resource allocation constraint; (<ref>) is the backhaul constraint for each ST, and B_C is the bandwidth of C band; (<ref>) is the power constraint for each BS/ST, which indicates the spectrum sharing is enabled; (<ref>) is the power constraint for the LEO satellite; (<ref>) means that each user can be associated with only one BS/ST, which indicates the service sharing is enabled, and (<ref>) indicates that x_i,k is a binary variable. Further, (<ref>) is referred to as the MBC, which ensures that each operator can obtain benefits. Specifically, U_z is the revenue operator z can obtain with inter-operator sharing, while U_z^0 is the revenue obtained without inter-operator sharing. The details of determining U_z^0 will be illustrated in Section <ref>.


In the above formulation, we do not take x_i,k into consideration in the left-hand-side (LHS) of constraint (<ref>) because R_i,k will automatically be 0 when x_i,k=0 at the optimality. To see this, consider the optimal beamforming and user association scheme of P1 for user k is x_i,k=0 while 𝐰_i,k^n≠0. Since x_i,k is considered in the objective (<ref>), by setting 𝐰_i,k^n=0, a feasible solution with higher objective value can be obtained. Then, according to (<ref>) and (<ref>), R_i,k will be 0, too.



§ CENTRALIZED OPTIMIZATION ALGORITHMS FOR WSRM


In this section, we leverage the block successive maximization (BSM) method <cit.> to solve P1, where the above problems are solved in three stages iteratively. 
In the first stage, the beamforming 𝐰 of BSs/STs and power 𝐩 of the multi-beam LEO satellite are optimized for the given user association scheme 𝐱 and time allocation 𝐭 using the SCA method. 
Then, 𝐰 and 𝐭 are optimized with fixed 𝐱 and 𝐩.
In the last stage, 𝐱 is optimized for the given beamforming design and resource allocation. 
Next, we also provide a method to obtain the revenue for each operator without inter-operator sharing (U_z^0).
Besides, we propose an initial point searching method for the SCA-based algorithm.


Finally, the proposed centralized algorithm can be implemented in the assigned control center.








 §.§ Resource Allocation and Beamforming Optimization for WSRM




P1 is a non-convex and mixed-integer problem, which is hard to tackle. To make the problem tractable, 

we introduce some auxiliary variables and rewrite P1 equivalently as

  
    𝐏1.1:   max_𝐩,𝐭,𝐰,𝐱,Γ,β,φ,ρ∑_i∈ℬ∑_k∈𝒰 x_i,kb_i,kr_i,k
       s.t. (<ref>), (<ref>), (<ref>), (<ref>), 
          Γ_i,k^n ≤|𝐡_i,k^n 𝐰_i,k^n|^2/β_i,k^n, 
          β_i,k^n ≥Ĩ_i,k^n + I_i,k^n + σ_t^2, 
          φ_i,k^n ≥|𝐡_i,k^n 𝐰_i,k^n|^2/ρ_i,k^n, 
          ρ_i,k^n ≤Ĩ_i,k^n + I_i,k^n + σ_t^2, 
          ∑_i∈ℬ∑_k∈𝒰x_i,kα_i,k^zr_i,k≥ U_z^0, 
          ∑_n∈𝒪∑_k∈𝒰 B_Clog_2(1+φ_i,k^n) ≤ C_ℓ,i, 
    
        
                  r_i,k=∑_n∈𝒪log_2(1+Γ_i,k^n),


where Γ≜{Γ_i,k^n}_n∈𝒪, i∈ℬ,k∈𝒰 represents the lower bound of SINR for user k associated with BS/ST i over band n; β≜{β_i,k^n}_n∈𝒪, i∈ℬ,k∈𝒰 represents the upper bound of total interference and noise received at user k; φ≜{φ_i,k^n}_n∈𝒪, i∈ℬ,k∈𝒰 represents the upper bound of SINR for user k associated with BS/ST i over band n; ρ≜{ρ_i,k^n}_n∈𝒪, i∈ℬ,k∈𝒰 represents the lower bound of interference and noise received at user k.


By introducing Γ and φ, the MBC (<ref>) and backhaul constraint (<ref>) can be satisfied at the optimality, respectively. Besides, it is obvious that the constraints from (<ref>) to (<ref>)  should hold with equality at the optimality. Thus, the equivalence between (<ref>) and (<ref>) is guaranteed.
P1.1 is still challenging to tackle due to the non-convex constraints (<ref>), (<ref>) and (<ref>). To make P1.1 more tractable, we propose to exploit the SCA method to deal with these non-convex constraints.

Note that the right-hand-side (RHS) of (<ref>) is a joint convex quadratic-over-linear function with respect to 𝐰_i,k^n and β_i,k^n. With this fact, we can use the SCA method to approximate the RHS of (<ref>) with a lower bound. The RHS of (<ref>) can be approximated with the first-order Taylor expression as

    |𝐡_i,k^n 𝐰_i,k^n|^2/β_i,k^n≥ g_1(𝐰_i,k^n, β_i,k^n),

where g_1(𝐰_i,k^n, β_i,k^n) is given in (<ref>) on top of the next page,

𝐰_i,k^n[τ-1], β_i,k^n[τ-1] are the iterative optimization variables obtained in iteration τ-1, and 𝐇_i,k^n≜(𝐡_i,k^n)^H𝐡_i,k^n. Note that the RHS of (<ref>) is also a convex function with respect to 𝐰, and thus its lower bound can be approximated as

    Ĩ_i,k^n + I_i,k^n + σ_t^2 ≥∑_(j,k^')≠(i,k)g_2(𝐰_j,k^'^n) + σ_t^2,

where g_2(𝐰_j,k^'^n) is given in (<ref>) on top of the next page.
The RHS of (<ref>) is a concave function with respect to φ_i,k^n and its upper bound can be obtained as

    log_2(1+φ_i,k^n)≤ g_3(φ_i,k^n)
       ≜log_2(1+φ_i,k^n[τ-1])+ φ_i,k^n-φ_i,k^n[τ-1]/ln2(1+φ_i,k^n[τ-1]).

The LHS of (<ref>) can be expressed as 
the difference of concave functions, and its lower bound is as follows:

    C_ℓ,i   ≥ C_ℓ,i(𝐩) 
       ≜ B_Kat_ℓ,i{log_2(∑_ℓ^'p_ℓ^'G_T(θ_ℓ^',i)G_R| f_i|^2+σ_s^2) 
        - C_ℓ,i^1 -  C_ℓ,i^2},

where 
C_ℓ,i^1 and C_ℓ,i^2 are given by

    C_ℓ,i^1     = log_2(∑_ℓ^'≠ℓp_ℓ^'[τ-1]G_T(θ_ℓ^',i)G_R| f_i|^2+σ_s^2), 
    
        C_ℓ,i^2     = ∑_ℓ^'≠ℓG_T(θ_ℓ^',i)G_R| f_i|^2/ln2(∑_ℓ^'≠ℓp_ℓ^'[τ-1]G_T(θ_ℓ^',i)G_R| f_i|^2+σ_s^2)
       ×(p_ℓ^'-p_ℓ^'[τ-1]).











With approximations (<ref>), (<ref>), (<ref>) and (<ref>), we can obtain the reformulated optimization problem at iteration τ of SCA as follows:

  
    𝐏1.2:   max_𝐩,𝐭,𝐰̃∑_i∈ℬ∑_k∈𝒰 x_i,kb_i,kr_i,k
       s.t. (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), 
          Γ_i,k^n ≤ g_1(𝐰_i,k^n, β_i,k^n), 
          ρ_i,k^n ≤∑_(j,k^')≠(i,k)g_2(𝐰_j,k^'^n)+σ_t^2, 
          ∑_n∈𝒪∑_k∈𝒰 B_Cg_3(φ_i,k^n) ≤ C_ℓ,i(𝐩).


For notational simplicity, we denote variables Γ, β,φ, ρ, and 𝐰 as 𝐰̃={𝐰,Γ, β,φ,ρ}, because those variables all depend on 𝐰.
The problem (<ref>) is still not convex because variables 𝐩,𝐭 and φ are coupled due to the constraint (<ref>). It can be observed that (<ref>) is convex with respect to either {𝐩, 𝐰̃} or {𝐭, 𝐰̃}.
Based on this observation, the BSM method is adopted to solve (<ref>). Specifically, {𝐩,𝐭, 𝐰̃} can be updated in the following manner: given {𝐭,𝐱}, optimize {𝐩, 𝐰̃}; and given {𝐩,𝐱}, optimize {𝐭, 𝐰̃}. 









 §.§ User Association Optimization for WSRM


When variables {𝐩,𝐭, 𝐰̃} are fixed, (<ref>) is an integer programming problem with respect to 𝐱, which is challenging to tackle. Therefore, we first relax the binary variable constraint (<ref>) into 0≤ x_i,k≤ 1. Then, to obtain a high-quality solution for the original integer programming problem, a penalty term ϱ(x_i,k^2-x_i,k) is introduced into the objective (<ref>) to force x_i,k to be binary, where ϱ is a positive penalty parameter. The reformulated problem is written as follows:

  
    𝐏1.3:   max_𝐱∑_i∈ℬ∑_k∈𝒰 x_i,kb_i,kR_i,k + ϱ∑_i∈ℬ∑_k∈𝒰 (x_i,k^2-x_i,k) 
       s.t. (<ref>), (<ref>),(<ref>),
              0 ≤ x_i,k≤ 1.


Considering the penalty term is a convex function of x_i,k, we approximate it with the SCA method, which is given by 

    (x_i,k^2-x_i,k)    ≥ f(x_i,k) 
       ≜ (2x_i,k[τ-1]-1)x_i,k
        - x^2_i,k[τ-1],

where x_i,k[τ-1] is the iterative optimization variable obtained in iteration τ-1.
With (<ref>), we can write the problem at iteration τ of the SCA as

    𝐏1.4:   max_𝐱∑_i∈ℬ∑_k∈𝒰 x_i,kb_i,kR_i,k + ϱ∑_i∈ℬ∑_k∈𝒰 f(x_i,k) 
       s.t. (<ref>), (<ref>).

P1.4 is a linear programming problem with respect to 𝐱 and can be easily solved. 





Besides, it can be verified that, when ϱ→∞, the solution 𝐱 of P1.4 will satisfy the binary constraint (<ref>) <cit.>. However, if ϱ is too large, the objective of (<ref>) will be dominated by the penalty term, and the original objective (i.e., WSR) will be diminished. To avoid this, we first initialize ϱ to a small value to find a good starting point and then gradually increase ϱ as: ϱ = qϱ, where q>1.
The overall algorithm is summarized in Algorithm <ref>.






 §.§ Obtaining Revenue without Inter-operator Sharing


U_z^0 can be determined by solving the WSRM problem without inter-operator sharing, which can be formulated as follows

  
    𝐏2:   max_𝐩,𝐭,𝐖,𝐱∑_i∈ℬ_z∑_k∈𝒰_z x_i,kb_i,kR_i,k
       s.t. (<ref>), (<ref>), (<ref>), (<ref>),  
    
              ∑_k∈𝒰_z∥𝐰_i,k^z∥^2 ≤ P_i, ∀ i ∈ℬ_z, 
              x_i,k = 
                   0  or  1,     i∈ℬ_z, k∈𝒰_z, 
    
                   0,    ,


where (<ref>) indicates there is no spectrum sharing, and (<ref>) indicates that each user can only be associated with the BSs/STs of the same operator. Further, R_i,k in (<ref>) is given by

    R_i,k = log_2(1 + γ_i,k^n),

where γ_i,k^n can be obtained by (<ref>) with I_i,k^n=0.

P2 can be decomposed into independent subproblems for different operators. 
Similar to P1, the subproblem for operator n can be equivalently reformulated as

    𝐏2.1:   max_𝐩,𝐭,𝐰̃,𝐱∑_i∈ℬ_z∑_k∈𝒰_z x_i,kb_i,kr_i,k
       s.t. (<ref>), (<ref>), (<ref>), (<ref>),(<ref>),
          (<ref>)-(<ref>).

It is clear that P2.1 can be regarded as a special case of P1.2.
Therefore,the algorithm proposed for P1.2 can be applied to solving P2.1. 
Different from P1.2, MBCs are not considered in P2.1, and there exists optimal 𝐱 of P2.1 are binary variables. Hence, we do not need to introduce the penalty term (<ref>) to the objective of (<ref>). 

If we solve 𝐱 directly to the global optimality, 𝐱 will not change during the algorithm and fix at 0 or 1. To find a good 𝐱, we exploit the gradient projection <cit.> for updating 𝐱, which is given by



    x_i,k = P_Ω_x(x_i,k+η b_i,kR_i,k),

where 
P_Ω_x(·) is the projection to the set Ω_x={x_i,k|(<ref>),(<ref>)}, and η is the step size. 
The overall algorithm is summarized in Algorithm <ref>.





 §.§ Finding Initial Points

The feasible initial points are important for the SCA-based algorithm. However, due to the backhaul constraint (<ref>) and MBC (<ref>), it is challenging to find feasible initial points for P1.2. To address this issue, we propose a penalty method to find the feasible initial points <cit.>. 
Specifically, we can obtain initial points by solving the following optimization problem

  
    𝐏3:   max_𝐩,𝐭,𝐰̃,𝐱,𝐬
        ∑_i∈ℬ∑_k∈𝒰 x_i,kb_i,kr_i,k - ξ (∑_z∈𝒪 s_u_z + ∑_i∈ℬ_S s_b_i) 
       s.t. (<ref>), (<ref>), (<ref>), (<ref>), (<ref>)
          ∑_i∈ℬ∑_k∈𝒰 x_i,kα_i,k^z r_i,k + s_u_z≥ U_z^0, 
          ∑_n∈𝒪∑_k∈𝒰 B_Cg_3(φ_i,k^n) ≤ C_ℓ,i(𝐩) + s_b_i, 
              s_u_z, s_b_i≥ 0,


where ξ is a positive penalty parameter and 𝐬≜{s_u_z,s_b_i}_z∈𝒪,i∈ℬ_S are slack variables. The basic idea of the above formulation is forcing the penalty terms to be zero by maximizing (<ref>). P3 can also be solved with the BSM method. Besides, the gradient projection can be exploited for updating 𝐱 and 𝐬, which is given by



    (x_i,k,s_u_z) = P_Ω_xs(x_i,k+η b_i,kR_i,k,s_u_z-ηξ),

where P_Ω_xs(·) is the projection to the set Ω_xs={x_i,k,s_u_z|(<ref>), (<ref>), (<ref>)}. 
The overall algorithm is summarized in Algorithm <ref>. 








 §.§ Convergence and Complexity Analysis

Algorithms proposed in this section have similar main process. As a result, we mainly analyze the convergence and complexity of Algorithm <ref>, and the results can be applied to Algorithm <ref> and <ref>.
The convergence of Algorithm <ref> is provided as follows.

  Algorithm <ref> is guaranteed to converge


  Please refer to Appendix <ref>






The main complexity of Algorithm <ref> is caused by solving (<ref>) and (<ref>) in the inner loop. The worst-case complexity of Algorithm <ref> depends on the number of variables, which can be upper bound as

    𝒪{I_out^1I_in^1(   (N_L + N_t(N_B+N_B^')(N_U+N_U^'))^4 
       + (N_B^' + N_t(N_B+N_B^')(N_U+N_U^'))^4 
       + ((N_B+N_B^')(N_U+N_U^'))^4 )}

where I_out^1 and I_in^1 are the number of outer and inner iterations, respectively.




§ DISTRIBUTED ALGORITHMS FOR WSRM

The centralized algorithm proposed in Section <ref> is challenging to implement in the real SAGIN.

In this section, we propose a distributed algorithm for the WSRM problem, which can be implemented in local centers independently. Specifically, we develop a distributed user association algorithm based on the Lagrangian dual decomposition and a distributed algorithm for beamforming design and resource allocation based on the consensus ADMM.



 §.§ Lagrangian Dual Decomposition-based User Association Design

With fixed {𝐩, 𝐭, 𝐰̃}, P1.1 can be reformulated as 

  
    𝐏1.5:   max_𝐱∑_i∈ℬ∑_k∈𝒰 x_i,kb_i,kR_i,k
       s.t. (<ref>), (<ref>), (<ref>),


It is clear that P1.5 is non-convex due to the binary variables 𝐱. Although the centralized algorithm can tackle this problem efficiently by relaxing binary variables, it is difficult to implement in the real SAGIN due to high complexity. 


Keeping that in mind, we propose a distributed user association algorithm based on the Lagrangian dual decomposition <cit.>. 

First, the Lagrangian function with respect to constraint (<ref>) is 

    L(𝐱, λ) = ∑_i∈ℬ∑_k∈𝒰 x_i,kb_i,kR_i,k + ∑_z∈𝒪λ_z(U_z - U_z^0),

where λ are dual variables introduced for (<ref>). Then the dual function is give as

    g(λ) =    max_𝐱 L(𝐱, λ)
       s.t. (<ref>), (<ref>).

Based on the dual function, the Lagrange dual problem of P1.5 can be written as

    min_λ g(λ).

Now, the solution to P1.5 can be obtained by solving (<ref>). Note that since P1.5 is non-convex and discrete in nature, solving (<ref>) is not equivalent to solving P1.5, and there may exist a duality gap. Nevertheless, good primal solutions can often be obtained by solving the dual problem <cit.>.

Compared with P1.3, (<ref>) can be decomposed into independent subproblems for each user and BS/ST and then solved in a distributed manner. Specifically, (<ref>) can be solved in two levels, i.e., the inner maximization problem at the user level and the outer minimization problem at the operator level. In the inner layer, 𝐱 is optimized to maximize the Lagrangian function (<ref>) with given λ. The maximization of the Lagrangian function (<ref>) has the following explicit analytic solution:

    x_ik^* =  1,    if  i = max_j (b_j,k + ∑_zλ_zα_j,k^z)R_j,k
    
        0,    otherwise.

We can see that optimal x_ik^* can be determined by each user k independently according to (<ref>).
On the other hand, the dual problem (<ref>) is convex concerning λ and thus can be solved with the subgradient method. With x_ik^*, the dual variables can be updated as

    λ_z = [λ_z - η(U_z - U_z^0)]^+,

where x^+=max{x,0}. Specifically, λ can be updated by each operator independently. The above algorithm is compactly written in Algorithm <ref>








 §.§ Consensus ADMM-based Resourca Allocation and Beamforming Design


Due to the inter-operator interference (<ref>) and the MBC (<ref>), Algorithm <ref> is required to be implemented in a centralized manner, which is challenging to realize in practice. To this end, we propose an algorithm based on the consensus ADMM <cit.> to tackle this problem distributedly.  


First, P1.2 can be equivalently rewritten as follows:

  
    𝐏1.2.1:   max_𝐩,𝐭,𝐰̅,
    Γ̇,θ̇,ψ̇∑_z∈𝒪∑_i∈ℬ_z∑_k∈𝒰 x_i,kb_i,kr_i,k^(z)
       s.t. (<ref>), (<ref>), (<ref>), 
              r_i,k^(z) = ∑_n∈𝒪log_2(1+Γ_i,k^n,(z)), 
          ∑_i∈ℬ_z∑_k∈𝒰 x_i,kα_i,k^zr_i,k^(z)≥ U_z^0, 
          Γ_i,k^n,(z)≤ g_1(𝐰_i,k^n,β_i,k^n), 
          β_i,k^n ≥Ĩ_i,k^n + ∑_m∈𝒪∖{z}θ_m,k^n,(z) + σ_t^2,  
          ρ_i,k^n ≤∑_i∈ℬ_z∑_k^'∈𝒰∖{ k}g_2(𝐰_i,k^'^n) +∑_m∈𝒪∖{z}ψ_m,k^n,(z), 
    
              θ_z,k^n,(z)≥∑_i∈ℬ_z∑_k^'∈𝒰∖{ k}|𝐡_i,k^n𝐰_i,k^'^n| ^2, 
    
              ψ_z,k^n,(z)≤∑_i∈ℬ_z∑_k^'∈𝒰∖{ k}g_2(𝐰_i,k^'^n), 
          Γ_i,k^n,(z) = Γ̇_i,k^n, 
    
              θ_z,k^n,(m) = θ̇_z,k^n, 
    
              ψ_z,k^n,(m) = ψ̇_z,k^n,


where 𝐰̅≜{𝐰̅_z}_z∈𝒪, and 𝐰̅_z≜{𝐰_i,k^n, β_i,k^n,φ_i,k^n,ρ_i,k^n,Γ_i,k^n,(z),θ_i,k^n, (z),ψ_i,k^n, (z)}_n∈𝒪,i∈ℬ_z,k∈𝒰 represents the local variables of operator z. Γ_i,k^n,(z), θ_i,k^n, (z) and ψ_i,k^n, (z) are newly introduced slack variables.


Specifically, Γ_i,k^n,(z) represents the SINR of user k associated with BS/ST i over band n, while θ_i,k^n, (z) and ψ_i,k^n, (z) denote the upper and lower bound for the inter-operator interference that operator z causes to user k over band n, respectively. 
The reason we introduce θ_i,k^n,(z) is to guarantee MBCs (<ref>), while that for ψ_i,k^n,(z) is to satisfy backhaul constraints (<ref>). 

Note that Γ_i,k^n,(z), θ_i,k^n,(z) and ψ_i,k^n,(z) are all obtained and optimized by operator z locally. 
Intuitively, Γ_i,k^n,(z) is a local version of Γ_i,k^n determined by operator z, while θ_i,k^n,(z) and ψ_i,k^n,(z) are the local versions of I_j,k^n,j∈ℬ∖ℬ_z. 
On the other hand, Γ̇≜{Γ̇_i,k^n}, θ̇_𝐢,𝐤^𝐧≜{θ̇_i,k^n} and ψ̇_𝐢,𝐤^𝐧≜{ψ̇_i,k^n} are the corresponding global versions of Γ_i,k^n, θ_i,k^n and ψ_i,k^n, which can be obtained through aggregating local variables across operators. The constraints (<ref>)-(<ref>) can make sure that operators finally reach a consensus on the value of SINR and inter-operator interference, guaranteeing the equivalence of (<ref>) and (<ref>). 


With the introduced slack variables, 



the constraints from (<ref>) to (<ref>) can be decomposed into independent local convex sets for each operator, and the local set for operator z is:

    𝒞_z = {   (𝐩,𝐭,𝐰̅_z ) |(<ref>), (<ref>), (<ref>), (<ref>)-(<ref>)
          
          
          
          
          
          
        },


Then, the subproblem of (<ref>) for operator z can be compactly expressed as

  
    𝐏1.2.2:   min_𝐩,𝐭,𝐰̅_z,
    Γ̇,θ̇,ψ̇∑_i∈ℬ_z∑_k∈𝒰 -x_i,kb_i,kr_i,k^(z)
       s.t. (<ref>)
          (𝐩,𝐭,𝐰̅_z) ∈𝒞_z,
          Γ_i,k^n,(z) = Γ̇_i,k^n, ∀ i∈ℬ,∀ k∈𝒰,
          θ_z,k^n,(m) = θ̇_z,k^n, ∀ m ∈𝒪,∀ i∈ℬ,∀ k∈𝒰,
          ψ_z,k^n,(m) = ψ̇_z,k^n, ∀ m ∈𝒪,∀ i∈ℬ,∀ k∈𝒰.


For notational simplicity, define Φ_z≜{𝐩,𝐭,𝐰̅_z} and ϝ_z≜{Γ_i,k^n,(z)}_n∈𝒪,i∈ℬ_z,k∈𝒰. Similarly, we define θ_z and ψ_z. Then the augmented Lagrangian of (<ref>) can be written in the following

  
    L(Φ_z,Γ̇,θ̇,ψ̇,ν_z)=
           ∑_i∈ℬ_z∑_k∈𝒰 -x_i,kb_i,kr_i,k^(z) + 
        
          𝕀_𝒞_z(Φ_z)
        + ν_z, Γ^T(ϝ_z-Γ̇) + c/2‖ϝ_z-Γ̇‖^2
        + ν_z, θ^T(θ_z-θ̇) + c/2‖θ_z-θ̇‖^2
        + ν_z, ψ^T(ψ_z-ψ̇) + c/2‖ψ_z-ψ̇‖^2,


where 𝕀_𝒞_z(Φ_z) is an indicator function defined as

    𝕀_𝒞_z(Φ_z) = 
        0,    if Φ_z∈𝒞_z,
    
        +∞,    otherwise,

ν_z≜{ν_z,Γ,ν_z,θ,ν_z,ψ} denotes the dual variables related to constraints (<ref>)-(<ref>) of operator z, and c is the penalty coefficient introduced for penalizing the violation of equality constraints.

In the following, we elaborate on the variable updates in ADMM. First, the local variables {𝐩,𝐰̅_z} or {𝐭,𝐰̅_z} in the (l+1)-th iteration can be updated as follows:

    (𝐩,𝐰̅_z)^(l+1)   = min_𝐩,𝐰̅_zL(Φ^(l),Γ̇^(l),θ̇^(l),ψ̇^(l),ν_z^(l)),
    
        (𝐭,𝐰̅_z)^(l+1)   = min_𝐩,𝐰̅_zL(Φ^(l),Γ̇^(l),θ̇^(l),ψ̇^(l),ν_z^(l)).

Next, the global variables can be updated as

    Γ̇^(l+1)   = 1/|𝒪|∑_z∈𝒪(ϝ^(l+1)_z + 1/cν_z,Γ^(l)), 
    θ̇^(l+1)   = 1/|𝒪|∑_z∈𝒪(θ^(l+1)_z + 1/cν_z,θ^(l)),
    ψ̇^(l+1)   = 1/|𝒪|∑_z∈𝒪(ψ^(l+1)_z + 1/cν_z,ψ^(l)),

where |𝒪| is the number of operators. 

Finally, dual variables can be updated as follows:

    ν_z, Γ^(l+1)   = ν_z, Γ^(l) + c(ϝ_z^(l+1)-Γ̇^(l+1)),
    ν_z, θ^(l+1)   = ν_z, θ^(l) + c(θ_z^(l+1)-θ̇^(l+1)),
    ν_z, ψ^(l+1)   = ν_z, ψ^(l) + c(ψ_z^(l+1)-ψ̇^(l+1)).

The algorithm for updating {𝐩[τ], 𝐰̅[τ]} or {𝐭[τ], 𝐰̅[τ]} is summarized in Algorithm <ref>.




Motivated by the fact observed in <cit.>, we adopt the early termination strategy in the ADMM-loop. However, the new objective value may be worse than the previous one due to limited iterations. In this case, we will increase I_ADMM and not accept the new results.


The overall algorithm 

is summarized in Algorithm <ref>. 















 §.§ Convergence and Complexity Analysis


The convergence of Algorithm <ref> is provided as follows.

  Algorithm <ref> is guaranteed to converge.


  Please refer to Appendix <ref>

The complexity of Algorithm <ref> is dominated by Algorithm <ref> (line 5 of Algorithm <ref>) and Algorithm <ref> (line 3 and 4 of Algorithm <ref>). The complexity of Algorithm <ref> is 𝒪(I_UA(N_B+N_B^')(N_U+N_U^')), in which I_UA is the iteration number of Algorithm <ref>. The main complexity of Algorithm <ref> is solving (<ref>), which is similar to Algorithm <ref>. As a result, the complexities for line 3 and 4 of Algorithm <ref> are 𝒪(I_ADMM(N_L+(N_B+N_B^')(N_U+N_U^'))^4) and 𝒪(I_ADMM(N_B^'+(N_B+N_B^')(N_U+N_U^'))^4), respectively. 








§ SIMULATION RESULTS


In this section, simulation results are provided to demonstrate the effectiveness of the proposed algorithms and the significance of the MBCs. The considered SAGIN is shown in Fig. <ref>. Specifically, GNO has N_B=2 BSs and N_U=3 users, while SNO has N_B^'=4 STs and N_U^'=7 users. Each BS/ST is equipped with N_t=4 antennas, and the maximum transmit power for each BS is 52 dBm, and for each ST, it is 49 dBm. Besides, we assume each operator has a B_C = 100 MHz dedicated C band. The parameters of the satellite-ground link are summarized in Table <ref>.
The channel for the ground-ground link is set as follows: the path-loss is modeled as 32.4 + 20log(f_C) + 30log(d), where f_C=3 is the carrier frequency in GHz and d is the distance in km. The Rayleigh fading is adopted to characterize the small-scale fading, and the noise power spectral density is -174 dBm/Hz. Other parameters are set as follows: ϱ=1×10^-4, q=50, ϵ=20, η=0.1, ϵ=10^-3, c=1.5, and b_i,k=1, ∀ i,k.





To show the effectiveness of the proposed algorithms, we consider the following baseline algorithms for comparison.

  
  * Closest Association and Equal Resource Allocation (CA-ERA): Each user is connected to the closest BS/ST, with equal power allocation and maximum ratio transmit (MRT) beamfroming employed by BSs/STs. The LEO satellite equally allocates power to beams and time to STs within a beam.
  
  * Closest Association and Optimized Power and Beamforming (CA-OPW): The power allocation of the LEO satellite and the beamforming of BSs/STs are optimized with the proposed algorithm.
  
  * Closest Association and Optimized Time and Beamforming (CA-OTW): The time allocation of each beam and the beamforming of BSs/STs are optimized with the proposed algorithm.










Firstly, we validate the convergence of proposed Algorithms.

In Fig. <ref>, we show the convergence of the proposed centralized Algorithm <ref> (`Cen.') and distributed Algorithm <ref> (`Dis.') when δ_G=δ_S=0.6. It can be seen that both algorithms converge finally, and the distributed algorithm achieves a very close performance to the centralized algorithm.
Then, the convergence of the Algorithm <ref> is demonstrated in Fig. <ref>. 


Secondly, we compare WSRs achieved by different approaches versus the transmit power of the LEO satellite in Fig. <ref>. 
It is clear that WSRs increase with P_Sat increasing, indicating that the impact of backhaul capacity on the WSR.
Then, we compare WSRs achieved by Algorithm <ref> and CA-ERA for the case without inter-operator sharing (`w/o. Sharing'). The proposed algorithm outperforms CR-ERA, demonstrating the effectiveness of Algorithm <ref>.
Besides, we illustrate the WSRs for different approaches in the case with inter-operator sharing (`w/ Sharing'). The distributed algorithm approaches the centralized algorithm and performs much better than other benchmarks. This points out that the necessity of jointly optimizing user association, resource allocation  and beamforming design. 



Moreover, the WSR achieved in the case with inter-operator sharing is much higher than that in the case without inter-operator sharing, indicating significant gains introduced by realizing SAGIN with inter-operator sharing. 



Thirdly, we investigate the WSRs achieved by different approaches versus the transmit power of STs in Fig. <ref>.
It can be seen that the WSR in the case with inter-operator sharing does not keep increasing as P_i increases. Particularly, when P_i is greater than 40 dBm, the WSR almost remains the same. This indicates that the backhaul capacity provided by the LEO satellite limits the WSR, which is consistent with the result obtained from Fig. <ref>. Additionally, by comparing WSRs achieved by different approaches in two cases, we can find that the proposed algorithms always realize the highest WSRs and the distributed algorithm performs close to that of the centralized algorithm.





Next, we present the WSR and operators' revenue versus sharing coefficients to investigate the effect of MBCs in Fig. <ref> and Fig. <ref>. Specifically, we consider two cases: i) the symbiotic case with MBCs (`w/ MBC'), and ii) the non-symbiotic case without MBCs (`w/o. MBC'). 
In Fig. <ref>, we demonstrate WSRs and revenue versus δ_G with δ_S=1. 
By comparing revenue in the two cases, we can find that each operator's revenue in the symbiotic case does not experience revenue loss after sharing spectrum and services. On the contrary, the GNO suffers revenue loss in the non-symbiotic case when δ_G is less than 0.2. This indicates that MBCs can enable operators to achieve mutual benefits.

Then, we compare the WSRs in the two cases. The WSRs in the symbiotic case are less than those in the non-symbiotic case when δ_G<0.2.

The reason is that when δ_G<0.2, the spectrum and service configuration maximizing the WSR makes GNO suffer revenue loss. Therefore, the GNO does not adopt this configuration in the symbiotic case.

It is worth noting that the WSRs achieved in the symbiotic case are close to those in the non-symbiotic case when δ_G>0.2, which demonstrates the proposed algorithms can efficiently maximize WSR and achieve mutual benefits.


Besides, the WSR and revenue of GNO both grow gradually as δ_G increases, while the revenue of the SNO increases first and then decreases when δ_G>0.3. This trend can be explained as follows. 


With δ_G increasing, the GNO can obtain more revenue from serving users of SNO. As a result, the GNO demonstrates a growing preference for inter-operator sharing, increasing the WSR and revenue of operators. However when δ_G>0.3,

the GNO takes away the most revenue from serving SNO's users, resulting in revenue loss of SNO.

Lastly, we study the WSR and revenue versus the sharing coefficient of SNO δ_S in Fig. <ref>. It can be seen that the MBC for SNO cannot be satisfied until δ_S≥0.8 in the non-symbiotic case. Besides, the WSR and revenue in the symbiotic case almost remain the same when δ_S<0.8. The reason is that the SNO cannot benefit enough from inter-operator sharing and lacks the motivation to share spectrum and services.
On the other hand, by comparing Fig. <ref> and Fig. <ref>, it can be observed that the GNO will obtain benefit when (δ_G≥0.3, δ_S=1) while the SNO gains when (δ_S≥0.8, δ_G=1). This result indicates a principle to achieve mutual benefits when inter-operator sharing is enabled: the operator with more resources, like the SNO in this paper, should get more compensation from others.

On the contrary, the GNO, who contributes less to the SAGIN, should not require much compensation from others.



















§ CONCLUSION


In this paper, we have studied a SAGIN consisting of a GNO and SNO, and inter-operator sharing has been adopted to construct such a large-scale network. 

To fully leverage the advantage of inter-operator sharing, both system performance and individual benefits should be considered.

Therefore, we have investigated the spectrum and service sharing in the SAGIN from a symbiotic communication perspective, which 

can instruct different operators to achieve mutual benefits. Specifically, we have studied a WSR maximization problem about jointly optimizing the user association, resource allocation, and beamforming design. Besides, we have introduced a sharing coefficient to characterize the revenue of each operator, based on which the MBC has been formulated to guarantee that revenue for all operators does not decline after inter-operator sharing. Then, we have proposed a centralized algorithm based on SCA to solve the WSR maximization problem. Since implementing the centralized is difficult in real networks, we have also developed a distributed algorithm based on Lagrangian dual decomposition and consensus ADMM. Finally, extensive simulation results have been provided to show the effectiveness of the proposed algorithms and that the distributed algorithm can approach the centralized algorithm.
Moreover, the results have indicated that the MBCs can enable the operators to realize symbiotic communication in the SAGIN.




§ PROOF OF PROPOSITION <REF>

Let F(𝐩[τ-1],𝐭[τ-1], 𝐱[τ-1], 𝐰̃[τ-1]) denote the original objective, and F_𝐩[τ-1](𝐩[τ-1],𝐭[τ-1], 𝐱[τ-1], 𝐰̃[τ-1]) denote the approximated objective with the SCA method at point 𝐩[τ-1]. For notational simplicity, we use F(𝐩[τ-1],𝐰̃[τ-1]) and F_𝐩[τ-1](𝐩[τ-1],𝐰̃[τ-1]) in the following. In step 4 of Algorithm <ref>, we have the following inequality
  
    F(𝐩[τ-1],𝐰̃[τ-1])    (a)= F_𝐩[τ-1](𝐩[τ-1],𝐰̃[τ-1]) 
       (b)≤ F_𝐩[τ-1](𝐩[τ],𝐰̃[τ]) 
       (c)≤ F(𝐩[τ],𝐰̃[τ])

where (a) comes from the fact that the first-order Taylor expressions used in (<ref>), (<ref>), and (<ref>) are tight at point 𝐩[τ-1]; (b) holds since (𝐩[τ],𝐰̃[τ]) are the optimal solutions to (<ref>) with given (𝐭[τ-1],𝐱[τ-1]); and (c) holds since F_𝐩[τ-1](𝐩[τ],𝐰̃[τ]) is a lower bound of F(𝐩[τ],𝐰̃[τ]).

Note that F_𝐩[τ](𝐩[τ],𝐰̃[τ]) is the input of step 5, and the inequality (<ref>) can be applied to analyzing step 5 and 6. Finally, we have F(𝐩[τ-1],𝐭[τ-1], 𝐱[τ-1], 𝐰̃[τ-1]) ≤ F(𝐩[τ],𝐭[τ], 𝐱[τ], 𝐰̃[τ]), which implies that Algorithm <ref> can generate a non-decreasing objective after each iteration. Moreover, the objective is constrained due to the power budget of BSs/STs and the satellite. As a result, Algorithm <ref> is guaranteed to converge.



§ PROOF OF PROPOSITION <REF>

In step 3 of Algorithm <ref>, we have the following inequality
  
    F_𝐩[τ-1](𝐩[τ-1],𝐰̃[τ-1])    (a)≤ F_𝐩[τ-1](𝐩[τ],𝐰̃[τ])

where (a) comes from the fact that the subproblem (<ref>) is convex with fixed (𝐭,𝐱) and consensus ADMM can converge to the optimal solution with enough iterations <cit.>. Besides, it is clear that (a) also holds when early termination strategy is adopted due to the Algorithm <ref> (step 7). The inequality (<ref>) is also applicable for step 4 of Algorithm <ref>. As a result, we have F(𝐩[τ-1],𝐭[τ-1], 𝐱[τ-1], 𝐰̃[τ-1]) ≤ F(𝐩[τ],𝐭[τ], 𝐱[τ-1], 𝐰̃[τ]).

In step 6 of Algorithm <ref>, we have the following inequality

    F_𝐱[τ-1](𝐱[τ-1])    (b)≤ F_𝐱[τ-1](𝐱[τ])≤ F(𝐱[τ])

where (b) holds since the subgradient method is guaranteed to converge to the optimal solution when the dual problem (<ref>) is convex <cit.>. As a result, we have F(𝐩[τ-1],𝐭[τ-1], 𝐱[τ-1], 𝐰̃[τ-1]) ≤ F(𝐩[τ],𝐭[τ], 𝐱[τ], 𝐰̃[τ]), which implies that Algorithm <ref> generates a non-decreasing objective after each iteration. Therefore, Algorithm <ref> is guaranteed to converge. 






IEEEtran ]