Resource Allocation in Heterogeneously-Distributed Joint Radar-Communications under Asynchronous Bayesian Tracking Framework

Optimal allocation of shared resources is key to deliver the promise of jointly operating radar and communications systems. In this paper, unlike prior works which examine synergistic access to resources in colocated joint radar-communications or among identical systems, we investigate this problem for a distributed system comprising heterogeneous radars and multi-tier communications. In particular, we focus on resource allocation in the context of multi-target tracking (MTT) while maintaining stable communication connections. By simultaneously allocating the available power, dwell time and shared bandwidth, we improve the MTT performance under a Bayesian tracking framework and guarantee the communications throughput. Our alternating allocation of heterogeneous resources (ANCHOR) approach solves the resulting nonconvex problem based on the alternating optimization method that monotonically improves the Bayesian Cram\'er-Rao bound. Numerical experiments demonstrate that ANCHOR significant improves the tracking error over two baseline allocations and stability under different target scenarios and radar-communications network distributions.

Abstract-Optimal allocation of shared resources is key to deliver the promise of jointly operating radar and communications systems. In this paper, unlike prior works which examine synergistic access to resources in colocated joint radar-communications or among identical systems, we investigate this problem for a distributed system comprising heterogeneous radars and multitier communications. In particular, we focus on resource allocation in the context of multi-target tracking (MTT) while maintaining stable communications connections. By simultaneously allocating the available power, dwell time and shared bandwidth, we improve the MTT performance under a Bayesian tracking framework and guarantee the communications throughput. Our alternating allocation of heterogeneous resources (ANCHOR) approach solves the resulting non-convex problem based on the alternating optimization method that monotonically improves the Bayesian Cramér-Rao bound. Numerical experiments demonstrate that ANCHOR significantly improves the tracking error over two baseline allocations and stability under different target scenarios and radar-communications network distributions.

I. INTRODUCTION
Spectrum sharing between radar and communications systems has lately garnered significant research interest in order to efficiently exploit the electromagnetic spectrum, a scarce natural resource [1,2]. While radars require wider bandwidths to deliver a high-resolution sensing performance [3], mobile communications are witnessing an ever-increasing demand for broader bandwidths in order to provide a designated qualityof-service (QoS) [4]. Spectrum-sharing solutions address these competing requirements either by opportunistic spectral access of broadcast systems by passive radars [5]; spectrally coexisting legacy systems that operate by mitigating the mutual interference with minimum design modifications at the cost of performance [6][7][8]; or spectrally co-designed joint radarcommunications (JRC) systems [9][10][11] that utilize a common transmit/receive hardware and waveform.
Their work is supported in part by ERC AGNOSTIC under grant EC/H2020/ERC2016ADG/742648 and in part by FNR CORE SPRINGER under grant C18/IS/12734677.
The conference precursor of this work was presented in the 2021 IEEE International Workshop on Signal Processing Advances in Wireless Communications (SPAWC).
A successful operation of these spectrum sharing systems relies on an efficient radio resource optimization [12][13][14]. In coexistence paradigms with single antenna or with a fixed beamforming, the resources may include total transmit power, transmit signal bandwidth, and transmission time slots [9]. Additionally, in phased array or multiple-input multiple-output (MIMO) systems, antennas need to be reserved for optimal sharing [15]. In single antenna JRC systems, the resource allocation objective is to optimize the transmit energy of the dual-purpose waveform based on the propagation channels of radar and communications users. In JRC systems employing a transmit antenna array, highly directional beamforming toward the radar surveillance area while also ensuring sidelobes on the users offers an interesting solution [16]. In JRC systems involving hybrid analog-digital transceiver architectures [17], antenna selection may also be viewed as a component of resource allocation objective.
While colocated JRC solutions continue to be comprehensively investigated, distributed JRC systems involving widely separated radars and communications transmitters remain relatively unexamined. Among prior works, the spectrum-sharing model in [18] comprises a widely distributed MIMO radar [19] but studies coexistence with a simplistic point-to-point MIMO communications. The distributed radar proposed in [20] exploits a communications waveform but operates in passive (receive-only) mode. The co-design framework in [21] develops new waveforms and precoders for a distributed but identical units of radar and communications. The displaced sensor imaging in [22] considers identical units of radar sensors on the same vehicular platform that coordinate timing via communications protocols.
There are very few resource allocation studies for identical or homogeneous distributed JRC systems [23]. In such a system, the objective is to optimize the transmission strategies (e.g., transmit power, duration, bandwidth) for identical transmit systems in order to minimize the target localization error in a dynamic scenario and maximize the communications capacity. In the context of distributed radar resource allocation for target localization and tracking, seminal works have appeared recently [24][25][26][27][28]. In [24], a joint beam and power allocation scheme was developed for a colocated MIMO radar network. In [25], target assignment and dwell time allocation were considered simultaneously in a phased array radar network. Subsequently, in [26], an adaptive power allocation was proposed with dynamic tracking threshold adjustment based arXiv:2108.10432v2 [eess.SP] 4 Mar 2022 on the asynchronous working manner of the radar networks. This was extended in [27] by considering heterogeneous radar (HetRad) network, i.e. comprising different types of radar units, for multi-target tracking (MTT).
However, aforementioned prior works do not consider spectrum sharing paradigm with wireless communications. Additionally, the distributed nature of radar deployment, strong likelihood of encountering heterogeneous radar units in practice, and an ever-increasing deployment of cellular communications potentially makes optimal resource allocation difficult in a dynamic scenario. In each of these cases, the resulting spectrum reuse leads to interference between the two systems and impacts the resource allocation strategy. In [23], resource allocation is considered for target localization excluding additional complexities of a heterogeneous radar network and multi-tier communications. Some recent state-ofthe-art works [29,30] on combining sensing functionalities in heterogeneous wireless networks (HetNets) explore the effect of varied communications coverage areas on the joint system performance. However, these studies do not explore resource allocation and do not, per se, qualify as JRC systems.
In this paper, we consider a heterogeneous mix of radar units and multi-tier communications to propose optimal allocation of resources such as power, dwell time, and bandwidth. Our approach optimizes the radar tracking performance while maintaining a reliable QoS in communications systems. Preliminary results of this work appeared in our conference publication [31] where multi-tier communications and various systems configurations were ignored. In this work, we further generalize the study in [31] to the following aspects: First, from system model perspective, we include target Doppler in the radar measurement model and estimate all key parameters, − range, angle-of-arrival (AOA), and Doppler velocity − of a target. Next, we model the communications system as a multi-tier cellular network. We also consider the allocation of the shared frequency bands besides power and dwell time in the joint HetRad-HetNet scenario. Finally, a novel algorithm based on the alternating optimization framework is developed to solve the resulting non-convex problem. For each subproblem, we further propose an iterative method to solve with low computational complexity and proven convergence analyses. Our proposed alternating allocation of heterogeneous resources (ANCHOR) algorithm for distributed resource allocation shows roughly 10 times improvement of tracking performance in the sense of root mean square error over the baseline allocations.
The rest of this paper is organized as follows. In the next section, we introduce the configurations and system model of the heterogeneous JRC network. We formulate the resource allocation problem in Section III and describe the proposed ANCHOR algorithm in Section IV. We validate our models and methods via numerical experiments in Section V. We conclude in Section VI.
Throughout the paper, we reserve boldface lowercase and boldface uppercase for vectors and matrices, respectively. We denote the transpose, conjugate, and Hermitian by (·) T , (·) * , and (·) H , respectively. The notation I N is the N × N identity matrix, and 1 N is a vector with all N elements being 1. x m i,q,k (ẏ m i,q,k ) Velocity along x (y) axis of target q at time t m i,q,k The Kronecker product is denoted by ⊗; || · || p is the p norm; and · 0 is the number of non-zero elements of the vector. The notation Tr {·} is the trace of the matrix, | · | is the cardinality of a set, E [·] is the statistical expectation function, ∇ is the gradient operator, and N (·, ·) represents the probability density function of the normal distribution. Table I summarizes some of the important notations in this paper.

