Optimal Spectrum Partitioning and Licensing in Tiered Access under Stochastic Market Models
OOptimal Spectrum Partitioning and Licensing inTiered Access under Stochastic Market Models
Gourav Saha, and Alhussein A. Abouzeid
Abstract —We consider the problem of partitioning a spectrumband into M channels of equal bandwidth, and then furtherassigning these M channels into P licensed channels and M − P unlicensed channels. Licensed channels can be accessed both forlicensed and opportunistic use following a tiered structure whichhas a higher priority for licensed use. Unlicensed channels canbe accessed only for opportunistic use. We address the followingquestion in this paper. Given a market setup, what values of M and P maximize the net spectrum utilization of the spectrumband? While this problem is of fundamental nature, it is highlyrelevant practically, e.g., in the context of partitioning the recentlyproposed Citizens Broadband Radio Service band. If M is toohigh or too low, it may decrease spectrum utilization due tolimited channel capacity or due to wastage of channel capacity,respectively. If P is too high (low), it will not incentivize thewireless operators who are primarily interested in unlicensedchannels (licensed channels) to join the market. These tradeoffsare captured in our optimization problem which manifests itselfas a two-stage Stackelberg game. We design an algorithm tosolve the Stackelberg game and hence find the optimal M and P . The algorithm design also involves an efficient Monte Carlointegrator to evaluate the expected value of the involved randomvariables like spectrum utilization and operators’ revenue. Wealso benchmark our algorithms using numerical simulations. Index Terms —Spectrum License, Opportunistic Spectrum Ac-cess, CBRS band, Stackelberg game, Iterated Removal of StrictlyDominated Strategies, Monte Carlo Integration, Optimization
I. I
NTRODUCTION
To support the ever growing wireless data traffic, the FederalCommunication Commission (FCC) released the underutilizedCitizens Broadband Radio Service (CBRS) band for shared usein 2015 [1]. CBRS band is a
M Hz federal spectrum bandfrom . GHz to . GHz . The
M Hz band is dividedinto channels of M Hz each. The shared use of theCBRS band follows an order of priority. Federal users have thehighest priority access to the channels. Out of the channels, are Priority Access Licenses (PALs). PAL licenses are soldthrough auctions and the lease duration of a PAL license mayrange between − years [1], [2], [3]. A PAL license holdercan use their channel only if federal users are not using it.The remaining out of the channels are reserved only foropportunistic use by General Authorized Access (GAA) users.Opportunistic channel allocation to GAA users can happen ata time scale of minutes to weeks. GAA users can use these channels if federal users are not using the channels. GAAusers can also use the PAL channels provided that neitherfederal users nor PAL license holders are using it.
G. Saha and A. A. Abouzeid are with the Department of Electrical, Com-puter and Systems Engineering, Rensselaer Polytechnic Institute, Troy, NY12180, USA; Email: [email protected], [email protected]. A preliminaryversion of this paper has been accepted for publication in WiOpt 2020.
Opportunistic Use D ec r ea s i ng P r i o r it y Licensed Use Net BandwidthUnlicensed channels, M − PLicensed Channels, P (similar to PAL channels) (similar to channels for GAA access)
Opportunistic Use
Figure 1. Pictorial representation of the tiered spectrum model underconsideration in this paper.
As mentioned in the previous paragraph, the CBRS band isdivided into M = 15 channels out of which there are P = 7 PAL licenses. But does M = 15 and P = 7 maximize the uti-lization of the CBRS band? In this paper, we are interested inthe following abstraction of this question whose application isnot limited to CBRS band. A net bandwidth is partitioned into M channels of equal bandwidth. These M channels are furtherdivided into P licensed channels (similar to PAL channels)and M − P unlicensed channels (similar to channels reservedfor GAA users). In this paper, the process of dividing the netbandwidth into M channels is called spectrum partitioning and the process of allocating these M channels as licensedand unlicensed channels is called spectrum licensing . Licensedchannels are used for both licensed use and opportunistic usewith the former having higher priority. Unlicensed channelsare reserved for opportunistic use only. This spectrum accessmodel is shown in Figure 1. The wireless operators earnrevenue by serving customer demands. A wireless operator isincentivized to join the market if the revenue which it can earnis above a desired threshold. For the given setup, what valueof M and P maximizes spectrum utilization where spectrumutilization is defined as the net amount of customer demandserved by the entire bandwidth?There are various factors that decide the optimal values of M and P . Some of these factors are as follows. If the numberof channels M increases, the bandwidth, and hence capacity,of each channel decreases. The capacity of each channelshould be large enough to accommodate a good portion of thecustomer demand of a wireless operator but not so large thatmost of the capacity of the channel is not utilized for majorityof the time. This suggests that M should not be too small ortoo large. If the number of licensed channels P is too high,there is a small number of unlicensed channels. Therefore,those operators who primarily rely on unlicensed channels toserve customer demands will not be able to generate enoughrevenue and hence will not be incentivized to join the market.Similarly, if the P is too low, wireless operators who primarilyrely on licensed channels to serve customer demands willnot be incentivized to join the market. P should be set suchthat enough operators join the market to ensure that the a r X i v : . [ c s . G T ] F e b ustomer demands served over the entire bandwidth is as highas possible. There may be other qualitative factors governingoptimal M and P . Therefore, in this paper, we design analgorithm to jointly optimize M and P such that spectrumutilization is maximized. A. Related Work
Variations of the spectrum partitioning and spectrum li-censing problems considered in this paper have been studiedseparately, but not jointly, in the spectrum sharing and relatedfields. There are a few works that have addressed problemssimilar to partitioning of a fixed bandwidth into optimalnumber of channels. In [4], the authors derive an analyticalexpression for the optimal number of channels such that thespatial density of transmission is maximized subject to a fixedlink transmission rate and packet error rate. Partitioning ofbandwidth in the presence of guard bands has been consideredin [5] where the authors used Stackelberg game formulation toanalyze how a spectrum holder should partition its bandwidthin order to maximize its revenue in spectrum auctions.The second problem studied in this paper essentially dealswith spectrum licensing. This has been widely studied in theliterature from various perspectives. Some works concentratedon minimizing the amount of bandwidth allocated to backupchannels (unlicensed channels in our case) while providing acertain level of guarantee to secondary users against channelpreemption [6]. There has also been research on overlay D2Dand cellular devices which studied optimal partitioning oforthogonal in-band spectrum to maximize the average through-out rates of cellular and D2D devices [7]. In [8], the authorsinvestigated whether to allocate an additional spectrum bandfor licensed or unlicensed use and concluded that the licenseduse is more favourable for maximizing the social surplus. Asimilar result has been shown in [9] which studied the effect ofadding an unlicensed spectrum band in a market consisting ofwireless operators with licensed channels. The authors showedthat if the amount of unlicensed spectrum band is below acertain limit, the the overall social welfare may decrease withthe increase in unlicensed spectrum band. The authors in [10],[11] studied the CBRS band for a market setup that consistsof Environmental Sensing Capability operators (ESCs) whosesole job is to monitor and report spectrum occupancy to thewireless operators. The authors analyzed how the ratio of thelicensed and unlicensed bands affects the market competitionbetween the ESC operators, the wireless operators, and theend users of the CBRS band. There is a line of work whichstudies spectrum partitioning for topics similar to licensed andunlicensed use using Stackelberg games; macro cells and smallcells [12], [13], long-term leasing market and short-term rentalmarket [14], and 4G cellular and Super Wifi services [15].Such a diverse body of work just on spectrum partitioningand licensing is justified because individual problem setupshave their own salient features and hence require their ownanalysis. Our problem setup considers jointly optimizing spec-trum partitioning and spectrum licensing, which has not beenconsidered in the existing literature. This problem is novelbecause of the combination of the following two reasons.
First ,our spectrum access model, like CBRS, is a combination of (a) Unlicensed spectrum access model. This is because M − P unlicensed channels are reserved specifically for opportunisticuse. (b) Primary-secondary spectrum access model. This isbecause P licensed channels can be used for opportunisticaccess following the priority hierarchy. Prior works like [9],[10], [11], which solved the spectrum licensing problem,did not simultaneously consider both of the spectrum accessmodels. Second , we consider a very generalized system modelin terms of the number of operators, their types, and theirheterogeneity. Such a setup leads to a scenario where theregulator has to decide M and P such that the right set ofwireless operators are incentivized to join the market. B. Contribution and Paper Organization
We now present an overall outline of the paper and, inthe process, discuss its main contributions. In Section II, wepresent a system model which can mathematically capture theeffect of the number of channels, M , and number of licensedchannels, P , on the spectrum utilization. The proposed systemmodel captures spectrum auctions using a simple stochasticmodel without going into complex game-theoretic formula-tions. Based on our reading of [1], [2], [3] and other literatureon CBRS band, it is not clear if there is a consensus in theliterature/policy about whether PAL license holders are alsoallowed to use channels opportunistically. So it is possiblethat PAL license holders may or may not be allowed touse channels opportunistically. Our model is general enoughto capture both of these cases. Our model also generalizeswell to various opportunistic access strategies like overlay orinterweave spectrum access [16].It is possible that a choice of values of M and P thatincentivizes one group of wireless operators may not incen-tivize another group. Therefore, a M and P which incentivizesall the wireless operators may not exist. This argument canbe exemplified by refering to [1], [2], [3] which shows alot of debate between the wireless operators concerning theparameters of the CBRS model. Even if it is possible to satisfyall the operators, it may not be optimal to do so in termsof maximizing the spectrum utilization. We capture this ideaby formulating our problem as a two-stage Stackelberg gamein Section III-A which forms the second contribution of thepaper. The Stackelberg game consists of the regulator (leader)and the wireless operators (followers). In Stage 1, the regulatorsets M and P to maximize spectrum utilization. In Stage 2,the wireless operators decide whether or not to join the marketbased on the M and P set by the regulator in Stage 1.In Section III-B, we design an algorithm to solve theStackelberg game and hence the optimal M and P whichmaximize spectrum utilization. We approach this in steps. Fewproperties associated with expected revenue of an operator arediscussed first. We show that when these properties hold, wecan design a polynomial time algorithm to solve Stage 2 of theStackelberg game. We finally solve Stage 1 of the Stackelberggame using a grid search approach to find the optimal M and P which maximizes spectrum utilization. To the best ofour knowledge, joint optimization of partitioning and tieredlicensing have not been considered in the existing relatediterature.
Hence, designing an algorithm for joint optimizationof M and P is the fourth contribution of the paper.The solve the Stackelberg game, we have to calculate theexpected revenue of an operator and expected spectrum utiliza-tion. The complex nature of the problem does not allow simpleanalytical formulas of these expected values. Even if suchanalytical formulas are possible, adapting them to changes insystem model can be time consuming. Therefore, we developa Monte Carlo integrator to evaluate these expected valuesin Section IV. Our choice of using a Monte Carlo integratorover deterministic numerical integration techniques is becauseour setup involves evaluation of high-dimensional integrals.Unlike deterministic numerical integration techniques, thecomputation time of Monte carlo integration does not scalewith dimension. One of the main bottlenecks of Monte Carlointegration is random sampling. While designing our MonteCarlo integrator, we reduced random sampling as much aspossible to make it more time efficient. Designing an efficientMonte Carlo integrator which can easily adapt to few changesin the system model is the third contribution of the paper.Finally, we use the algorithms designed in Sections III-Band IV to obtain important numerical results in Section Vwhich constitutes the final contribution of the paper. Ournumerical results show how optimal values of M and P varywith market parameters.II. S YSTEM M ODEL
In this section, we discuss individual components of oursystem model in Sections II-A to II-C. The list of importantnotations is included in Table I. We introduce three settheoretic notations. Consider two sets A and B . The operation A (cid:83) B implies the union A and B . The operation A\B is asset which consists of all those elements in A which are not in B . A singleton set consisting of element a is denoted by { a } . A. Channel Model
A net bandwidth of
W hz is divided into M channels ofequal bandwidth WM . Out of the M channels, P channels arelicensed channels while the remaining M − P channels areunlicensed channels. In our model, time is divided into slotswhere t ∈ Z + denotes the t th time slot. Licensed channelsare allocated for prioritized licensed use and opportunistic usewhile unlicensed channels are allocated only for opportunisticuse. Allocation of licensed channels for licensed use happensthrough auctions. These auctions occur every T ≥ T is the lease duration. An entire lease duration is calledan “epoch”. Epoch γ is from time slot ( γ − T + 1 to γT .Allocation of licensed channels and unlicensed channels foropportunistic use occur every time slot.Those operators who are allocated licensed channels forlicensed use are called Tier-1 operators while those who arenot are called
Tier-2 operators . In our model, an operatorcan be allocated atmost one licensed channel for licensed usein an epoch, i.e. spectrum cap is one. Similar assumptionshas been made in prior works like [17]. Spectrum cap ofone ensures fairness by allocating the licensed channels toas many operators as possible. A Tier-1 operator can alsouse opportunistic channels to serve its customer demand in
Table IA
TABLE OF IMPORTANT NOTATIONS . Notation Description t , γ t th time slot and epoch γ resp. M , P Number of channels and number of licensed channels resp. φ Indicator variable which decides if a Tier-1 operator canalso use channels opportunistically. We have, φ ∈ { , } . D Capacity of the entire bandwidth for licensed access. α L , α U Interference parameter associated with opportunistic use oflicensed and unlicensed channels resp. S CL , S CU Set of candidate licensed and unlicensed operators resp. S L , S U Set of interested licensed and unlicensed operators resp. S Set of interested operators; S = S L (cid:83) S U . T ( γ ) Set of Tier-1 operators in epoch γ , i.e. set of interestedlicensed operators who won licensed channels in epoch γ . T ( γ ) Set of Tier-2 operators in epoch γ ; T ( γ ) = S\T ( γ ) . x k ( t ) Customer demand of the k th operator in t th time slot. µ θk , σ θk Mean and standard deviation resp. of Gaussian randomvariable θ k ( t ) where, x k ( t ) = max (0 , θ k ( t )) . (cid:101) x k,i ( t ) Amount of demand served by the k th operator in t th timeslot if it is a Tier- i operator where i ∈ { , } . (cid:101) x k,a ( t ) Amount of demand served by the k th operator in t th timeslot when spectrum access type is a where a ∈ { lc, op } . X k,a ( γ ) Net demand served by the k th operator in epoch γ whenspectrum access type is a where a ∈ { lc, op } . R k,a ( γ ) Revenue of k th operator in epoch γ using spectrum accessof type a where a ∈ { lc, op } . R k,i ( γ ) Revenue of k th operator in epoch γ if it is a Tier- i operator. V k ( γ ) Bid of the k th operator, where k ∈ S L , in epoch γ . h k ( · ) A function associated with k th operator which maps themean of X k,a ( γ ) to the mean of R k,a ( γ ) . µ Xk,a , σ
Xk,a
Mean and standard deviation of X k,a ( γ ) . µ Rk,a , σ
Rk,a
Mean and standard deviation of R k,a ( γ ) . We have, µ Rk,a = h k (cid:16) µ Xk,a (cid:17) . ρ k Correlation coefficent between X k,a ( γ ) and R k,a ( γ ) . ω k Correlation coefficent between V k ( γ ) and R k,lc ( γ ) . λ k Minimum revenue requirement of the k th operator. ξ k The tuple (cid:16) µ θk , σ θk , h k ( · ) , σ Rk,a , ρ k , ρ k , λ k (cid:17) associ-ated with the k th operator. ( x ) + ( x ) + = max (0 , x ) case the bandwidth of the allocated licensed channel is notsufficient. Tier-2 operators uses channels opportunsitically.Tier-1 operators may also use channels opportunistically. Let φ ∈ { , } be an indicator variable which decides if Tier-1 operators can use channels opportunistically. If φ = 1 ,then Tier-1 operators can use channels opportunistically andotherwise they cannot.The capacity of a channel/bandwidth is the maximum unitsof customer demand that can be served using that chan-nel/bandwidth in a time slot. Let D be the capacity of theentire bandwidth of W hz when used for licensed access.As the entire bandwidth is partitioned into M channels, eachlicensed channel has a capacity of DM when used for licenseduse while the unlicensed channels has a capacity of α U DM where α U ∈ [0 , is the interference parameter associated withunlicensed channels for opportunistic use. Licensed channelscan also be used for opportunistic use following the priorityhierarchy shown in Figure 1. Let the customer demand of aTier-1 operator be d units. It will use its licensed channelto serve its customer demand. The remaining capacity of thelicensed channel which can be utilized for opportunsitic useis given by the function C ( d, α L ) where α L ∈ [0 , is theinterference parameter associated with licensed channels forpportunistic use. The expression for C ( d, α L ) depends onthe opportunistic spectrum access (OSA) strategy: overlay orinterweave [16]. For overlay spectrum access , C ( d, α L ) = α L (cid:0) DM − d (cid:1) + , where ( x ) + = max (0 , x ) . For interweavespectrum access , C ( d, α L ) is equal to if d > and equalto α L DM if d = 0 . Parameters α L and α U capture the lowerefficiency of opportunistic use as compared to licensed use [1].In general, we expect α L ≤ α U . This may happen because thetransmission power cap for opportunistic use maybe lower forlicensed channels compared to unlicensed channels in order toprotect Tier-1 operators from harmful interference. B. Operators, their Demand and Revenue Model
The market consists of the candidate licensed operators denoted by S CL and the candidate unlicensed operators de-noted by S CU where S CL and S CU are disjoint sets. A candidatelicensed operator is a Tier-1 operator in those epochs in whichit is allocated a licensed channel in the auction and a Tier-2operator in those epochs in which it is not allocated a licensedchannel. A candidate unlicensed operator is always a Tier-2operator. A candidate operator has to invest in infrastructuredevelopment if it wants to join the market.All the candidate operators have to invest in infrastructuredevelopment to join the market. In order to generate returnon infrastructure cost and the cost of leasing a channel, acandidate operator wants to earn a minimum expected revenuein an epoch. Let λ k be the minimum expected revenue (MER)of the k th operator. A subset of candidate licensed andunlicensed operators are interested in joining the market ifthe value of M and P set by the regulator is such that theexpected revenue of the operator in an epoch is greater than itsMER. The set of interested licensed operators and interestedunlicensed operators are denoted by S L and S U respectively.We have S L ⊆ S CL and S U ⊆ S CU . The set of operators, (cid:0) S CL − S L (cid:1) (cid:83) (cid:0) S CU − S U (cid:1) , does not join the market. A can-didate licensed/unlicensed operator gets to decide whether tojoin or not join the market only once. An operator gets toparticipate in auctions for licensed channels or to use channelsopportunistically only if it decides join the market.In our model, every operator has a separate pool of cus-tomers each with its own stochastic demands, i.e. we donot model price competition between operators to attract acommon pool of customers. Consider the t th time slot ofepoch γ . The customer demand, or simply demand, of the k th operator in the t th time slot is x k ( t ) . In our model, x k ( t ) = max (0 , θ k ( t )) where θ k ( t ) are iid Gaussian ran-dom variable with mean µ θk and standard deviation σ θk , i.e. θ k ( t ) ∼ N (cid:16) µ θk , (cid:0) σ θk (cid:1) (cid:17) , ∀ t . The k th operator may be able toserve only a fraction of the customer demand. Let (cid:101) x k, ( t ) and (cid:101) x k, ( t ) denote the amount of customer demand served by the k th operator if it is a Tier-1 and a Tier-2 operator respectivelyin epoch γ . We have (cid:101) x k, ( t ) ≤ x k ( t ) and (cid:101) x k, ( t ) ≤ x k ( t ) . (cid:101) x k, ( t ) and (cid:101) x k, ( t ) can be expressed as follows (cid:101) x k, ( t ) = (cid:101) x k,lc ( t ) + (cid:101) x k,op ( t ) (1) (cid:101) x k, ( t ) = (cid:101) x k,op ( t ) (2) All the iid random variables used throughout the paper are identical withrespect to time slot, t , or epoch, γ , and not with respect to operator index k . where (cid:101) x k,lc ( t ) = min (cid:0) x k ( t ) , DM (cid:1) . The term (cid:101) x k,lc ( t ) in (1) isthe amount of customer demand served a Tier-1 operator usingthe channel allocated to it for licensed use. The term (cid:101) x k,op ( t ) in (1) and (2) is the demand served by an operator by usingchannels opportunistically. It will be shown in Section II-Cthat (cid:101) x k,op ( t ) is a iid random variable. Also, if φ = 0 , thena Tier-1 operator cannot use channels opportunistically andhence (cid:101) x k,op ( t ) ≡ in (1). In (1) and (2), (cid:101) x k, ( t ) and (cid:101) x k, ( t ) can be expressed as a time invariant function of iid randomvariables x k ( t ) and (cid:101) x k,op ( t ) . Therefore, (cid:101) x k, ( t ) and (cid:101) x k, ( t ) are iid random variables as well. Let (cid:101) µ xk,a and (cid:101) σ xk,a denote themean and standard deviation of (cid:101) x k,a ( t ) respectively. We have, (cid:101) µ xk,lc = DM (cid:90) ϑf θk ( ϑ ) dϑ + DM ∞ (cid:90) DM f θk ( ϑ ) dϑ (3) (cid:0)(cid:101) σ xk,lc (cid:1) = DM (cid:90) ϑ f θk ( ϑ ) dϑ + (cid:18) DM (cid:19) ∞ (cid:90) DM f θk ( ϑ ) dϑ − (cid:0)(cid:101) µ xk,lc (cid:1) (4)where f θk ( ϑ ) is the probability density function of θ k ( t ) .In general, an analytical expression for (cid:101) µ xk,op and (cid:101) σ xk,op isnot possible because of the complex nature of opportunisticspectrum allocation algorithm. We have designed a MonteCarlo integrator which can compute (cid:101) µ xk,op in Section IV.Throughout the rest of the paper we will use the subscript k, i , where i ∈ { , } , to denote variables associated with k th operator when its is a Tier- i operator. Also, we will usethe subscript k, a , where a ∈ { lc, op } , to denote variablesassociated with k th operator when access type is licensed ( a = lc ) or opportunistic ( a = op ). Let X k,a ( γ ) denote the netdemand served by the k th operator in epoch γ when accesstype is a . Mathematically, X k,a ( γ ) = γT (cid:88) t =( γ − T +1 (cid:101) x k,a ( t ) ; a ∈ { lc, op } (5)Since (cid:101) x k,a ( t ) is iid random variable and the lease duration T is quite large in practice, X k,a ( γ ) can be approximated asa Gaussian random variable using Central Limit Theorem [18,Chapter 8]. The mean, µ Xk,a , and standard deviation, σ Xk,a , of X k,a ( γ ) are given by µ Xk,a = (cid:101) µ xk,a T ; σ Xk,a = (cid:101) σ xk,a √ T (6)To this end we have, X k,a ( γ ) ∼ N (cid:18) µ Xk,a , (cid:16) σ Xk,a (cid:17) (cid:19) , ∀ γ . Remark 1: Gaussian nature of X k,a ( γ ) . X k,a ( γ ) is alwaysa positive quantity because net demand served is alwayspositive. But, we approximated X k,a ( γ ) as a Gaussian randomvariable and hence the approximated X k,a ( γ ) can be negative.However, the probability of X k,a ( γ ) being negative is P [ X k,a ( γ ) <
0] = 12 (cid:32) erf (cid:32) − (cid:101) µ xk,a √ T √ (cid:101) σ xk,a (cid:33)(cid:33) where erf ( · ) is the error function. For all practical setup, T islarge enough that P [ X k,a ( γ ) < is very small. The use ofGaussian model for non-negative random variables have beenused in prior works like [19].n operator generates revenue by serving customer demand.Let R k,a ( γ ) denote the revenue earned by the k th operatorin epoch γ when access type is a . We model R k,a ( γ ) as arandom variable which follows the stochastic model (cid:20) X k,a ( γ ) R k,a ( γ ) (cid:21) ∼ N (cid:20) µ Xk,a µ Rk,a (cid:21) , (cid:16) σ Xk,a (cid:17) ρ k σ Xk,a σ Rk,a ρ k σ Xk,a σ Rk,a (cid:16) σ Rk,a (cid:17) (7)for all γ where µ Rk,a = h k (cid:16) µ Xk,a (cid:17) . According to (7), the netdemand served and the net revenue earned in epoch γ arejointly Gaussian. The mean of R k,a ( γ ) is µ Rk,a = h k (cid:16) µ Xk,a (cid:17) where h k (cid:16) µ Xk,a (cid:17) is a monotonic increasing function of themean demand served by the k th operator in an epoch, µ Xk,a .The standard deviation of R k,a ( γ ) is σ Rk,a which can be usedto capture the effect of exogeneous stochastic processes likemarket dynamics on R k,a ( γ ) . The relative change between R k,a ( γ ) and X k,a ( γ ) is captured with correlation coefficent ρ k ∈ [0 , . It captures how much a deviation of X k,a ( γ ) around its mean µ Xk,a will effect the deviation of R k,a ( γ ) around its mean h k (cid:16) µ Xk,a (cid:17) . A monotonic increasing function, h k ( · ) , and a positive correlation coefficient, ρ k , are intuitivebecause from a statistical standpoint it implies that an operatorwho serves more customer demand generates higher revenue.Let R k,i ( γ ) denote the revenue earned by the k th operatorif it is a Tier- i operator in epoch γ . Tier-1 serves customerdemand using both licensed and opportunistic access whileTier-2 operator serves its customer demand using opportunisticaccess only. Hence, R k, ( γ ) = R k,lc ( γ ) + R k,op ( γ ) (8) R k, ( γ ) = R k,op ( γ ) (9)Notice that since R k,lc ( γ ) and R k,op ( γ ) are Gaussian, R k, ( γ ) and R k, ( γ ) are Gaussian as well. C. Spectrum Allocation Model
Licensed channels are allocated to the set of interested li-censed osperators, S L , through spectrum auctions. The auctionfor epoch γ happens at time slot ( γ − T + 1 . The set ofinterested licensed operators bids for licensed channels. Let V k ( γ ) be the bid of the k th operator in epoch γ . Intuitively,the bid of the k th operator should depend on its valuation ofa licensed channel. The value of a licensed channel to the k th operator in epoch γ is R k,lc ( γ ) , the revenue it can generatein an epoch using the licensed channel. We capture thisdependence between V k ( γ ) and R k,lc ( γ ) using a correlationcoefficient ω k ∈ [0 , between them. In our model, V k ( γ ) and R k,lc ( γ ) are jointly Gaussian. We have, (cid:20) V k ( γ ) R k,lc ( γ ) (cid:21) ∼ N (cid:20) µ Rk,lc µ Rk,lc (cid:21) , (cid:16) σ Rk,lc (cid:17) ω k (cid:16) σ Rk,lc (cid:17) ω k (cid:16) σ Rk,lc (cid:17) (cid:16) σ Rk,lc (cid:17) (10)for all γ . In (10), µ Rk,lc = h k (cid:16)(cid:101) µ xk,lc T (cid:17) (refer to (6) and (7)).Using a stochastic model like (10) to capture the relationbetween V k ( γ ) and R k,lc ( γ ) leads to a more generic system Spectrum Auctionfor Epoch 2Spectrum Auctionfor Epoch 1
Epoch 1 Epoch 2
10 11 S CL S L S U S CU S U S U T (1) T (2) T (1) T ( ) T ( ) T (2) Figure 2. Pictorial representation of the set of candiate licensed operators S CL ,candiate unlicensed operators S CU , interested licensed operators S L , interestedunlicensed operators S U , the set of interested licensed operators who areallocated (not allocated) licensed channels in an epoch T ( γ ) ( T ( γ ) ) andthe set of Tier-2 operators in an epoch T ( γ ) . Note that T ( γ ) , T ( γ ) and T ( γ ) are not same for epochs and . model because we can abstract away from the exact bidestimation strategies of the operators which may rely onvarious market externalities.Given that there are P licensed channels and the spectrumcap is one, the interested licensed operators with the P highestbids V k ( γ ) are allocated one licensed channel each in epoch γ . The operators who are allocated a licensed channel haveto pay a price that is determined by the regulator. Sucha pricing model has been used in prior works [15]. Let T ( γ ) ⊆ S L denote the set of interested licensed operatorswho are allocated licensed channels in epoch γ . Similarly, T ( γ ) = S L \T ( γ ) are the set of interested licensed operatorswho are not allocated licensed channels in epoch γ . Theoperators in T ( γ ) serve their customer demand as Tier-1operators in epoch γ . On the other hand, operators in T ( γ ) serve their customer demand as Tier-2 operators in epoch γ . Itis to be noted that T ( γ ) and T ( γ ) are random sets as theyget decided by the bids V k ( γ ) which are random variables.The set of Tier-2 operators in epoch γ is T ( γ ) = T ( γ ) (cid:83) S U ,i.e., interested unlicensed operators and interested licensedoperators who are not allocated a licensed channel in epoch γ . Unlike the sets S L and S U which are decided once, sets T ( γ ) , T ( γ ) , and T ( γ ) are decided in the beginning ofevery epoch. A pictorial representation of all the importantsets discussed till now is shown in Figure 2. Figure 2 alsoshows T ( γ ) , T ( γ ) , and T ( γ ) varies with epoch γ .Opportunistic spectrum allocation happens in every timeslot to all the Tier-2 operators. Tier-1 operators may or maynot participate in opportunistic spectrum access depending onwhether φ is one or zero. In order to capture both thesecases under a single mathematical abstraction, we modifythe demand of Tier-1 and Tier-2 operators. Let x k ( t ) be themodified demand of the k th operator which needs to be servedusing OSA. For time slot t of epoch γ , x k ( t ) = (cid:40) φ · (cid:0) x k ( t ) − DM (cid:1) + ; k ∈ T ( γ ) x k ( t ) ; k ∈ T ( γ ) (11)ccording to (11), for a Tier-2 operator, its entire demand x k ( t ) needs to be served using OSA. For Tier-1 operators, theexcess demand which could not be satisfied with licensed useis (cid:0) x k ( t ) − DM (cid:1) + . If φ = 1 , this excess demand has to servedusing OSA. If φ = 0 , then x k ( t ) = 0 for Tier-1 operatorsimplying that they don’t participate in OSA.Opportunistic channel capacity in time slot t of epoch γ is D O ( t ) = α U (cid:16) M − (cid:101) P (cid:17) DM + (cid:80) k ∈T ( γ ) C ( x k ( t ) , α L ) (12)where (cid:101) P = min ( |S L | , P ) . In (12), the first term is the netchannel capacity of unlicensed channels and the second termis the net remaining channel capacity of the licensed channels.The variable (cid:101) P is used to capture edge cases where thenumber of licensed channels is more than number of interestedlicensed operators. In such cases, the remaining P − |S L | channels which are not allocated to licensed operators areused as unlicensed channels. The expression for C ( x k ( t ) , α L ) depends on the OSA strategy (overlay or interweave) and hasbeen discussed in Section II-A. As our model is inspired bythe CBRS band, we have to ensure that opportunistic spectrumallocation is fair [20]. One approach to ensure fair allocationand to avoid wastage of channel capacity is to use a max-min fair algorithm, like the famous Waterfilling algorithm. Adetailed exposition of max-min fairness can be found in [21],[22]. In this section, we present the Waterfilling algorithm,explain its working with an example and qualitatively justifyhow it ensures fairness and avoids wastage of channel capacity.Waterfilling algorithm will be used for opportunistic channelallocation throughout this paper.Algorithm 1 is the psuedocode of Waterfilling algorithm.Let S denote the set of interested operators, i.e. S = S L (cid:83) S C .The union of Tier-1 and Tier-2 operators, T ( γ ) (cid:83) T ( γ ) , isequal to S . The inputs to Algorithm 1 are the opportunisticchannel capacity, D O ( t ) , and the modified demands of all theinterested operators, { x k ( t ) } k ∈S . The output of Algorithm1 is the opportunistic channel capacity allocated to all theinterested operators, { (cid:101) x k,op ( t ) } k ∈S . (cid:101) x k,op ( t ) is also equal tothe demand served by the operators using OSA. We use thefollowing example to explain Algorithm 1: the set of interestedoperators is S = { , , , , } , their corresponding modifieddemand is { , , , , } , and the opportunistic channel capac-ity D O ( t ) = 17 . The example is shown in Figure 3.Waterfilling algorithm allocates channel capacity to the setof interested operators in ascending order of their modifieddemand ( lines 1 and ). The sorted list of modified demandis { , , , , } and the operator index κ ( j ) correspondingto position j = 1 , , , , of the sorted list is , , , , respectively. In line 2 , unused opportunistic channelcapacity C = 17 and the remaining number of interestedoperators who needs to be allocated channel capacity N S = 5 .Inside the for loop, the algorithm reserves equal portion ofunused opportunistic channel capacity C for the remaining N S interested operators. This is done in line 4 where amaximum channel capacity of CN S is reserved for the κ ( j ) th operator. This step ensures fairness of Waterfilling algorithm.The channel capacity allocated to the κ ( j ) th operator isthe minimum of its modified demand (the required channel Algorithm 1:
Waterfilling Algorithm for OpportunisticChannel Allocation
Input: D O ( t ) , { x k ( t ) } k ∈S Output: { (cid:101) x k,op ( t ) } k ∈S Sort the list { x k ( t ) } k ∈S in ascending order of x k ( t ) .Let κ ( j ) denote the operator index corresponding to the j th position of the sorted list. Set unused opportunistic channel capacity C = D O ( t ) and the remaining number of interested operators toallocate channel capacity N S = |S| . for j ← to |S| do Set (cid:101) x κ ( j ) ,op ( t ) = min (cid:16) x κ ( j ) ( t ) , CN S (cid:17) . Set C = C − (cid:101) x κ ( j ) ,op ( t ) and N S = N S − . Operator Index D e m and Figure 3. Pictorial representation of the example for Waterfilling algorithm. capacity) and the maximum reserved channel capacity of CN S .Accordingly, C and N S are updated in line 5 . In our example,for j = 1 , (cid:101) x ,op ( t ) = min (cid:0) , (cid:1) = 2 and hence the updated C = 17 − and N S = 4 . For j = 2 , (cid:101) x ,op ( t ) =min (cid:0) , (cid:1) = 3 and hence the updated C = 15 − and N S = 3 . For j = 3 , (cid:101) x ,op ( t ) = min (cid:0) , (cid:1) = 4 and hencethe updated C = 12 − and N S = 2 . The loop continuesand finally (cid:101) x ,op ( t ) = (cid:101) x ,op ( t ) = 4 . Waterfilling algorithm prevents wastage of channel capacity by allocating no morethan required channel capacity in line 4 . This ensures that theunused opportunistic channel capacity C is as high as possiblefor the operators with higher customer demand.We end this section by proving that the output (cid:101) x k,op ( t ) ofAlgorithm 1 are iid random variables. By refering to (12) and(11), we can conclude that D O ( t ) and x k ( t ) , which formsthe input to Algorithm 1, are the ouptuts of time invariantfunctions of iid random variables x k ( t ) and V k ( γ ) ( V k ( γ ) decides the random set T ( γ ) in (12)). This implies that D O ( t ) and x k ( t ) are iid random variables as well. Also notethat except the inputs D O ( t ) and x k ( t ) , Algorithm 1 is notdependent on time t . Therefore, Algorithm 1 can be expressedas a time-invariant function of iid random variables D O ( t ) and x k ( t ) . This directly implies that the output (cid:101) x k,op ( t ) ofAlgorithm 1 are iid random variables. Remark 2: Generality of the OSA model.
We want tohighlight that our OSA model is very general for three reasons.
First , any opportunistic channel allocation algorithm can beused as long as (cid:101) x k,op ( t ) are iid random variables. Second ,the parameter φ helps us capture cases where Tier-1 operatorscan/cannot participate in OSA. Third , our model can captureoth overlay and interweave OSA strategy.III. O
PTIMIZATION P ROBLEM
We start this section by formulating the optimization prob-lem for joint spectrum partition as a two- stage StackelbergGame in Section III-A. In the process of formulating theStackelberg Game, we introduce two functions.
First , is therevenue function of an operator which captures the expectedrevenue of an operator in an epoch.
Second , is the objectivefunction which is proportional to spectrum utilization of all theinterested operators in the market. We then develop efficientalgorithms to solve the two stages of the Stackelberg Gamein Section III-B and hence find the optimal M and P whichmaximizes spectrum utilization. In sections III-A and III-B, weassume complete information games. This assumption leads tonotational simplicity. Also, this assumption does not affect thetechnical contribution of the paper as the overall approach canbe easily extended to incomplete information games. This isdiscussed in Section III-C. A. Stackelberg Game Formulation
In this section, we formulate the optimal spectrum partition-ing problem as a two-stage Stackelberg game between the reg-ulator and the wireless operators. The k th operator can be com-pletely characterised by seven parameters which can be rep-resented as a tuple ξ k = (cid:16) µ θk , σ θk , h k ( · ) , σ Rk,a , ρ k , ω k , λ k (cid:17) .In sections III-A and III-B, we assume complete informationgames, i.e. an operator and the regulator knows ξ k of all theoperators. The player in stage-1 of the Stackelberg game is theregulator whose decision variables are M and P . The payoffof the regulator is the expected spectrum utilization over aperiod of Γ ≥ epochs which is given by Q = E (cid:34) Γ (cid:80) γ =1 γT (cid:80) t =( γ − T +1 ( Q lc ( γ, t ) + Q op ( γ, t )) (cid:35) where,(13) Q lc ( γ, t ) = (cid:80) k ∈T ( γ ) (cid:101) x k,lc ( t ) (14) Q op ( γ, t ) = (cid:80) k ∈S (cid:101) x k,op ( t ) (15)In (13), the inner summation is to calculate spectrumutilization over all time slots in an epoch while the outer sum-mation is to calculate spectrum utilization over all the epochs. Q lc ( γ, t ) and Q op ( γ, t ) are the net spectrum utilization intime slot t of epoch γ by using licensed and opportunisticspectrum access respectively. The regulator wants to maximize Q . Using linearity of expectation , we can rewrite (13) as Q = Γ (cid:88) γ =1 γT (cid:88) t =( γ − T +1 E [ Q lc ( γ, t ) + Q op ( γ, t )] (16)We now prove that E [ Q lc ( γ, t ) + Q op ( γ, t )] is not a func-tion of γ and t . Based on (14) and (15), (cid:101) x k,lc ( t ) , (cid:101) x k,op ( t ) ,and T ( γ ) are the only random variables in the expressions of Q lc ( γ, t ) and Q op ( γ, t ) . As discussed in previous sections, (cid:101) x k,lc ( t ) and (cid:101) x k,op ( t ) are iid random variables. T ( γ ) is afunction of bids V k ( γ ) of the operators. Since, V k ( γ ) is an iid random variable, so is T ( γ ) . This discussion implies that Q lc ( γ, t )+ Q op ( γ, t ) itself is an iid random variable and henceits expectation is independent of γ and t . Infact, it is a functionof M , P , S L and S U . Let, U ( M, P, S L , S U ) = E [ Q lc ( γ, t ) + Q op ( γ, t )] (17)Substituting (17) in (16) we get, Q = Γ T U ( M, P, S L , S U ) (18)Equation 18 shows that maximizing Q is same asmaximizing U ( M, P, S L , S U ) . Therefore, we will use U ( M, P, S L , S U ) as the payoff function of the regulator in therest of the paper. U ( M, P, S L , S U ) is also called the objectivefunction as it is a direct measure of spectrum utilization whichwe are trying to maximize in this paper.The players in stage-2 of the Stackelberg game are thecandidate licensed operators, S CL , and candidate unlicensedoperators, S CU . The decision variables of the Stage-2 gameare the set of interested licensed operators, S L , and the setof interested unlicensed operators, S U . The k th operator isinterested in joining the market only if the expected revenuein an epoch is greater than λ k . The expected revenue in anepoch of an interested licensed or unlicensed operator is givenby the revenue function . The formula for revenue functionif different for interested licensed operators and interestedunlicensed operators. The revenue function of an interestedlicensed operator, i.e. k ∈ S L , is R k ( M, P, S L , S U )= E [ R k, ( γ ) | E γk ] P [ E γk ] + E (cid:104) R k, ( γ ) | E γk (cid:105) P (cid:104) E γk (cid:105) (19) = E [ R k,lc ( γ ) | E γk ] P [ E γk ]+ E [ R k,op ( γ ) | E γk ] P [ E γk ]+ E (cid:104) R k,op ( γ ) | E γk (cid:105) P (cid:104) E γk (cid:105) (20) = E [ R k,lc ( γ ) | E γk ] P [ E γk ] + h k (cid:0) µ Xk,op (cid:1) (21)where P [ Z ] denotes the probability of event Z , E γk ( E γk ) is theevent that k ∈ T ( γ ) ( k ∈ T ( γ ) ). In (19), E [ R k,i ( γ ) | E γk ] is the expected revenue of the k th operator if it is a Tier- i operator in epoch γ . Finally, (19) is obtained using the lawof total expectation . Equation 20 is obtained by substituting R k, ( γ ) = R k,lc ( γ ) + R k,op ( γ ) (refer to (8)). Equation 21 isobtained by noticing that the sum of the second and the thirdterm of (20) is equal to E [ R k,op ( γ )] which in turn in equal to h k (cid:16) µ Xk,op (cid:17) according to (7). Similar to the objective function,the revenue function of an interested licensed operator is alsonot a function of epoch γ . This is because the statisticalproperties of the involved random variables R k, ( γ ) and R k, ( γ ) are independent of γ .If the k th operator is an interested unlicensed operator, i.e. k ∈ S U , it is always a Tier-2 operator. Hence, its expectedrevenue in an epoch is R k ( M, P, S L , S U ) = E [ R k, ( γ )] = h k (cid:0) µ Xk,op (cid:1) (22)here E [ R k, ( γ )] = h k (cid:16) µ Xk,op (cid:17) because R k, ( γ ) = R k,op ( γ ) (refer to (9)) and E [ R k,op ( γ )] = h k (cid:16) µ Xk,op (cid:17) .Payoff function of an operator who is interested in joiningthe market either as a licensed or an unlicensed operator is π k ( M, P, S L , S U ) = R k ( M, P, S L , S U ) − λ k (23)where R k ( M, P, S L , S U ) is given by (19) if k ∈ S L and by(22) if k ∈ S U . If an operator does not join the market, itspayoff is zero. An operator decides to enter the market onlyif its payoff π k ( M, P, S L , S U ) is strictly greater than zero.With (23) as the payoff function, Stage-2 game may havemultiple Nash Equilibriums which complicates the analysis.This can be simplified if we assume that the operators are pessimistic in nature. By doing so, we can get an uniquesolution of the Stage-2 game. Pessimistic models to addressthe issue of multiple Nash Equilibriums have been consideredin prior works like [23], [24], [25]. One simple approachto model pessimistic decision making strategy is to use theconcept of dominant strategy , i.e. an operator decides to jointhe market if and only if joining the market is its optimalstrategy irrespective of whether other operators decides to jointhe market. However, in this paper, we model a pessimisticoperators’ decision making strategy using iterated eliminationof strictly dominated strategies (IESDS) [26]. Compared todominant strategy, IESDS is a less pessimistic decision makingstrategy because more operators will join the market.IESDS can be explained as follows. IESDS consists ofiterations. Consider the first iteration which is the originalStage-2 game. We iterate through all the candidate licensedand unlicensed operators to check if either joining the marketor not joining the market is a dominant strategy for any of theoperators. The operators whose dominant strategy is to join(not join) the market will join (not join) the market irrespectiveof other operators’ decisions. This reduces the size of theStage-2 game as it effectively consist of those operators whocould not decide whether to join (not join) the market in thefirst iteration. Such operators are called confused operators in this paper. These confused operators who did not have adominant strategy in the original Stage-2 game may have adominant strategy in the reduced Stage-2 game. Therefore,in the second iteration, we find the dominant strategy ofthe confused operators in the reduced Stage-2 game. Suchiterations continue until convergence, which happens when thereduced Stage-2 game does not have any dominant strategy.It is possible that there are ’confused operators’ even afterconvergence. Those operators will not join the market becausein our model, the operators are pessimistic in nature. B. Solution of the Stackelberg Game
In this subsection, we use the process of backward induction[27] to solve the Stackelberg Game formulated in SectionIII-A. To apply backward induction, we first solve Stage-2of the game followed by Stage-1. The following properties ofthe revenue function (as given by (19) and (22)) are crucialin designing an efficient algorithm to solve Stage-2 of theStackelberg Game.
Property 1. R k ( M, P, S L , S U ) is monotonic decreasingin S L , i.e. R k ( M, P, S L , S U ) ≥ R k ( M, P, S L (cid:83) { a } , S U ) where a / ∈ S L and a ∈ S CL . Property 2. R k ( M, P, S L , S U ) is monotonic decreasingin S U , i.e. R k ( M, P, S L , S U ) ≥ R k ( M, P, S L , S U (cid:83) { a } ) where a / ∈ S U and a ∈ S CU . We have verified these properties numerically using theMonte Carlo integrator which will described in Section IV.These properties can be justified as follows. Property 1 statesthat as the set of interested licensed operators, S L , increases,the revenue function of both the licensed and the unlicensedoperators decreases. The revenue function of a licensed op-erator decreases with increase in S L because the operatorhas to compete with more operators in the spectrum auctionsto get a channel. This reduces the operator’s probability ofwinning spectrum auctions which in turn decreases its revenuefunction as it can effectively serve fewer customer demand.The revenue function of an unlicensed operator also decreaseswith increase in S L . This happens because with increase in S L , there is an increase in the number of operators interestedin opportunistic channel access. This reduces the share ofopportunistic channels for an unlicensed operator. Therefore,its revenue decreases as it can serve fewer customer demand.Property 2 states that as the set of interested unlicensedoperators, S U , increases, the revenue function of both thelicensed and the unlicensed operators decreases. This happensbecause with increase in S U , the share of opportunistic channeldecreases for a licensed or an unlicensed operator. This in turndecreases its revenue function.The psuedocode to solve Stage-2 of the Stackelberg gameis given in Algorithm 2. The inputs of Algorithm 2 areclearly described in Table I. Let S L ( M, P ) and S U ( M, P ) denote the set of interested licensed and unlicensed operatorsif the entire bandwidth is divided into M channels out ofwhich P are licensed channels. S L ( M, P ) and S U ( M, P ) are the outputs of Algorithm 2. As mentioned in SectionIII-A, S L ( M, P ) and S U ( M, P ) are decided by the operatorsbased on IESDS. Algorithm 2 uses Properties 1 and 2 tocompute S L ( M, P ) and S U ( M, P ) in polynomial time whenan operator’s decision making strategy to join/not join themarket is based on IESDS.Let (cid:98) X l and (cid:101) X l , where (cid:98) X l , (cid:101) X l ∈ S CL , denote the set oflicensed operators who are sure to join the market and theset of confused licensed operators respectively till the l th iteration. Note that (cid:98) X l and (cid:101) X l are disjoint sets and the set S CL \ (cid:16) (cid:98) X l (cid:83) (cid:101) X l (cid:17) consists of those licensed operators who aresure not to join the market till the l th iteration. Similarly, (cid:98) Y l and (cid:101) Y l , where (cid:98) Y l , (cid:101) Y l ∈ S CU , denote the set of unlicensedoperators who decided to join the market and the set ofconfused unlicensed operators till the l th iteration respectively.We will now explain the working of Algorithm 2. Algorithm2 starts with iteration . Initially, none of the operators are surewhether to join the market or not; all of them are confused.Hence, in iteration , we initialize (cid:98) X = ∅ , (cid:101) X = S CL , (cid:98) Y l = ∅ and (cid:101) Y = S CU ( line 1 ). The while loop in lines 3-19 finds (cid:98) X l , (cid:101) X l , (cid:98) Y l and (cid:101) Y l for the l th iteration given (cid:98) X l − , (cid:101) X l − , (cid:98) Y l − and lgorithm 2: Algorithm to solve Stage 2 of theStackelberg Game for Joint Spectrum Partitioning
Input: M , P , T , D , α L , α U , S CL , S CU and ξ k ; ∀ k ∈ S CL (cid:83) S CU Output: S L ( M, P ) and S U ( M, P ) Set (cid:98) X = ∅ , (cid:101) X = S CL , (cid:98) Y = ∅ and (cid:101) Y = S CU . Set converged = F alse and l = 0 . while not ( converged ) do Set converged = T rue and l = l + 1 . Set (cid:98) X l = (cid:98) X l − , (cid:101) X l = (cid:101) X l − , (cid:98) Y l = (cid:98) Y l − and (cid:101) Y l = (cid:101) Y l − . for k in (cid:101) X l − do if R k (cid:16) M, P, (cid:98) X l − (cid:83) (cid:101) X l − , (cid:98) Y l − (cid:83) (cid:101) Y l − (cid:17) > λ k then Set (cid:98) X l = (cid:98) X l (cid:83) { k } and (cid:101) X l = (cid:101) X l \ { k } . Set converged = F alse . else if R k (cid:16) M, P, (cid:98) X l − (cid:83) { k } , (cid:98) Y l − (cid:17) ≤ λ k then Set (cid:101) X l = (cid:101) X l \ { k } . Set converged = F alse . for k in (cid:101) Y l − do if R k (cid:16) M, P, (cid:98) X l − (cid:83) (cid:101) X l − , (cid:98) Y l − (cid:83) (cid:101) Y l − (cid:17) > λ k then Set (cid:98) Y l = (cid:98) Y l (cid:83) { k } and (cid:101) Y l = (cid:101) Y l \ { k } . Set converged = F alse . else if R k (cid:16) M, P, (cid:98) X l − , (cid:98) Y l − (cid:83) { k } (cid:17) ≤ λ k then Set (cid:101) Y l = (cid:101) Y l \ { k } . Set converged = F alse . Set S L ( M, P ) = (cid:98) X l and S U ( M, P ) = (cid:98) Y l (cid:101) Y l − of the ( l − th iteration. Since the operators in sets (cid:98) X l − and (cid:98) Y l − will surely join the market, we initialize (cid:98) X l and (cid:98) Y l to (cid:98) X l − and (cid:98) Y l − respectively in the beginning of the l th iteration( line 2 ). The set of confused licensed and unlicensed operators, (cid:101) X l and (cid:101) Y l , are initialized to (cid:101) X l − and (cid:101) Y l − respectively inthe beginning of the l th iteration ( line 2 ). In the for loop in lines 6 - 12 , we check if any licensed operator in set (cid:101) X l − issure to either join or not join the market. Similarly, in the forloop in lines 13 - 19 , we check if any unlicensed operator inset (cid:101) Y l − is sure to either join or not join the market.We will now explain the working of the for loop in lines6-12. The largest possible set of interested licensed operatorsin the (cid:101) l th iteration, for (cid:101) l ≥ l is (cid:98) X l − (cid:83) (cid:101) X l − . This is becausethe operators in set S CL \ (cid:16) (cid:98) X l − (cid:83) (cid:101) X l − (cid:17) are sure not to jointhe market till the ( l − th iteration. Similarly, the largestpossible set of interested unlicensed operators in the (cid:101) l th iteration, for (cid:101) l ≥ l is (cid:98) Y l − (cid:83) (cid:101) Y l − . Therefore, accordingto Properties 1 and 2, the minimum revenue of the k th operator, where k ∈ (cid:101) X l − , in the (cid:101) l th iteration, for (cid:101) l ≥ l is R k (cid:16) M, P, (cid:98) X l − (cid:83) (cid:101) X l − , (cid:98) Y l − (cid:83) (cid:101) Y l − (cid:17) . So if R k (cid:16) M, P, (cid:98) X l − (cid:91) (cid:101) X l − , (cid:98) Y l − (cid:91) (cid:101) Y l − (cid:17) > λ k then joining the market becomes the dominant strategy ofthe k th operator in the l th iteration. Therefore, in line 8 , weremove the k th operator from the set of confused licensedoperators and add it to the set of licensed operators who are sure to join the market. If the k th operator, where k ∈ (cid:101) X l − ,joins the market, then the smallest possible set of interestedlicensed and unlicensed operators in the (cid:101) l th iteration, for (cid:101) l ≥ l are (cid:98) X l − (cid:83) { k } and (cid:98) Y l − respectively. Therefore, accordingto Properties 1 and 2, the maximum revenue of the k th operator, where k ∈ (cid:101) X l − , in the (cid:101) l th iteration, for (cid:101) l ≥ l is R k (cid:16) M, P, (cid:98) X l − (cid:83) { k } , (cid:98) Y l − (cid:17) . So if R k (cid:16) M, P, (cid:98) X l − (cid:91) { k } , (cid:98) Y l − (cid:17) ≤ λ k then not joining the market becomes the dominant strategyof the k th operator in the l th iteration. Therefore, in line 11 ,we remove the k th operator from the set of confused licensedoperators but we do not add it to the set of licensed operatorswho are sure to join the market. The for loop in lines 13-19work in a similar way to decide if any unlicensed operator inset (cid:101) Y l − is sure to either join or not join the market.The variable converged which is declared in line 2 andupdated in lines 9, 12, 16, 19 decides when the while loopterminates. This can be explained as follows. Say that fewof the confused operators in the l th iteration decides to notjoin the market, i.e. if statements in lines 10 or are true .In this case, converged is set to f alse and hence the whileloop continues. Since few of the operators decides not tojoin the market in the l th iteration, then due to Properties1 and 2, the revenue function of the remaining confusedoperators in the ( l + 1) th iteration is more compared to theircorresponding values in the l th iteration. Therefore, it ispossible that for some of these confused operators, joiningthe market becomes the dominant strategy in the ( l + 1) th iteration. The opposite happens when few of the confusedoperators in the l th iteration decides to join the market. Thisdiscussion captures the fundamental idea behind IESDS.Say that after the end of the l tho iteration, (cid:98) X l o = (cid:98) X l o − , (cid:101) X l o = (cid:101) X l o − , (cid:98) Y l o = (cid:98) Y l o − and (cid:101) Y l o = (cid:101) Y l o − . This happenswhen if statements in lines 7, 10, 14, 17 are all f alse . Whenthis happens, converged is true after the end of the l tho iteration and hence the while loop terminates. This is becauseif (cid:98) X l o = (cid:98) X l o − , (cid:101) X l o = (cid:101) X l o − , (cid:98) Y l o = (cid:98) Y l o − and (cid:101) Y l o = (cid:101) Y l o − ,then the value of the revenue function in lines 7, 10, 14,17 in the ( l o + 1) th iteration is same as that in l tho iteration.Therefore, the if statements in lines 7, 10, 14 and 17 will be f alse in the ( l o + 1) th iteration just like the l tho iteration. Thisargument suggests that (cid:98) X l = (cid:98) X l o , (cid:101) X l = (cid:101) X l o , (cid:98) Y l = (cid:98) Y l o and (cid:101) Y l = (cid:101) Y l o for all l ≥ l o and hence convergence in (cid:98) X l , (cid:101) X l , (cid:98) Y l and (cid:101) Y l has been achieved. After convergence is achieved, thereare three kinds of operators. First , the operators in sets (cid:98) X l o and (cid:98) Y l o who are sure that they should join the market. Second , theoperators in sets (cid:98) X l o \ (cid:101) X l o and (cid:98) Y l o \ (cid:101) Y l o who are sure that theyshould not join the market. Third , the ’confused’ operators insets (cid:101) X l o and (cid:101) Y l o . Since our model assumes that the operatorsare pessimistic, confused operators will not join the market.Hence, the set of interested licensed and unlicensed operatorsare (cid:101) X l o and (cid:98) Y l o respectively where l o is the last iteration ofAlgorithm 3 ( line 20 ). Proposition 1.
Time complexity of Algorithm 2 is O (cid:0) N (cid:1) where N = (cid:12)(cid:12) S CL (cid:12)(cid:12) + (cid:12)(cid:12) S CU (cid:12)(cid:12) .roof: The while loop continues until none of the con-fused operators of an iteration have a dominant strategy. Sucha condition is possible at most N times because there are only N candidate operators. Hence, the while loop is executed atmost N times. For a given iteration of the while loop, the innerfor loop in lines 6-12 is executed at most (cid:12)(cid:12) S CL (cid:12)(cid:12) times and thatin lines 13-19 is executed at most (cid:12)(cid:12) S CU (cid:12)(cid:12) times. Therefore, theinner for loops runs at most (cid:12)(cid:12) S CL (cid:12)(cid:12) + (cid:12)(cid:12) S CU (cid:12)(cid:12) = N times. Thisshows that the time complexity of Algorithm 2 is O (cid:0) N (cid:1) .This completes the proof. Remark 3: Efficiency of Algorithm
2. Algorithm 2 usesProperties 1 and 2 to decide whether a confused operator willjoin the market or not by computing its revenue function forthe largest/smallest set of interested operators. Without theseproperties, we have to compute the revenue function for anexponential number of set of interested operators to decidewhether a confused operator will join the market or not.
Remark 4: Comparison with dominant strategy . Only the st iteration of Algorithm 2 is required to find the dominantstrategies of the operators. It is for this reason that the set ofinterested operators, S L ( M, P ) and S U ( M, P ) , will alwaysbe larger if operators’ decision making strategy is based onIESDS rather than dominant strategy.Given that S L ( M, P ) and S U ( M, P ) are the solutions ofthe Stage-2 game, the objective function in (17) can be re-written as (cid:101) U ( M, P ) = U ( M, P, S L ( M, P ) , S U ( M, P )) (24)In Stage-1, the regulator chooses M and P to maximize (cid:101) U ( M, P ) . Let the optimal solution be M ∗ and P ∗ , the optimalvalue of the objective function be U ∗ , where U ∗ = (cid:101) U ( M ∗ , P ∗ ) ,and the optimal set of interested licensed and unlicensedoperators be S ∗ L and S ∗ U , where S ∗ L = S L ( M ∗ , P ∗ ) and S ∗ U = S U ( M ∗ , P ∗ ) . M ∗ and P ∗ are found by performinga grid-search. The grid search is detailed in Algorithm 3. Asshown in lines 2 and 3 of Algorithm 3, the grid search isperformed from M = 1 to a certain M max and from P = 0 to min (cid:0)(cid:12)(cid:12) S CL (cid:12)(cid:12) , M (cid:1) . Note that since spectrum cap is one, thenumber of licensed channels should be lesser than the numberof candidate licensed operators, (cid:12)(cid:12) S CL (cid:12)(cid:12) . Time complexity ofAlgorithm 3 is O (cid:0) M max (cid:12)(cid:12) S CL (cid:12)(cid:12)(cid:1) . C. Incomplete Information Stackelberg Game
In this section, we discuss the generalization to the caseof incomplete information games where the j th operator hasa point estimate (cid:98) ξ jk of the tuple ξ k of the k th operator. Let S C = S CL (cid:83) S CU . We have j ∈ { } (cid:83) S C and k ∈ S C , where j = 0 is the index of the regulator. Also, (cid:98) ξ kk = ξ k becausethe k th operator knows its own tuple ξ k . We now discuss thesteps involved in finding the optimal value of M , P , and theobjective function in this incomplete information setting.The regulator decides the optimal values of M and P . Whilein the complete information case, the regulator exactly knows ξ k ; ∀ k ∈ S C , it now only has an estimate (cid:98) ξ k ; ∀ k ∈ S C in theincomplete information case. Therefore, as far as deciding theoptimal value of M and P is concerned, the regulator simplyuses (cid:98) ξ k ; ∀ k ∈ S C as inputs to Algorithms 2 and 3. Let M ∗ and P ∗ denote the optimal values of M and P . Algorithm 3:
Algorithm to solve Stage 1 of theStackelberg Game for Joint Spectrum Partitioning
Input: T , D , α L , α U , S CL , S CU , and ξ k ; ∀ k ∈ S CL (cid:83) S CU Output: M ∗ , P ∗ , S ∗ L , S ∗ U , and U ∗ Set U ∗ = −∞ . for M ← to M max do for P ← to min (cid:0)(cid:12)(cid:12) S CL (cid:12)(cid:12) , M (cid:1) do Call Algorithm 2 to get the set of interestedlicensed and unlicensed operators, S L ( M, P ) and S U ( M, P ) respectively, for current M and P . Set (cid:101) U = U ( M, P, S L ( M, P ) , S U ( M, P )) . if (cid:101) U > U ∗ then Set M ∗ = M , P ∗ = M , S ∗ L = S L ( M, P ) , S ∗ U = S U ( M, P ) , and U ∗ = (cid:101) U . S L ( M ∗ , P ∗ ) and S U ( M ∗ , P ∗ ) obtained using (cid:98) ξ k ; ∀ k ∈ S C as an input to Algorithm 2 may not be the true set of interestedlicensed and unlicensed operators corresponding to M ∗ and P ∗ . This is because the estimate (cid:98) ξ jk varies between the regu-lator and the operators. Let A jL = S L (cid:18) M ∗ , P ∗ , (cid:110)(cid:98) ξ jk (cid:111) k ∈S C (cid:19) and A jU = S U (cid:18) M ∗ , P ∗ , (cid:110)(cid:98) ξ jk (cid:111) k ∈S C (cid:19) denote the outputs ofAlgorithm 2 with M ∗ , P ∗ , and (cid:98) ξ jk ; ∀ k ∈ S C as its inputs. A jL and A jU are the set of interested licensed and unlicensedoperators corresponding to M ∗ and P ∗ according to the j th operator. The j th operator is interested in joining the marketif and only if j ∈ A jL or j ∈ A jU . Let S trueL ( M ∗ , P ∗ ) and S trueU ( M ∗ , P ∗ ) denote the true set of interested licensedoperators and unlicensed operators respectively. We have, S trueL ( M ∗ , P ∗ ) = (cid:110) j ∈ S CL : j ∈ A jL (cid:111) S trueU ( M ∗ , P ∗ ) = (cid:110) j ∈ S CU : j ∈ A jU (cid:111) To compute S trueL ( M ∗ , P ∗ ) and S trueU ( M ∗ , P ∗ ) , Algo-rithm 2 has to be called (cid:12)(cid:12) S CL (cid:12)(cid:12) + (cid:12)(cid:12) S CU (cid:12)(cid:12) times, one correspondingto each of the (cid:12)(cid:12) S CL (cid:12)(cid:12) + (cid:12)(cid:12) S CU (cid:12)(cid:12) operators. Finally, the true valueof the objective corresponding to M ∗ and P ∗ is U true = U (cid:0) M ∗ , P ∗ , S trueL ( M ∗ , P ∗ ) , S trueU ( M ∗ , P ∗ ) (cid:1) Please note that U true is implicitly a function of the tuples { ξ k } k ∈S C as well. While calculating U true , { ξ k } k ∈S C shouldbe used and not any of its point estimates.IV. M ONTE C ARLO I NTEGRATOR D ESIGN
Algorithms 2 and 3 rely on the computation of the ob-jective function, U ( M, P, S L , S U ) , and the revenue function, R k ( M, P, S L , S U ) . In this section, we design an efficientMonte Carlo integrator to compute these two functions. Thesefunctions are the mean of certain random variables. MonteCarlo integrator estimates mean of a random variable bycalculating the sample mean of the random variable. Considera random variable Z ∼ F Z , where F Z is the probabilitydistribution of Z . Let the mean and the standard deviation of be µ Z and σ Z respectively. The following recursive formulacan be used to compute the sample mean of Z , (cid:98) z r = ( r − (cid:98) z r − + z r r (25)where r is the number of samples, z r is the r th sample of Z and (cid:98) z r is the sample mean of Z calculated over the first r samples. (cid:98) z r is an estimate of µ Z . Note that (cid:98) z r itself is arandom variable with mean µ Z and standard deviation σ Z √ r .According to Chebyshev’s inequality , the probability that (cid:98) z r is within a ∆ bound of µ Z is lower bounded as follows P [ | (cid:98) z r − µ Z | ≤ ∆] ≥ − σ Z r ∆ (26)We want to design a Monte Carlo integrator whose maxi-mum acceptable percentage error in (cid:98) z r is β with a minimumprobability of β . β and β captures the “goodness” ofestimate (cid:98) z r ; a lower β and a higher β implies a betterestimate. To achieve this we substitute ∆ = β µ Z in (26)which makes the RHS of (26) equal to − σ Z rβ µ Z . So wehave to recursively calculate (cid:98) z r until σ Z ≤ rβ µ Z (1 − β ) (27)Inequality 27 can be used as one of the stopping criteriafor the Monte Carlo integrator. However, we don’t know µ Z and σ Z of (27); infact we want to calculate µ Z . One possibleheuristic would be to use the sample mean and the samplevariance in place of µ Z and σ Z respectively. Sample meancan be calculated using (25). Sample variance δz r can becomputed using the following recursive formula [28], δz r = ( r − r δz r − + ( r − (cid:0)(cid:98) z r − (cid:98) z r − (cid:1) (28)To summarize, µ Z is estimated by recursively calculatingthe sample mean using (25) until the following stoppingcriteria is reached, δz r ≤ rβ ( (cid:98) z r ) (1 − β ) (29)Now, we discuss all the sample means which we have tocalculate in order to estimate the objective and the revenuefunctions. By refering to (14), (15) and (17), we can saythat the objective function is the expected value of the netdemand served by all the interested operators in one time slotusing either licensed or opportunistic access. Let (cid:98) U r denotethe sample mean over r samples of the net demand served byall the interested operators in one time slot. Equation 21 showsthat the revenue function of the k th licensed operator consistsof two terms. The first term in (21) is the expected valueof the k th licensed operator’s revenue in an epoch generatedusing licensed access. Let (cid:98) R rk,lc denote the sample mean over r samples of the k th licensed operator’s revenue in an epochgenerated using licensed access. The second term in (21) isthe expected value of the k th licensed operator’s revenue in anepoch generated using opportunistic access. This value is equalto h k (cid:16) µ Xk,op (cid:17) according to (7). But µ Xk,op = (cid:101) µ xk,op T (referto (6)) and hence h k (cid:16) µ Xk,op (cid:17) = h k (cid:16)(cid:101) µ xk,op T (cid:17) . (cid:101) µ xk,op is theexpected value of the demand served by the k th operator in atime slot using opportunistic spectrum access. Let (cid:98) U rk,op be the estimate of (cid:101) µ xk,op over r samples. Finally, (cid:98) R rk,lc + h k (cid:16) (cid:98) U rk,op T (cid:17) is the estimate of the k th licensed operator’s revenue functionover r samples. According to (22), the estimate of the k th unlicensed operator’s revenue function over r samples is h k (cid:16) (cid:98) U rk,op T (cid:17) . To this end, we have to calculate the samplemeans (cid:98) U r , (cid:98) U rk,op and (cid:98) R rk,lc to estimate the objective and therevenue function.We now present a proposition, which helps in generatingrandom samples efficiently for the Monte Carlo integrator. Proposition 2.
Define, ϕ k = DM (cid:90) ϑ f θk ( ϑ ) dϑ + DM ∞ (cid:90) DM ϑf θk ( ϑ ) dϑ − µ θk (cid:101) µ xk,lc (30) where f θk ( ϑ ) is the probability density function of θ k ( t ) .Then, θ k ( t ) , R k,lc ( γ ) and V k ( γ ) are jointly Gaussian randomvariables with joint probability distribution, θ k ( t ) R k,lc ( γ ) V k ( γ ) ∼ N ( ψ k , Σ k ) (31) for all γ and for all t ∈ [( γ − T + 1 , γT ] where, ψ k = (cid:2) µ θk µ Rk,lc µ Rk,lc (cid:3) T (32) Σ k = (cid:0) σ θk (cid:1) ρ k σ Rk,lc σ Xk,lc ϕ k ω k ρ k σ Rk,lc σ Xk,lc ϕ k ρ k σ Rk,lc σ Xk,lc ϕ k (cid:16) σ Rk,lc (cid:17) ω k (cid:16) σ Rk,lc (cid:17) ω k ρ k σ Rk,lc σ Xk,lc ϕ k ω k (cid:16) σ Rk,lc (cid:17) (cid:16) σ Rk,lc (cid:17) (33) Proof:
Please refer to Appendix A for the proof.Recall that in (30), (cid:101) µ xk,lc is given by (3). In (32), µ Rk,lc = h k (cid:16)(cid:101) µ xk,lc T (cid:17) (refer to (6) and (7)). In (33), σ Xk,lc = (cid:101) σ xk,lc √ T where (cid:101) σ xk,lc is given by (4).The pseudocode for the Monte Carlo integrator is givenin Algorithm 4. The sample means (cid:98) U r , (cid:98) U rk,op and (cid:98) R rk,lc areinitialized to zero for r = 0 ( line 1 ). Inside the while loop, (cid:98) U r , (cid:98) U rk,op and (cid:98) R rk,lc are computed recursively until a stoppingcriteria. We discuss the stopping criteria later in this section.In line 5 , the r th sample of θ k ( t ) , R k,lc ( γ ) and V k ( γ ) are generated for all the licensed operator according to theprobability distribution given by (31). We have dropped the γ and t inside the paranthesis for notational simplicity. Similarly,in line 6 , θ k ( t ) is generated for all the unlicensed operator. θ k ( t ) follows the probability distribution N (cid:16) µ θk , (cid:0) σ θk (cid:1) (cid:17) (re-fer to Section II-B). The r th sample of θ k ( t ) , R k,lc ( γ ) and V k ( γ ) are denoted by θ rk , R rk,lc and V rk respectively. Thecustomer demand of the k th operator for the r th sample is x rk = max (0 , θ rk ) . Tier-1 and Tier-2 operators for the r th sample are decided in line 7 . Licensed operators with the P highest bids, V rk , are the Tier-1 operators for the r th sample. T r denotes the set of Tier-1 operators for the r th sample. Theremaining operators, S\T r , are the Tier-2 operators for the r th sample. T r denotes the set of Tier-2 operators for the r th sample. In lines 8-10 , demand served by the operators usinglicensed and opportunistic spectrum access are calculated. lgorithm 4: Monte Carlo Integrator
Input: M , P , T , D , φ , α L , α U , S L , S U , and ξ k ; ∀ k ∈ S L (cid:83) S U Output: U ( M, P, S L , S U ) and R k ( M, P, S L , S U ) ; ∀ k ∈ S L (cid:83) S U Set (cid:98) U = 0 , (cid:98) U k,op = 0 ; ∀ k ∈ S , and (cid:98) R k,lc = 0 ; ∀ k ∈ S L . Set stop = F alse and r = 0 . while not ( stop ) do Set r = r + 1 . For all k in S L , sample θ rk , R rk,lc and V rk fromprobability distribution (31). Set x rk = max (0 , θ rk ) . For all k in S U , sample θ rk from the probabilitydistribution N (cid:16) µ θk , (cid:0) σ θk (cid:1) (cid:17) . Set x rk = max (0 , θ rk ) . Sort the list { V rk } k ∈S L in descending order of V rk . Let T r be the subset of operators in S L with the P highestvalues of V rk . T r are the Tier-1 operators. T r = S\T r are the Tier-2 operators. Demand served by Tier-1 operators using licensedspectrum access are (cid:101) x rk,lc = min (cid:0) x rk , DM (cid:1) ; ∀ k ∈ T r . Calculate modified demand, x rk , and opportunisticchannel capacity, D rO , using (11) and (12) respectively. Call Algorithm 1 to get (cid:110)(cid:101) x rk,op (cid:111) k ∈S , the demandserved by operators using opportunistic spectrum ac-cess. The input to Algorithm 1 are D rO and { x rk } k ∈S . Set (cid:98) U r = ( r − (cid:98) U r − + (cid:80) k ∈T r (cid:101) x rk,lc + (cid:80) k ∈S (cid:101) x rk,op r . Set (cid:98) R rk,lc = ( r − (cid:98) R r − k,lc + R rk,lc r for all k in T r . Set (cid:98) R rk,lc = ( r − (cid:98) R r − k,lc +0 r for all k in S L \T r . Set (cid:98) U rk,op = ( r − (cid:98) U r − k,op + (cid:101) x rk,op r for all k in S . if r = 1 then Set (cid:98) U = 0 , (cid:98) U k,op = 0 ; ∀ k ∈ S , and (cid:98) R k,lc = 0 ; ∀ k ∈ S L . else Set δU r = ( r − δU r − ( r − + r ( (cid:98) U r − (cid:98) U r − ) δU rk,op = ( r − δU r − k,op ( r − + r ( (cid:98) U rk,op − (cid:98) U r − k,op ) ; ∀ k ∈ S δR rk,lc = ( r − δR r − k,lc ( r − + r ( (cid:98) R rk,lc − (cid:98) R r − k,lc ) ; ∀ k ∈ S L Set stop = Stop (cid:16) r, δU r , δU rk,op , δR rk,lc , (cid:98) U r , (cid:98) U rk,op , (cid:98) R rk,lc (cid:17) . Set U ( M, P, S L , S U ) = (cid:98) U r R k ( M, P, S L , S U ) = (cid:98) R rk,lc + h k (cid:16) (cid:98) U rk,op T (cid:17) ; ∀ k ∈ S L R k ( M, P, S L , S U ) = h k (cid:16) (cid:98) U rk,op T (cid:17) ; ∀ k ∈ S U .Demand served by operators using licensed and opportunisticspectrum access for the r th sample are denoted using (cid:101) x rk,lc and (cid:101) x rk,op respectively.The sample means (cid:98) U r , (cid:98) U rk,op and (cid:98) R rk,lc are calculatedin lines 11-14 using recursive formulas analogous to (25).The formula to update (cid:98) U r is shown in line 11 . The term (cid:80) k ∈T r (cid:101) x rk,lc + (cid:80) k ∈S (cid:101) x rk,op is the net demand served by all theoperators in a time slot for the r th sample. The formula toupdate (cid:98) R rk,lc is shown in lines 12 and 13 . If the k th licensedoperator is a Tier-1 operator for the r th sample, then it earnsa revenue of R rk,lc in an epoch using licensed spectrum access( line 12 ). But if the k th licensed operator is a Tier-2 operatorfor the r th sample, then it earns a revenue of using licensedspectrum access ( line 13 ). (cid:98) U rk,op is updated in line 14 . Theoperators serves (cid:101) x rk,op customer demand using opportunisticspectrum access ( line 14 ). The sample variance correspondingto sample means (cid:98) U r , (cid:98) U rk,op and (cid:98) R rk,lc are calculated in lines15-18 . These variances are initialized to zero for the st sample( line 16 ) and updated using recursive formulas similar to (28)for r > ( line 18 ). In line 19 , the Stop ( · ) function decideswhether to stop the Monte Carlo integrator. The stoppingcriteria is based on (29). The Stop ( · ) function returns T rue if and only if all the following conditions are met, δU r ≤ rβ (cid:16) (cid:98) U r (cid:17) (1 − β ) (34) δU rk,op ≤ rβ (cid:16) (cid:98) U rk,op (cid:17) (1 − β ) ; ∀ k ∈ S (35) δR rk,lc ≤ rβ (cid:16) (cid:98) R rk,lc (cid:17) (1 − β ) ; ∀ k ∈ S L (36) r ≥ r min (37)The last condition r ≥ r min ensures that the Monte Carlointegrator samples the mean over atleast r min samples. Wehave used r min = 10000 , β = 1 and β = 0 . unlessstates otherwise. Finally, the estimated values of the objectivefunction and the revenue function are set in line 20 accordingto what we have discussed before in this section (refer to theparagraph before Proposition 2).V. N UMERICAL R ESULTS
In this section, we conduct numerical simulations to bench-mark the algorithms developed in the previous sections. Wealso explore how the optimal solution M ∗ and P ∗ varieswith interference parameters. Throughout this section, eachtime slot has a duration of one week and lease duration oflicensed channels is one year. Hence, T = 52 . In all oursimulations we have: (i) h k (cid:16) µ Xk,a (cid:17) = a k µ Xk,a where a k > . (ii) σ Rk,a = η k h k (cid:16) µ Xk,a (cid:17) where η k > is the coefficientof variation of R k,a ( γ ) . (iii) λ k = Λ k · (cid:0) a k µ θk T (cid:1) where Λ k ∈ [0 , and the term a k µ θk T is the mean revenue ofthe k th operator in an epoch if it can serve all its customerdemand in every time slot. (iv) The maximum capacity ofthe entire bandwidth D is a fraction υ of the sum of µ θk of all the candidate operators, i.e. D = υ (cid:80) k ∈S C µ θk where S C = S CL (cid:83) S CU . Given our choice of h k (cid:16) µ Xk,a (cid:17) , σ Rk,a , and λ k , the tuple ξ k is equivalent to (cid:0) µ θk , σ θk , a k , η k , ρ k , ω k , Λ k (cid:1) in this section.
10 20 30 40 50 60 70 80 90 10000.10.20.30.40.50.60.70.80.91 0 10 20 30 40 50 60 70 8000.10.20.30.40.50.60.70.80.91 (b)(a)
Figure 4. Cumulative distribution function of the percentage increase inobjective function, ∆ U ∗ , for four different types of opportunistic spectrumaccess when the sub-optimal algorithm is: (a) optimizing M while holding P fixed. (b) optimizing P while holding M fixed. A. Benefit of joint optimization of M and P In our first numerical simulation, we analyze the increasein spectrum utilization that one can obtain using joint opti-mization of M and P when compared to optimizing M whileholding P fixed and vice-versa. Our numerical setup is asfollows. There are four candidate licensed operators and nocandidate unlicensed operator. There are parameters whichcompletely defines a market setting: µ θk , σ θk , a k , η k , ρ k , ω k , Λ k , υ , α L , and α U . We generate such market settingsby randomly selecting these parameters from uniformdistributions each of which is associated with a certain range.The range of the parameters µ θk , σ θk , a k , η k , ρ k , ω k , and Λ k for all the operators are [0 . , . , [0 . , . , [0 . , . , [0 . , . , [0 . , . , [0 . , . , and [0 . , . respectively.The range of υ , α L , and α U are [0 . , . , [0 . , . , and [0 . , . respectively. While generating α L and α U , weensure that α L ≤ α U .The optimal value of the objective function correspondingto Algorithm 3 is U ∗ . We compare Algorithm 3 with a sub-optimal algorithm. Let the optimal value of the objectivefunction corresponding to a sub-optimal algorithm be (cid:98) U ∗ .The percentage increase in the objective function is ∆ U ∗ = U ∗ − (cid:98) U ∗ D · . The reason for having D in the denominatoris as follows. The objective function given by (17) is themean demand served by all the operators in one time slotwhich cannot be greater than D , the maximum capacity ofthe entire bandwidth. Hence, U ∗ , (cid:98) U ∗ ≤ D which implies that U ∗ − (cid:98) U ∗ ≤ D . We compute ∆ U ∗ for sub-optimal algorithmsand plot the cumulative distribution function (CDF) of ∆ U ∗ in Figure 4. Recall that φ can be or , and the OSA strategycan be either interweave or overlay. So there are four possiblecombination of OSA. For a given sub-optimal algorithm, wecompute CDFs for all the four combinations.We consider two sub-optimal algorithms. For the first algo-rithm, P is fixed and M is optimized. An intuitive choice of P is the number of candidate licensed operators. In that way,every candidate licensed operators wins a licensed channel inevery epoch. For the second algorithm, M is fixed and P isoptimized. We set M = (cid:4) Dϑ (cid:5) where (cid:98)·(cid:99) is the floor functionand ϑ = |S C | (cid:80) k ∈S C µ θk is the sample mean of the mean of anoperator’s customer demand. This choice of M is to ensurethat the bandwidth ϑ of a licensed channel is neither too high (a) (b) Figure 5. (a) Plots showing the effect of interference parameter α for amarket containing only candidate licensed operators on the optimal numberof channel, M ∗ , and the optimal number of licensed channels, P ∗ . (b) Plotsshowing the effect of interference parameter of licensed channel α L for amarket containing both candidate licensed operators and unlicensed operatorson the ratio of the bandwidth allocated for unlicensed channels, M ∗ − P ∗ M ∗ . that most of it is wasted and neither too low that a licensedoperator has to reject most of its cushtomer demand.In Figure 4, a lower value of CDF for a given ∆ U ∗ implies that the difference in spectrum utilization betweenjoint optimization and the sub-optimal algorithm is higher.By comparing Figures 4.a and 4.b we can say that joint op-timization leads to more improvement in spectrum utilizationwhen P is fixed rather than when M is fixed. Based on Figure4.a, we can say that when P is fixed, joint optimization leadsto more improvement in spectrum utilization for: (i) overlaystrategy than interweave strategy when φ is fixed. (ii) φ = 0 than φ = 1 . Based on Figure 4.b, we can say that when M is fixed, joint optimization leads to more improvementin spectrum utilization for interweave strategy than overlaystrategy when φ is fixed. We don’t observe any such systematictrend for φ when M is fixed. B. Effect of interference parameters
Our second numerical simulation is to study the effect ofinterference parameter on optimal solution. We consider twosimulation setups. The first simulation setup is as follows.For this setup, α L = α U = α . There are candidatelicensed operators and no candidate unlicensed operators.We consider a homogeneous market setting. The minimumrevenue requirement λ k is set to zero for all the operatorswhich ensures that all the operators join the market. Theremaining parameters of the market are: µ θk = 1 , σ θk = 0 . , a k = 1 , η k = 0 . , ρ k = 0 . , and ω k = 0 . for all k ’s.Also, υ = 0 . . We study how M ∗ and P ∗ varies with α .The simulation result is shown in Figure 5.a. Since there areno candidate unlicensed operators, it is intuitive that there areno unlicensed channels, i.e. M ∗ = P ∗ . Figure 5.a shows that M ∗ decreases with increase in α . This can be explained asfollows. If M is low, the bandwidth, and hence the capacityof each licensed channel is high. Therefore, a licensed operatorcan serve more customer demand using the allocated licensedchannel thereby increasing spectrum utilization. But if M istoo low, only few of the licensed operators are allocatedthe licensed channels in an epoch. The remaining operatorswho uses channels opportunistically as Tier-2 operators. Theefficiency of opportunistic access is decided by α . If α islow, it is better to have fewer Tier-2 operators in an epochbecause opportunistic spectrum access is inefficient. This canbe ensured with a higher M so that there are more Tier-1operators in every epoch. Figure 6. Cumulative distribution function of the percentage increase inobjective function, ∆ U ∗ , when the sub-optimal algorithm is to choose thevalue of M and P that maximize the number of interested operators. In our second simulation setup, we include candidate unli-censed operators. The simulation setup is similar to the firstsetup but differs in the following ways.
First , out of the operators, four are candidate licensed operators and four arecandidate unlicensed operators. Second , the interference pa-rameters α L and α U are not same. We set α U = 0 . and vary α L from to . . We study how the ratio of the bandwidthallocated for unlicensed channels characterized by the ratio M ∗ − P ∗ M ∗ changes with α L . This is shown in Figure 5.b. Unlikethe previous simulation setup, the current simulation setuphas candidate unlicensed operators. Therefore, we expect thatthere will be unlicensed channels dedicated for the candidateunlicensed operators. But the question is: what portion ofthe bandwidth should be allocated for unlicensed channels?If α L is high, most of the bandwidth can be reserved forlicensed channels because even if the Tier-1 operators are notusing the licesnsed channels, the Tier-2 operators can use theremaining capacity of the licensed channels efficiently. But as α L decreases, the opportunistic access of licensed channelsbecomes inefficient. Therefore, it is better to reserve higherportion of the bandwidth for unlicensed channels. C. Market competition vs Spectrum Utilization
For most markets, an increase in competition improves thesocial welfare. In our setup, we use the number of interestedoperators, |S L | + |S U | , as the measure of market competitionand spectrum utilization as the measure of social welfare. Inthis numerical simulation, we show that there exists marketsetups where an increase in |S L | + |S U | decreases spectrumutilization. The simulation setup and the definition of ∆ U ∗ issimilar to Section V-A but differs in the following ways. First, in this setup we have three candidate licensed operators andthree candidate unlicensed operators.
Second, the sub-optimalalgorithm in this setup finds M and P that maximize |S L | + |S U | instead of the objective function defined in equation (17).If there are multiple values of M and P that maximize |S L | + |S U | , we choose the ones that maximize the objective functiondefined in equation (17).The simulation result is shown in Figure 6 where weplot the CDF of ∆ U ∗ for market setups. To establishour claim that maximizing |S L | + |S U | doesn’t necessarilymaximize spectrum utilization, we want to find market setupswhere ∆ U ∗ is strictly greater than . We can see that for (1 − . · . of the market setups, ∆ U ∗ > . This establishes our claim that there are market setups, how-ever few, where maximizing |S L | + |S U | doesn’t necessarilymaximize spectrum utilization. However, for these . of themarket setups, ∆ U ∗ is upper bounded by implying onlya marginal improvement in spectrum utilization.VI. C ONCLUSION
In this paper, we designed an optimization algorithm topartition a bandwidth into channels and further decide thenumber of licensed channels in order to maximize spectrumutilization. The access to this bandwidth is governed by atiered spectrum access model inspired by the CBRS band.We first propose a system model which accurately capturesvarious aspects of the tiered spectrum access model. Basedon this model, we formulate our optimization problem as atwo-staged Stackelberg game and then designed algorithms tosolve the Stackelberg game. Finally, we get numerical resultsto benchmark our algorithm and to also study certain optimaltrends of spectrum partitioning and licensing as a function ofinterference parameters.There can be various directions for future research relatedto generalization of the Stackelberg Game model.
First , is tocapture collusion between operators in Stage-2 of the Stack-elberg Game.
Second , in our current model, every operator isassumed to be equally pessimistic. It would be interesting toassociate each operator with a degree of pessimism.A
PPENDIX AP ROOF OF P ROPOSITION t represents the t th time slot of γ th epoch, i.e. t ∈ [( γ − T + 1 , γT ] . We start by derivingthe stochastic model between θ k ( t ) and X k,lc ( γ ) , the netdemand served by the k th operator in γ th epoch using licensedspectrum access. For Gaussian random variables θ k ( t ) and X k,lc ( γ ) , their joint Gaussian distribution is characterised bytheir mean and their covariance matrix. The mean and the vari-ance of θ k ( t ) are µ θk and (cid:0) σ θk (cid:1) respectively. The mean and thevariance of X k,lc ( γ ) are µ Xk,lc and (cid:16) σ Xk,lc (cid:17) respectively. So,to find the joint Gaussian distribution of θ k ( t ) and X k,lc ( γ ) ,all we have to derive is the covariance of θ k ( t ) and X k,lc ( γ ) .To do so, we are going to first recall few notations fromSection II-B. We have, x k ( t ) = max (0 , θ k ( t )) , (cid:101) x k,lc ( t ) =min (cid:0) x k ( t ) , DM (cid:1) and X k,lc ( γ ) = γT (cid:80) t =( γ − T +1 (cid:101) x k,lc ( t ) . Define, X k,lc ( γ, t ) = t − (cid:88) v =( γ − T +1 (cid:101) x k,lc ( v ) + γT (cid:88) v = t +1 (cid:101) x k,lc ( v ) (38)where, (cid:101) x k,lc ( v ) = min (cid:0) max (0 , θ k ( v )) , DM (cid:1) . The mean of X k,lc ( γ, t ) is E (cid:2) X k,lc ( γ, t ) (cid:3) = (cid:101) µ xk,lc ( T − . The covari-ance of θ k ( t ) and X k,lc ( γ ) is ϕ k = E (cid:2)(cid:0) θ k ( t ) − µ θk (cid:1) (cid:0) X k,lc ( γ ) − µ Xk,lc (cid:1)(cid:3) = E (cid:2)(cid:0) θ k ( t ) − µ θk (cid:1) (cid:0) X k,lc ( γ, t ) − (cid:101) µ xk,lc ( T − (cid:1) + (cid:0) θ k ( t ) − µ θk (cid:1) (cid:0)(cid:101) x k,lc ( t ) − (cid:101) µ xk,lc (cid:1)(cid:3) = E (cid:2)(cid:0) θ k ( t ) − µ θk (cid:1) (cid:0) X k,lc ( γ, t ) − (cid:101) µ xk,lc ( T − (cid:1)(cid:3) + E (cid:2)(cid:0) θ k ( t ) − µ θk (cid:1) (cid:0)(cid:101) x k,lc ( t ) − (cid:101) µ xk,lc (cid:1)(cid:3) (39)n (39), θ k ( t ) and X k,lc ( γ, t ) are independent randomvariables. This is because θ k ( t ) are iid random variables andaccording to (38), X k,lc ( γ, t ) is not directly dependent on θ k ( t ) . Therefore, the first term of (39) is zero. We have, ϕ k = E (cid:2)(cid:0) θ k ( t ) − µ θk (cid:1) (cid:0)(cid:101) x k,lc ( t ) − (cid:101) µ xk,lc (cid:1)(cid:3) = E [ θ k ( t ) (cid:101) x k,lc ( t )] − µ θk (cid:101) µ xk,lc = E (cid:20) θ k ( t ) min (cid:18) max (0 , θ k ( t )) , DM (cid:19)(cid:21) − µ θk (cid:101) µ xk,lc = DM (cid:90) ϑ f θk ( ϑ ) dϑ + DM ∞ (cid:90) DM ϑf θk ( ϑ ) dϑ − µ θk (cid:101) µ xk,lc (40)where, f θk ( ϑ ) is the probability density function of θ k ( t ) . Tothis end, we have derived the joint Gaussian distribution of θ k ( t ) and X k,lc ( γ ) . We have, (cid:20) θ k ( t ) X k,lc ( γ ) (cid:21) ∼ N (cid:20) µ θk µ Xk,lc (cid:21) , (cid:0) σ θk (cid:1) ϕ k ϕ k (cid:16) σ Xk,lc (cid:17) (41)The joint Gaussian distribution of θ k ( t ) and X k,lc ( γ ) isgiven by (41), that of R k,lc ( γ ) and X k,lc ( γ ) is given by (7)and that of R k,lc ( γ ) and V k ( γ ) is given by (10). We want toderive the joint Gaussian distribution of θ k ( t ) , R k,lc ( γ ) and V k ( γ ) . Just like θ k ( t ) and X k,lc ( γ ) , all we need to know toderive the joint Gaussian distribution of θ k ( t ) , R k,lc ( γ ) and V k ( γ ) are their mean and covariance matrix. We have alreadydiscussed the mean and variance of θ k ( t ) in the beginning ofthis section. The mean and variance of both R k,lc ( γ ) and V k ( γ ) are µ Rk,lc and (cid:16) σ Rk,lc (cid:17) respectively (refer to (10)).According to (10), the covariance of R k,lc ( γ ) and V k ( γ ) is ω k (cid:16) σ Rk,lc (cid:17) . All we have to derive is the covariance of twoset of random variables. First, θ k ( t ) and R k,lc ( γ ) and second, θ k ( t ) and V k ( γ ) . The following proposition is important forthis derivation. Proposition 3.
Consider jointly Gaussian variables A and B , (cid:20) AB (cid:21) ∼ N (cid:18)(cid:20) µ A µ B (cid:21) , (cid:20) σ A σ AB σ AB σ B (cid:21)(cid:19) (42) where µ A and σ A are the mean and standard deviation of A respectively, µ B and σ B are the mean and standard deviationof B respectively and ρ AB is the correlation coefficient of A and B . Let L ( A | B = b ) denote the conditional distributionof A given B = b . Then, L ( A | B = b ) is a Gaussian distri-bution with mean µ A + σ AB σ B ( b − µ B ) and standard deviation (cid:114) σ A − σ AB σ B . Mathematically, L ( A | B = b ) = N (cid:18) µ A + σ AB σ B ( b − µ B ) , (cid:18) σ A − σ AB σ B (cid:19)(cid:19) (43) Proof:
The proof of a generalized version of this propo-sition can be found in [29, Chapter 3].The covariance of θ k ( t ) and R k,lc ( γ ) is ϕ k = E (cid:2)(cid:0) θ k ( t ) − µ θk (cid:1) (cid:0) R k,lc ( γ ) − µ Rk,lc (cid:1)(cid:3) (44) Let f Xk ( φ ) be the probability density function of X k,lc ( γ ) .Define, h ( φ ) = E (cid:104)(cid:0) θ k ( t ) − µ θk (cid:1) (cid:16) R k,lc ( γ ) − µ Rk,lc (cid:17) | X k,lc ( γ ) = φ (cid:105) (45)where h ( φ ) is the covariance of θ k ( t ) and R k,lc ( γ ) condi-tioned on X k,lc ( γ ) = φ . Using the law of total expectation ,(44) can be equivalently written as ϕ k = ∞ (cid:90) −∞ h ( φ ) f Xk ( φ ) dφ (46)Based on (41) and (7), if X k,lc ( γ ) is given, R k,lc ( γ ) and θ k ( t ) are independent of each other. Therefore, (45) can beequivalently written as h ( φ ) = E (cid:2)(cid:0) θ k ( t ) − µ θk (cid:1) | X k,lc ( γ ) = φ (cid:3) · E (cid:2)(cid:0) R k,lc ( γ ) − µ Rk,lc (cid:1) | X k,lc ( γ ) = φ (cid:3) = (cid:0) E [ θ k ( t ) | X k,lc ( γ ) = φ ] − µ θk (cid:1) · (cid:0) E [ R k,lc ( γ ) | X k,lc ( γ ) = φ ] − µ Rk,lc (cid:1) (47)Using (41) and Proposition 3, E [ θ k ( t ) | X k,lc ( γ ) = φ ] = µ θk + ϕ k (cid:16) σ Xk,lc (cid:17) (cid:0) φ − µ Xk,lc (cid:1) (48)Using (7) and Proposition 3, E [ R k,lc ( γ ) | X k,lc ( γ ) = φ ] = µ Rk,lc + ρ k σ Xk,lc σ Rk,lc (cid:16) σ Xk,lc (cid:17) (cid:0) φ − µ Xk,lc (cid:1) (49)Substituting (47), (48) and (49) in (46) we get, ϕ k = ρ k σ Rk,lc ϕ k (cid:16) σ Xk,lc (cid:17) ∞ (cid:90) −∞ (cid:0) φ − µ Xk,lc (cid:1) f Xk ( φ ) dφ (50)The integral in (50) is the variance of X k,lc ( γ ) which isequal to (cid:16) σ Xk,lc (cid:17) . Therefore, (50) is equal to ϕ k = ρ k σ Rk,lc ϕ k (cid:16) σ Xk,lc (cid:17) (cid:0) σ Xk,lc (cid:1) = ρ k σ Rk,lc σ Xk,lc ϕ k (51)Now, we have to derive the covariance of θ k ( t ) and V k ( γ ) .This derivation is same as the derivation for the covarianceof θ k ( t ) and R k,lc ( γ ) and has been skipped for brevity.We simply state that the covariance of θ k ( t ) and V k ( γ ) is ω k ρ k σ Rk,lc σ Xk,lc ϕ k . Finally, based on our derivation, the mean, ψ k , and the covariance matrix, Σ k , of θ k ( t ) , R k,lc ( γ ) and V k ( γ ) are given by (32) and (33) respectively. Hence, thejoint probability distribution of θ k ( t ) , R k,lc ( γ ) and V k ( γ ) isgiven by (31). This completes the proof.R EFERENCES[1] Federal Communications Commission, “Amendment of thecommission’s rules with regard to commercial operationsin the 3550-3650 mhz band,” 2015. [Online]. Available:https://docs.fcc.gov/public/attachments/FCC-15-47A1.pdf[2] ——, “Amendment of the commission’s rules with regard to commercialoperations in the 3550-3650 mhz band,” 2016. [Online]. Available:https://docs.fcc.gov/public/attachments/FCC-16-55A1.pdf3] ——, “Promoting investment in the 3550-3700 mhz band,”2017. [Online]. Available: https://docs.fcc.gov/public/attachments/FCC-17-134A1.pdf[4] N. Jindal, J. G. Andrews, and S. Weber, “Bandwidth partitioningin decentralized wireless networks,”
IEEE Transactions on WirelessCommunications , vol. 7, no. 12, pp. 5408–5419, 2008.[5] X. Feng, P. Lin, and Q. Zhang, “Flexauc: Serving dynamic demands ina spectrum trading market with flexible auction,”
IEEE Transactions onWireless Communications , vol. 14, no. 2, pp. 821–830, 2014.[6] Q. Liang, H.-W. Lee, and E. Modiano, “Robust design of spectrum-sharing networks,”
IEEE Transactions on Mobile Computing , vol. 18,no. 8, pp. 1924–1937, 2019.[7] X. Lin and J. G. Andrews, “Optimal spectrum partition and modeselection in device-to-device overlaid cellular networks,” in . IEEE, 2013, pp.1837–1842.[8] C. Bazelon, “Licensed or unlicensed: The economic considerations inincremental spectrum allocations,”
IEEE Communications Magazine ,vol. 47, no. 3, pp. 110–116, 2009.[9] T. Nguyen, H. Zhou, R. A. Berry, M. L. Honig, and R. Vohra, “Theimpact of additional unlicensed spectrum on wireless services compe-tition,” in . Citeseer, 2011, pp. 146–155.[10] A. Ghosh and R. Berry, “Competition with three-tier spectrum accessand spectrum monitoring,” in
Proceedings of the Twentieth ACM In-ternational Symposium on Mobile Ad Hoc Networking and Computing ,2019, pp. 241–250.[11] ——, “Entry and investment in cbrs shared spectrum,” in . IEEE, 2020, pp. 1–8.[12] C. Chen, R. A. Berry, M. L. Honig, and V. G. Subramanian, “Compet-itive resource allocation in hetnets: The impact of small-cell spectrumconstraints and investment costs,”
IEEE Transactions on CognitiveCommunications and Networking , vol. 3, no. 3, pp. 478–490, 2017.[13] K. Zhu, E. Hossain, and D. Niyato, “Pricing, spectrum sharing, andservice selection in two-tier small cell networks: A hierarchical dynamicgame approach,”
IEEE Transactions on Mobile Computing , vol. 13,no. 8, pp. 1843–1856, 2013.[14] X. Feng, Q. Zhang, and B. Li, “Head: A hybrid spectrum tradingframework for qos-aware secondary users,” in
Dynamic Spectrum AccessNetworks (DYSPAN), 2014 IEEE International Symposium on . IEEE,2014, pp. 489–497. [15] Y. Chen, L. Duan, J. Huang, and Q. Zhang, “Balancing income and userutility in spectrum allocation,”
IEEE Transactions on Mobile Computing ,vol. 14, no. 12, pp. 2460–2473, 2015.[16] M. R. Hassan, G. C. Karmakar, J. Kamruzzaman, and B. Srinivasan,“Exclusive use spectrum access trading models in cognitive radio net-works: A survey,”
IEEE Communications Surveys & Tutorials , vol. 19,no. 4, pp. 2192–2231, 2017.[17] Y. Chen, P. Lin, and Q. Zhang, “Lotus: Location-aware online truthfuldouble auction for dynamic spectrum access,” in .IEEE, 2014, pp. 510–518.[18] S. Ross,
A First Course in Probability , 8th ed. Pearson, 2009.[19] A. Asheralieva, “Prediction based bandwidth allocation for cognitivelte network,” in . IEEE, 2013, pp. 801–806.[20] A. Sahoo, “Fair resource allocation in the citizens broadband radioservice band,” in . IEEE, 2017, pp. 1–2.[21] D. P. Bertsekas, R. G. Gallager, and P. Humblet,
Data networks .Prentice-Hall International New Jersey, 1992, vol. 2.[22] B. Radunovic and J.-Y. Le Boudec, “A unified framework for max-minand min-max fairness with applications,”
IEEE/ACM Transactions onnetworking , vol. 15, no. 5, pp. 1073–1083, 2007.[23] Y. Tohidi, M. R. Hesamzadeh, and F. Regairaz, “Sequential coordinationof transmission expansion planning with strategic generation invest-ments,”
IEEE Transactions on Power Systems , vol. 32, no. 4, pp. 2521–2534, 2016.[24] I. Gilboa,
Theory of decision under uncertainty . Cambridge universitypress, 2009, vol. 45.[25] E. Hasan and F. D. Galiana, “Electricity markets cleared by merit order -part ii: strategic offers and market power,”
IEEE Transactions on PowerSystems , vol. 23, no. 2, pp. 372–379, 2008.[26] K. Leyton-Brown and Y. Shoham, “Essentials of game theory: Aconcise multidisciplinary introduction,”
Synthesis lectures on artificialintelligence and machine learning , vol. 2, no. 1, pp. 1–88, 2008.[27] A. Mas-Colell, M. D. Whinston, J. R. Green et al. , Microeconomictheory . Oxford university press New York, 1995, vol. 1.[28] B. Welford, “Note on a method for calculating corrected sums of squaresand products,”
Technometrics , vol. 4, no. 3, pp. 419–420, 1962.[29] M. L. Eaton, “Multivariate statistics: a vector space approach.”