II. NETWORK CONFIGURATION AND SYSTEM MODEL
Consider a heterogeneously-distributed radar and communications network (HRCN) (Figure 1), which employs different types of radars that simultaneously track multiple targets. The radars operate in the same spectrum as a multi-tier communications system that comprises a macro base station (BS) and multiple micro and macro users. The HRCN is equipped with a fusion center for signalling, synchronization, and resource allocation. Here, an efficient use of spectrum is possible provided the resources are well allocated amongst radar and communications entities. The HRCN is embedded with the following heterogeneities:  Heterogeneity in network architecture The HRCN has both radar and cellular network components. The radar network comprises MIMO, phased array, and mechanically scanning radars that are heterogeneously distributed while tracking the same target. The multi-tier cellular communications network has macro and micro BSs along with multiple macro users. Heterogeneity in allocated resources Both tracking and communications QoS are affected by respective signalto-interference-plus-noise ratios (SINRs), which further depend on the transmit power, radar's dwell time, and bandwidth usage. To mitigate the mutual interference and enhance the overall performance, these heterogeneous resources must be judiciously allocated. In the following, we describe each of these HRCN configurations in detail.

A. Configurations of the HRCN
The HRCN is a complex system with several operating units of radars and communications that require numerous system parameters to be optimized for multiple resources.
a) Radar set-up: Consider N radars in the HRCN to track Q widely separated point-like targets in the surrounding environment. In particular, assume N c colocated MIMO radars, N p phased array radars, and N m mechanical scanning radars with N = N c + N p + N m , located at different places with the coordinates {(x i , y i )} N i=1 . Define the sets ϕ c {1, . . . , N c }, ϕ p {N c + 1, . . . , N c + N p }, and ϕ m {N c + N p + 1, . . . , N } to index these radars. The operational scheme of each of the radars are as follows: Colocated MIMO radar For i ∈ ϕ c , the radar i is capable of pointing multiple beams to illuminate multiple targets simultaneously. Herein, the multiple targets are illuminated for the same dwell time but with different transmit powers. Hence, the revisit time intervals for the Q targets are the same. Phased array radar For i ∈ ϕ p , the radar i steers the beam to illuminate multiple targets sequentially by configuring the phases of the transmit array. In this scheme, each target is illuminated with identical transmit power but with a different dwell time. Hence, the revisit time intervals for the Q targets are usually different. Mechanical scanning radar For radars indexed by ϕ m , the system parameters are usually fixed. The radar scans mechanically and illuminates all the targets sequentially with identical power and dwell time (and hence identical revisit times).
The scanning operations of each radar are depicted in Figure 2, where M i,q,k is the number of measurements collected from target q by radar i during the k-th fusion interval 1 (from t k to t k+1 ); t m i,q,k is the arrival time of the m-th measurement for m = 1, . . . , M i,q,k ; P m i,q,k and T m i,q,k denote, respectively, the transmit power and dwell time corresponding to the measurement at time t m i,q,k . Assume that the fusion interval is a constant denoted by T 0 . The inter-measurement duration ( Fig. 2a) at ith MIMO radar for the qth target (revisit intervals) is denoted by a constant T i,q . However, P m i,q,k varies for q = 1, . . . Q. In case of phased-array radars (Fig. 2b), T i,q are different but P m i,q,k is identical for q = 1, . . . Q. Finally, for a mechanically scanning radar (Fig. 2c), both T i,q and P m i,q,k are the same for q = 1, . . . Q.
b) Wireless communications configuration: Consider a multi-tier heterogeneous network (HetNet) of communications ( Fig. 1) that consists of a macrocell (primary network) with multiple macro users and a microcell (secondary network). The primary network of the macro BS serves J macro users and one micro BS. The secondary network is served by the micro BS with L micro users. The co-existence of the multitier networks, frequency reuse, and precoding has been wellstudied in literature [34]. For the HRCN resource allocation, we make the following model to simplify the problem. Communications downlinks The coverage areas of the macro and micro BSs are different. When reusing the frequencies, we assume that the macro downlinks cause interfere to only L micro users. Similarly, the micro downlinks interfere within the coverage area of only micro BS. This is especially true when the micro network is power limited. SINR of macro user The macro BS employs orthogonal frequency-division multiple access (OFDMA) for downlink and allocates non-interfering resource blocks. As a result, the mutual interference among the J macro downlinks is non-existent. The SINR of the j-th macro user is 2 1 Fusion interval [32,33] is the period during which multiple sensors measure and estimate independently. The target state is determined by fusing all radar estimates. 2 The radar illumination is irregular and non-continuous in a fusion interval. Hence, the commonly used transient SINR is unavailable. Instead, we use the average SINR as the communications metric. Additionally, the averaging period may be replaced by the airtime Tc of the communications transmission, which has the same expression but with a scaling constant Tc/T 0 .  where β j and P j c are the channel gain and allocated power to the j-th macro entity, respectively; α r j,i is the interference channel gain from the i-th radar; and σ 2 c,j is the noise power at the receiver front-end of jth user. SINR of micro user To improve spectrum utilization, the micro BS may reuse the same band of the macro BS [35] and then apply the (non-interfering) OFDMA resource blocks amongst the L micro users. This leads to negligible mutual interference among the L micro users. Compared to the SINR of macro entities, the SINR for the l-th micro user becomes where δ j P j c T 0 represents the interference from the j-th macro user of the same frequency; δ j is the interference channel gain; and the other parameters are similar to their counterparts in (1). Codebook We assume Gaussian codebook for all transmissions. c) Configurations for shared bandwidth: Assume B to be the total bandwidth available to the HRCN. The sub-carrier interval of OFDMA ∆f and the spectrum is divided into F sub-channels with B = ∆f (F − 1). The allocated bandwidth is configured as follows: Bandwidth configuration for communications To serve the J macro entities, the communications network requires where F j c is the number of sub-channels for the j-th user. A selection vector f c j,k ∈ R F , comprising of binary entries, represents the selected sub-channels for the j-th downlink in the k-th fusion interval. Similar operations are ascribed to the L micro users. The bandwidth allocation (i.e. f c j,k J j=1 ) may be reassigned for each fusion interval to meet the changing demands of the wireless communications. Bandwidth configuration for radars The bandwidth required by the i-th radar is B i = ∆f · F i r . Similar to communications, a selection vector f r i,k ∈ R F represents the selected sub-channels for the i-th radar during the kth fusion interval. In this work, we fix the allocation  to N radar bands (i.e. the occupied sub-channels) implying a higher priority to radar performance than communications. This also simplifies the problem formulation because only communications bands need to be assigned. However, this may result in limited exploitation of the opportunities offered by a flexible HRCN. Fig. 3 illustrates the bandwidth shared by radars and communications, where some communications sub-channels partially/completely overlap with those used by radars. d) HRCN interference: During a fusion interval, multiple radar measurements coexist with J communications downlinks. Consequently, we have the following cases of interference within the HRCN: Intra-system interference As mentioned earlier, there is no intra-tier interference arising from the non-interfering OFDMA resource-block allocation. But spectrum reuse results in cross-tier interference to the micro users. For N radars, we do not impose any limitation on the bandwidth overlap. The interference among the N radars is assumed to be negligible because all the constituent radars adopt directional beams. Inter-system interference This refers to the mutual interference between radars and communications. The coefficients accounting for the interference from the jth macro downlink to the i-th radar and from the ith radar to the j-th downlink are denoted by α c i,j and α r j,i , respectively. Both α c i,j and α r j,i are determined by the spectral overlap. The in-band interference from communications to radar is proportional to the overlap [36]. Thus, i is a constant unrelated to the bandwidth. The mutual interference between the i-th radar and the l-th micro downlink is analogously defined.

B. Target tracking model
During the k-th fusion interval, the state of tar- denote target's position and Doppler velocity, respectively, in Cartesian coordinate system. The state T . The state transition and measurement model is where s m i,q,k represents the state of target q to radar i at the represents the process noise, w m i,q,k ∼ N 0, Σ m i,q,k represents the measurement error, and is the measurement function given by (4) with R m i,q,k , θ m i,q,k and ν m i,q,k being the range, AoA and velocity, respectively.
and σ 2 ν m i,q,k are the lower bounds on the MSE error of the corresponding parameters in the subscripts. As per [38,39], where T m i,q,k is the dwell time 4 of radar i on target q, α c i,j accounts for the interference coefficient from the j-th 3 If the channel is frequency-selective [37], then α c where the vectorα c i,j represents frequencyselectivity. 4 The dwell time is defined as the product of number of pulses for each measurement and pulsewidth [39]. communications downlink to the i-th radar, σ 2 r,i is the variance of the noise at the receiver of radar i, η m i,q,k is the radar crosssection of the target q to radar i at the receive time t m i,q,k , ζ i and B i are the transmit signal bandwidth and the 3dB receive beam width of radar i, respectively, and c R , c θ and c ν are unrelated constants. Therefore, we rewrite where C m i,q,k is the matrix that contains the above-mentioned parameters.

III. PROBLEM FORMULATION
The HRCN configurations presented above indicate that the constituent radars operate in an asynchronous manner. To circumvent the demanding synchronization requirement, we adopt the concept of composition measure (CM), i.e., the CM estimate will serve as the input measurements to the tracking filter [40]. Thus, the estimation accuracy of the CM will determine the tracking performance. It can be gauged from (3)-(6) that the Cramér-Rao bound (CRB) of CM would depend on resources allocated to different entities. Further, the overall filter estimates also depend on these allocations. Here, contrary to [27], our HRCN performance involves the additional dimension of communications which influences resource allocation and impacts CRB derivations.

A. Composition measure and Bayesian CRB
A composition measure (CM)ŝ q t k+1 for each target q at the fusion time t k+1 is obtained based on all the measurements during the k-th fusion interval y q k . The CMŝ q t k+1 is an estimate of the true state s q t k+1 and shown to be statistically efficient even for small sampling sizes [41]. In particular, The measurements are independent and the probability density function where s m This requires solving a nonlinear least-square problem and the iterated least squares algorithm [42] is usually applied to obtain the ML estimate. Following [27], we approximate the CM covariance matrix by the CRB of s q The CMŝ q t k+1 and its approximate covariance matrix are the input to the tracking filter. We denote the filtered estimate bys q t k+1 , which is function ofŝ q t k+1 . The mean-squared error (MSE) of s q t k+1 satisfies the inequality where B s q t k+1 −1 is the Bayesian CRB matrix, and B s q t k+1 is the corresponding FIM that is adapted from [43] by further incorporating the interference from communications system. These steps produce, whereĤ m i,q,k is the Jacobian of the measurement ats q t k+1|k that is the predicted estimate of the true target state.

B. Resource allocation
To avoid frequent reconfigurations of the system, the allocated power and dwell time from radar i on for target q are constant over multiple measurements of a fusion interval, i.e., Recall that the subchannel selection vectors of j-th user and i-th radar are, respectively, where f r i,k is predefined and fixed over fusion intervals. The vectors f c j,k satisfy the following constraints, Based on the above definitions, the bandwidth overlap of the i-th radar and the j-th user is simply f r The FIM of the CMs becomes Define the following metric of estimation accuracy where Tr ΛB −1 s q t k+1 Λ T accounts for the total MSE error of s q t k+1 . Therefore, maximizing this metric will lead to reduction in the MSE for each target estimate.
It is natural to consider both macro and micro throughput in resource allocation for HRCN system. However, this results in a centralized solution to the multi-tier cellular network, which is computationally expensive because of the high crosstier signaling overhead [35]. Therefore, assuming that the microcell resource is allocated independently within the micro BS, we consider the joint allocation over the radars and the macro BS of the multi-tier cellular network. Without loss of generality, our proposed algorithm is also applicable when the micro user throughput is included with appropriate constraints because of the linearity of interference terms in the denominator of (2). This leads to a two-stage optimization procedure, which could be one-shot or used in a loop; we omit this here and focus on only macroscopic problem.
Define the optimization variables z k = P 1,1,k , . . . , P Nc,Q,k , T Nc+1,1,k , . . . , T Nc+Np,Q,k , The resource allocation requires solving the optimization In (26), the first constraint represents the throughput requirements on the J macro downlinks. The second and third constraints, respectively, account for radar's power and dwell time budgets. The last two constraints represent the budgets of communications power and bandwidth, respectively. This is a general formulation that subsumes the problem in [27] as a special case when only radar's power and dwell time are considered. For the sake of simplicity, we considered only the macro users. However, the problem can be easily extended to micro users and, importantly, our proposed ANCHOR algorithm is still applicable. Overall, the problem is highly challenging because the objective function as well as some constraints are non-convex; variables z k and F c k are coupled; and constraints on F c k imply a mixed-integer programming.

IV. OPTIMIZATION METHOD FOR RESOURCE ALLOCATION
We now develop an algorithm to solve (26) of the k-th fusion interval and, for notational simplicity, we omit the subscript k hereafter. We employ the alternating optimization to decouple the optimization for z and F c ; this corresponds to an alternating update of the binary frequency allocation and the continuous power and dwell time allocation.

A. Frequency allocation
For a fixed z at the -th iteration, the subproblem of optimizing F c in problem (26) is Considering the special structure of the constraint set F c specified in (19) (clearly, N f is an integer). Thus, problem (27) is equivalently recast into where B q in (28) is recast as Denote the optimal solution of (27) and (29) by f c j = s j 1 F c min and s j , respectively. The problem in (29) can be interpreted as an assignment problem with a nonlinear objective function. The constraints in (29) indicate that the problem is to assign an exclusive location to the only non-zero element, i.e. 1, for each s j so that the objective value is maximized. For the classical assignment problem, Hungarian method [44] is a classical approach that yields the global optimal solution with a cubic complexity. However, it can only solve for linear objective function. Although many subsequent works extend it to the nonlinear cases, most of them are heuristic and only available to bilinear, quadratic, biquadratic, and cubic objectives [45]. The objective function of problem (29) does not belong to any of the above-mentioned classes. Consequently, prior approaches to the nonlinear assignment problems are inapplicable here.
Since the number of feasible assignments is finite, an exhaustive search over the whole feasible space can always yield the optimal solution in polynomial time. However, for large N f , this approach suffers from heavy computations. Here, simulated annealing method offers a trade-off between computational load and optimality while solving (29). In traditional simulated annealing, computational time and complexity are dictated by the exploration/exploitation mechanism and the termination criteria potentially leading to global optimum. However, in the current setting, fewer and faster computations are preferred and hence, in the sequel, we adapt the classical simulated annealing to account for these constraints. a) State Transition: Note that in problem (29), each s j should satisfy the constraints The correspondingŝ j ∈ R N f ×1 , ∀j = 1, . . . , J is where [·] n represents the n-th element of a vector. Thus, a feasible s j is generated by selecting one nonzero position from s j . Through this generation approach, the original searching space is reduced. Accordingly, we propose state transform (i.e. moving from the current {s j } to a neighbor) in Algorithm 1.
This problem-specific state transition enables faster neighbourhood discovery within the feasible set. Randomly selectn ∈ N as the nonzero location of sĵ 5: N ← N \n 6: if s T j s j = 0, ∀j ∈ J \ĵ then 7: Exit 8: end if 9: end while 10: Definej be the index associated with s T j sj = 0 11:ĵ ←j and go to Line 3 b) Optimality requirement: Denote the state transition in Algorithm 1 by H (·) and combine it with the simulated annealing algorithm summarized in Algorithm 2. Note that the simulated annealing , as a global optimization approach, has the theoretical guarantee of convergence in the probabilistic sense [46]. However, the convergence cannot be guaranteed for a specific run/instance. Algorithm 2 is a nested component of the alternating optimization framework. So, it is viable to terminate the iterations of the algorithm as long as the objective function of problem (27) improves over the previous iteration. This mitigates the demanding requirement on the simulated-annealing-based algorithm to find a global optimum while offering the designer a flexibility in the selection of the termination criterion.

Algorithm 2 Global optimal search for frequency allocation
if ∆ > 0 then To allocate power and time-slots, we adopt the alternating ascent-descent method.
1) Maximin reformulation: For a fixed F c +1 , the subproblem of optimizing z is The following Lemma 2 casts (32) to a more convenient form.
Proof: See Appendix B. Clearly, through this reformulation, we get rid of the inverse operation on B q .
2) Alternating descent-ascent method: For problem (33), we resort to the alternating ascent-decent method, which has been widely applied to solve maximin problems 5 and has several variants or extensions [47,48]. The convergence of these methods can be guaranteed under some mild conditions. For more details, we refer the interested reader to [49,50].
At the n-th iteration of the alternating ascent-descent method, given a fixed z n , the inner minimization problem of {V q } has the closed form solution where B ,n q is obtained by substituting F c +1 and z n into (21).
where ω n i,q = Tr Λ V n+1 q TC i,q Λ V n+1 q ≥ 0, and the constant term is unrelated to z. Thus, the maximization problem of z becomes Rewrite the objective function in (36) as The constraints on z of problem (36) expressed compactly are and Therefore, problem (36) is equivalently recast into which maximizes a sum of linear fractional function over a convex set.
Several methods have been proposed in the literature to solve this fractional programming problem [51][52][53][54]. However, in order to be within the framework of the alternating ascentdescent method, an incremental update on the variable z is necessary. In other words, a global solution to problem (48) might violate the alternating ascent-descent method and its theoretical guarantees. Therefore, within the alternating ascent-descent framework, the update rule of z is simply z n+1 = P Z (z n + η∇ z f (z n )), where η is the stepsize parameter, and P Z (·) represents the projection to the convex set Z = {z|Az ≤ b, z ≥ 0}, which actually involves solving a convex quadratic problem.

3) Computational complexity and convergence analysis:
The proposed method for the power and time allocation is summarized in Algorithm 3. The overall computational complexity of Algorithm 3 is linear with the number of iterations.
Algorithm 3 Alternating ascent-descent method to power and dwell time allocation  (38), (39), (41), (42) and (43) 8: is also linear with O (QN ). The update rule of z n+1 involves the projection operator P Z (·), which is, in fact, to solve a convex quadratic problem. The MOSEK optimization package will reformulate the problem into epigraph form by introducing an additional slack variable. Consequently, there is a second order cone constraint and the remaining constraints are linear. The computational complexity of solving the reformulated problem is therefore upper bound by O (N c Q + N p Q + J) 3.5 , the same order as second-order cone programming. The convergence of Algorithm 3 is stated by the following theorem.
Lemma 3. Assume {z n } to be the generated sequence of Algorithm 3. Then, every limit point of the sequence {z n }is the stationary point of problem (32).
The convergence of the combined algorithm ANCHOR to solve problem (26) is stated in the following Lemma 4.   update F c +1 via Algorithm 2 5: update z +1 via Algorithm 3 6: ← + 1 7: until convergence

C. Algorithmic workflow
The ANCHOR algorithm is applied to allocate resources for a fusion interval. It is expected that the optimized resource allocation will improve the accuracy of the CMŝ q t k+1 and its covariance matrix J −1 y q k s q t k+1 , which further enhances the tracking performance based on Kalman filter. At the k-th fusion interval,s q t k|k is the filtered estimate of the target state s q t k of the last fusion interval and assumed to be known, and s q t k+1 is the CM estimated based on the measurements y q k . Based on the system model as equation (3), we have wheres q t k+1|k ands q t k+1|k+1 are the predicted and filtered estimate of the true target state, respectively, and C q k+1|k and C q k+1|k+1 are the corresponding covariance matrices, and K q k+1 is Kalman gain with The filtered estimates q t k+1|k+1 is the prior knowledge of the Kalman filter for the next fusion interval. Note that the initial states q t 1|1 and its covariance matrix C q 1|1 are usually estimated. The workflow of the resource allocation for the HRCN is summarized in Fig. 4. We observe that the entire resource allocation procedure proceeds in an iterative closed-loop as the Kalman filtering process. As the fusion interval elapses, the resource allocation iteratively improves the Bayesian CRB matrix of Kalman filtering thereby leading to an enhanced tracking performance in the sense of smaller estimation covariance.

V. NUMERICAL EXPERIMENTS
We conducted extensive numerical experiments to evaluate the performance of our proposed resource allocation algorithm ANCHOR for the HRCN scenario. Throughout the experiments, unless specified otherwise, we employ the following parameter settings: 1) HRCN scenario settings: We consider a typical outdoor scenario (Fig. 5), where target 1 and 2 start from the locations (−2 km, −4 km) and (4 km, 2 km) and move with the velocities (50 m/s, 50 m/s) and (−25 m/s, −50 m/s), respectively. We deploy heterogeneously-distributed radars to track the two moving targets. Specifically, we set N cr = 3 colocated MIMO radars (MMRs), N pr = 3 phased array radars (PARs) and N mr = 2 mechanical scanning radars (MSRs), and their locations are randomly generated in the region of interest. The duration of a single fusion interval is T 0 = 10 s. With respect to the staring time t = 0 of a fusion interval, the initial and revisit time of these radars are provided in Table II. For the HetNet, we consider only the downlink transmission from a BS to J = 6 devices, which are randomly placed within the BS coverage area. For the shared bandwidth between radars and communications, assume that the total bandwidth is B = 400 MHz available, for example, in X-band and higher spectral bands) and the unit interval is f = 4 MHz (i.e. B/ f = 100 frequencies to be assigned). Fig. 6 demonstrates pre-assigned radar frequencies, where each block represents a 40 MHz frequency band. 2) Algorithmic parameter settings: We set T max = 1000, T min = 0.1 and δT = 1. For the alternating ascentdescent method to power and dwell time allocation, both power and dwell time are initialized using the uniform allocation (i.e. P i,q,k = P i total /N c , P j c,k = P c total /J and T i,q,k = T i total /N p ), and the stopping criterion is the increase of the objective value is less than 1 × 10 −4 . The ANCHOR algorithm stops when the objective value increase is less than 1 × 10 −4 or the number of iterations is beyond 100.
3) System parameter settings: All interference channel gains 6 and noise power are obtained by squaring the value generated from the Gaussian distribution. In the Kalman filter, justified by the fact that the tracking is along a series of fusion intervals, the initial states q t 1|1 is set as the true target state with Gaussian noise of standard deviation 5% of the true value. The corresponding initial covariance 6 We assume block fading for the interference channel gain, wherein it is a constant over a fusion interval. However, it may be appropriately changed depending on the actual operating band [55][56][57]. This also applies to other interference channel gains.  matrix C q 1|1 is usually set to the estimated covariance in the previous fusion interval. However, this information is not available to us and hence we set it to be 10-times scaled version of the covariance matrix used in the state model.

A. ANCHOR performance within a fusion interval
We first explore the performance of ANCHOR within a fusion interval and subsequently consider its performance across intervals. a) Convergence and optimized objective function: We assess the performance in the k-th fusion interval (i.e. from t k to t k+1 ). Note that the Bayesian CRB matrix B k and the Kalman estimation covariance matrix C k|k should be available as the inputs to the Bayesian tracking scheme in the k-th interval. We initialize them with an estimate ( i.e. an amplified version of the state covariance matrix), as is typical in Kalman filtering given that they will converge quickly after the initial phase [58]. Recall that the tracking performance depends on the improvement on the Bayesian CRB, which is evaluated by the objective function Q q=1 Λ T . Therefore, we compare the improvements of the objective value under three resource allocation schemes for the k-th fusion interval in Figure 7. For random allocation, we randomly generate the allocated resources satisfying the constraints of problem (26). The uniform allocation distributes the resources equally. The optimized allocation is designed by our proposed algorithm ANCHOR, in which the uniform allocation serves as initialization. Given ANCHOR is based on the alternating optimization method solving two subproblems, we also demonstrate the two nested iterations of the ANCHOR for the first outer iteration. We note that, based on the two nested iterations, the ANCHOR can increase the objective value monotonically as the outer iteration index increases. Compared to the other two allocations, the converged objective value of the optimized allocation is much larger. b) Radar resource allocation: Fig. 8 demonstrates the allocated radar resources for three methods (random, uniform, and ANCHOR). All values are normalized by the counterparts of the uniform allocation. Note that for the PARs with the optimized allocation, only target 1 is tracked by using radar 5 and 6. In fact, the designed dwell time T 4,1,k , T 4,2,k , T 5,2,k and T 6,2,k are close to zero numerically, which physically means that the corresponding radars stop scanning at some given time. While it is intuitive that the tracking performance is improved by increasing the power and dwell time as they potentially increase the SINR of the receiving measurements. However, in an interference scenario, we cannot increase all resources while satisfying various constraints and maintaining lower mutual interference. Thus, a trade-off is achieved via the optimized allocation (Fig. 8).  Fig. 9 shows the power and frequency allocation for the communications downlink transmissions, where the throughput margin is defined as the difference between the achieved and required throughput. Note that the constraint parameters are set to satisfy the uniform allocation in the simulations. Consequently, the throughput margins of the uniform allocation in Fig. 9a are all zero among the users. Interestingly, all users in the random allocation have margins to their throughput thresholds. Compared to the random allocation which only needs to satisfy the constraints of the optimization problem, the ANCHOR allocation optimizes the Bayesian CRB while being feasible to the constraint set. In practice, if the throughput margin is necessary, it is always achieved equivalently by increasing the throughput threshold for the optimized allocation probably at cost of a slight degeneration on tracking performance. d) State estimation performance: For the three allocations, we compare the tracking performance in Table III, where the calibrated root mean square error (CRMSE) is with N t being the number of the Monte Carlo trials, and the estimated target state is averaged over the N t states. The true target states at time t k+1 are given by T . Although the averaged target states are estimated well for all the three allocations, the optimized allocation achieves the smallest CRMSE on average, which reflects the estimation stability over the 500 trials. e) Performance across scenarios: To explore the statistical stability of the proposed allocation algorithm, we demonstrate the CRMSE performance over 15 different testing scenarios. In each scenario, the radars are randomly placed to track two targets with randomly generated initial states. The BS is also randomly placed and the communications users are within the BS coverage. The system parameters are initialized correspondingly by the same distributions used in the previous experiments. Fig. 10 shows the CRMSE of a single fusion interval, where the initial B k|k and C k|k are randomly generated for different scenarios. Each CRMSE value is averaged   (a) Allocation of communications power (left) and throughput margin (the difference between the achieved and required throughput) (right).  over 200 Monte Carlo trials. Clearly, over these scenarios, the average CRMSEs for the optimized allocation are roughly between 100 m and 800 m, which are further improved in the subsequent fusion intervals.

B. ANCHOR performance across fusion intervals
We conducted the experiments over 10 consecutive fusion intervals indexing from k = 1 to k = 10. Based on the allocated resources, we update B q k+1 Q q=1 and C k+1|k+1 Q q=1 of the k-th fusion interval as the inputs to the next interval. To evaluate the efficacy of the consecutive tracking scheme, we initialized B 1 and C 1|1 with a significantly amplified version of the state covariance matrix as a rough guess on the initial state. The tracking estimation is conducted based on (50) and the performance is averaged over 200 Monte Carlo trials. The trials are independent and each of them run the 10 consecutive fusion intervals completely. Fig. 11 demonstrates the averaged values over multiple fusion intervals. Starting with the same values, the optimized allocation improves the objective value significantly compared to the other two allocations. Besides, for the random allocation, it is expected that the objective value fluctuates over the different intervals.
In terms of the average CRB, the average CRMSE performance over 200 Monte Carlo trials is demonstrated in Fig. 11. The CRMSE values of the three allocations at the first fusion interval is extremely large and close to each other. Recall that we set B 1 and C 1|1 with large values to represent a rough initial guess, and thereby we infer that the first tracking performance is heavily influenced by this "bad" initial guess  and the improvement brought by the optimized allocation is negligible. Although the first tracking performance is not usable, all three CRMSE curves monotonically decrease with the fusion interval. This validates the efficacy of the Bayesian tracking scheme. The optimized allocation achieves the smallest CRMSE over all fusion intervals. Fig. 12 shows the tracking performance by inspecting the average deviations of location and velocity over all the intervals. Each average deviation is obtained by calculating the RMSE similar to (52) without matrix Λ. Compared to the uniform and random allocations, the optimized allocation consistently performs the best over all fusion intervals for both the location and velocity derivations; this is consistent with the result in Fig. 11 indicating the average CRMSE of the optimized allocation as being the lowest. In addition, it is interesting to note that the velocity derivation has a relatively large fluctuation at the second fusion interval compared to the monotonicity of the location derivation. To explain this phenomenon, recall that each radar can only observe the radial velocity in our system model and thereby an accurate estimation on target velocity relies on the information fusion from all radars. Each radar is expected to provide some "unique" information about the velocity, which requires a proper placement of the radars with respect to the tracking target. Otherwise, the velocity ambiguity cannot be reduced substantially and consequently leads to unsatisfying velocity estimation.

VI. SUMMARY
We considered spectrum sharing between heterogeneouslydistributed JRC with the goal of tracking multiple radar targets while maintaining the throughput levels for the communications downlinks. Within the asynchronous multi-target tracking framework, we proposed a Bayesian CRB-based metric for optimization, which further depends on the system resources (i.e. radar and communications power, dwell time, and shared bandwidth). The resulting resource allocation problem was non-convex and involved discrete and continuous variables. We solved this through our proposed ANCHOR algorithm based on the alternating optimization framework with guaranteed monotonicity. The frequency and power-time were allocated alternately. Our numerical experiments illustrated the key algorithmic and system aspects of resource allocation. We demonstrated that, compared to the trivial uniform and random allocations, the system performance is significantly improved by properly allocating the accessible heterogeneous resources of both radar and communications through ANCHOR algorithm. This study is helpful in meeting challenges of nextgeneration wireless systems where heterogeneous networks are envisaged such as different radio access technologies (RATs) and cell-free massive MIMO systems [59].
where H m i,q,k is the Jacobian matrix of h s m i,q,k , i on s q t k+1 based on s q t k+1 = f s m i,q,k , t k+1 − t m i,q,k . We have which is further decoupled into Q subproblems as minimize V q Tr V T qΛ T BqΛV q subject to Tr (V q ) = 1.
Problem (56) has the Lagrangian L (V q , λ q ) = Tr V T qΛ T B qΛ V q + λ q (Tr (V q ) − 1). Following the Lagrangian method yields ∂L ∂V q = 2Λ T BqΛV q + λqI = 0 ∂L ∂λq = Tr (V q ) − 1 = 0, (57) which has the solution Substituting the expression of V q into the objective function of problem (33), we have which is the objective function of problem (26). This completes the proof.

C. Proof of Lemma 3
Define w (z, V q ) = Tr V T qΛ T B qΛ V q andw (z) = min {V q } Q q=1 w (z, V q ), which is the objective function of problem (33). The function w (z, {V q }) is differentiable in z and the feasible set Z is convex. Furthermore, w (z, V q ) is strongly convex in V q and the set V q = {V q |Tr (V q ) = 1, V q 0} is compact. Thus, by applying the Danskin's theorem [60], we have ∇w (z) = Q q=1 ∂ V q w (z, V q ) | V q =V q , where V q = arg min Tr(V q )=1,V q 0 Tr V T qΛ T B qΛ V q . From this perspective, the update rules of Algorithm 3 can be merged into z n+1 = P Z (z n − η∇w (z n )). By applying [61,Theorem 31] with = 0 (due to the availability of the optimal V ), we arrive at the conclusion that the sequence {z n } will converge to the stationary point ofw (z) for z ∈ Z.

D. Proof of Theorem 4
At the -th iteration, the objective value of problem (26) is g (z , F c ). We have g (z , F c ) ≤ g z , F c +1 ≤ g z +1 , F c +1 , where the first inequality holds because of Algorithm 2, and the second inequality holds because Algorithm 3 converges to a stationary point. Thus, the sequence {g (z , F c )} is non-decreasing. Since the objective function g (z, F c ) is upper-bounded, the non-decreasing {g (z , F c )} will converge to a finite value.