A DFRC BS is serving $K$ single-antenna communication users and sensing $M$ radar targets as shown in Fig. 1. The BS is equipped with a fully-digital uniform linear array (ULA) with $N_{t}$ transmit antennas. Without loss of generality, the space between the adjacent elements for both the transmit array is $d_{0}=\frac{\lambda}{2}$, where $\lambda$ denotes the wavelength. The users and the targets are assumed to be sufficiently far from the BS, and thus the target can be viewed as points [6].

We assume $K+M$ baseband data streams transmitted from the BS for communication and sensing, denoted as $\mathbf{s}\in\mathbb{C}^{(K+M)\times 1}$, where $K$ streams are dedicated for communication, and all signals can be utilized for sensing. The transmitted signal $\mathbf{s}$ satisfies independent complex Gaussian distribution such that $\mathbb{E}\{\mathbf{s}\mathbf{s}^{H}\}=\mathbf{I}_{K+M}$. Then, the transmitted signal can be written as $\mathbf{x}=\mathbf{Ws}$, where $\mathbf{W}$ denotes the full digital transmit beamformer, which can be defined as,

The received signal at user $k$ can be represented as

where $\mathbf{h}_{k}\in\mathbb{C}^{N_{t}\times 1}$ denotes the real channel between the BS and user $k$, and $n_{k}\sim\mathcal{CN}(0,\sigma_{k}^{2})$ denotes the additive white Gaussian noise (AWGN) at user $k$. Accordingly, the received SINR at user $k$ can be derived as

with the received data rate at user $k$ obtained as

As for sensing, we denote $\theta_{m},\,m\in\{1,\dots,M\}$ as the direction of angle (DoA) of target $m$. The communication and sensing data streams can be simultaneously employed for sensing the targets. Then, the beampattern gain at target $m$ can be written as

where $\mathbf{a}(\theta_{m})=[1,e^{-j\frac{2\pi}{\lambda}d_{0}\sin(\theta_{m})},\dots,e^{-j\frac{2\pi}{\lambda}(N_{t}-1)d_{0}\sin(\theta_{m})}]$ is the steering vector towards angle $\theta_{m}$, $\mathbf{R}_{\mathbf{w}}=\mathbb{E}\{\mathbf{x}\mathbf{x}^{H}\}=\mathbf{W}\mathbf{W}^{H}$ is the covariance matrix of the transmitted signal.

The real channel $\mathbf{h}_{k}$ can be modeled as

where $\hat{\mathbf{h}}_{k}$ denotes the estimated channel, and $\Delta\mathbf{h}_{k}$ is the CSI error following the bounded error model as [16]

Because we only focus on the angle of direction (AoD) of the point-like target under far-field assumption, consider the deviation of the target AoD, and the steering vector of target $m$ can be modeled as

where $\theta_{m,min}$ and $\theta_{m,max}$ denote the minimum and the maximum values of possible $\theta_{m}$, respectively.

Under the existence of CSI and AoD uncertainty in the ISAC system, we aim to propose a framework to maximize the weighted sum of worst-case user data rate and the beampattern gain, by optimizing the digital beamforming such that

where $\rho$ is the weight allocated to communication function. Constraint (9a) denotes the power budget at the BS. The main challenges to solve this problem are i) the non-smooth objective function due to max-min problem; ii) the intractable CSI and angle uncertainty; iii) the non-convex expression of $r_{k}$ and $P(\theta_{m})$ as the objective function.

Focusing on communication, by introducing auxiliary variable $\delta_{k}$, (P1-a) can be further transformed as

where (12a) is equivalent to

where $\bm{\beta}=[\beta_{1},\dots,\beta_{K}]^{T}$ and $\bm{\nu}=[\nu_{1},\dots,\nu_{K}]^{T}$ are slack variables. The main challenge introduced by the communication measurement is to tackle the above three non-convex constraints.

Recalling (6), (13) can be rewritten as

where the quadratic expression with respect to variable $\mathbf{w}_{k}$ makes the constraint non-convex. Thus, this problem can be solved iteratively by successive convex approximation (SCA) method. The lower bound of the left-hand-side can be approximated by its first Taylor expansion, and (16) can be further transformed as

where $\mathbf{\Lambda}_{k}$, $\mathbf{b}_{k}$ and $c_{k}$ can be obtained by

where $\mathbf{w}_{k}^{(t)}$ denotes the results obtained at the last iteration.

Next, the S-Procedure is employed to deal with the bounded CSI error model, and (17) is equivalent to

where $\bm{\varphi}=[\varphi_{1},\dots,\varphi_{K}]$, $\varphi_{k}\geq 0$ are slack variables.

As for (14), first we apply Schur complement to transform (14) as

where $\mathbf{W}_{\bar{k}}$ is stacked as

Then, (22) can be further equivalently transformed as [17]

where $\bm{\xi}=[\xi_{1},\dots,\xi_{K}]$, $\xi_{k}\geq 0$ are the slack variables.

To address (15), we first rewrite it as

where the two slack variables at the right-hand-side of the inequality are decoupled by rewriting it as a difference of convex function. We further linearly approximate its upper bound, and then (15) is obtained as

Based on the above derivations, problem (P1-a) can be reformulated as

Focusing on (P1-b), to make the sensing beampattern gain more tractable, (5) can be rewritten as

where $\mathbf{A}(\theta_{m})=\mathbf{a}(\theta_{m})\mathbf{a}(\theta_{m})^{H}$. Because we only focus on sensing in this stage, temporarily ignoring communication metric and constant $(1-\rho)$, problem (P1-b) can be reduced as

where $\tilde{\mathbf{A}}(\theta_{m})=\{\mathbf{a}(\theta_{m})\mathbf{a}^{H}(\theta_{m})|\theta_{m}\in[\theta_{m,min},\theta_{m,max}]\},\forall m$. To make it tractable, we construct a convex hull of $\tilde{\mathbf{A}}(\theta_{m})$ as [18]

where $\mu_{m,s}$ is the weight, $S_{k}$ is the total number of samples, and $\mathbf{A}_{s}(\theta_{m})$ is the $s$-th sample. By focusing on the sum of beampattern gain in the objective function, and ignoring the constraints discussed in the last part, we have the following proposition.

Proposition 1: To maximize the the sum of minimization over the set $\tilde{\mathbf{A}}(\theta_{m})$ is equivalent to that over $\mathbf{B}_{m}$, which can be represented as

Proof: The proof can be referred to [18].

Thus, problem (P1-b1) is equivalent to

Then, we will optimize $\mathbf{W}$ and $\{\mu_{m,s}\}$ in two layers.

In the outer layer, assuming $\mathbf{W}$ is fixed, the reverse H$\ddot{\text{o}}$lder inequality can be applied as

where the equality holds if and only if $\frac{\mu_{m,1}^{\frac{1}{2}}}{\left(\text{tr}\left(\mathbf{A}_{s}(\theta_{m})\mathbf{R}_{\mathbf{w}}\right)\right)^{-1}}=\dots=\frac{\mu_{m,S}^{\frac{1}{2}}}{\left(\text{tr}\left(\mathbf{A}_{S}(\theta_{m})\mathbf{R}_{\mathbf{w}}\right)\right)^{-1}}$. Together with $\sum_{s=1}^{S}\mu_{m,s}=1,0\leq\mu_{m,s}\leq 1$, the $\mu_{m,s}$ that results in the worst case beampattern gain can be computed as

where $(t)$ denotes the optimal value obtained at the last iteration.

In the inner layer, with fixed $\{\mu_{m,s}\}$, (P1-b2) can be reformulated as

where $\bar{\mathbf{B}}_{m}=\sum_{s=1}^{S}\mu_{m,s}^{*(t)}\mathbf{A}_{s}(\theta_{m})$. To solve this non-convex problem, SCA method can be utilized again to approximate the objective function by its first Taylor expansion as

Then with the approximation, (P1-b3) is transformed to be a simple convex problem. Because (P1-a2) about communication is unrelated to updating $\mu_{m,s}$, we update beamforming matrix $\mathbf{W}$ in the inner layer, by combining the communication part solving (P1-a2) in section III.A with sensing part solving the transformed (P1-b3) as

which is an semidefinite programming (SDP) problem that can be efficiently solved by the existing toolbox such as CVX.

Overall, problem (P2) can be solved by two layers, where $\{\mu_{m,s}\}$ is optimized in the outer layer, and $\mathbf{W}$ is updated in the inner layer. The overall algorithm is summarized as Algorithm 1 and Algorithm 2.

The computational complexity of the proposed robust algorithm mainly depends on the complexity of optimizing $\mathbf{W}$ in Algorithm 1 and updating $\bm{\mu}$ in Algorithm 2. In the inner loop, solving problem (P2) contains variables $n$ of order $N_{t}M+(N_{t}+5)K$, $K+1$ quadratic constraints and $K$ linear matrix inequality constraints with sizes $N_{t}+1$ and $N_{t}+K+M$, respectively. Thus, the total computational complexity of solving problem (P2) is $\mathcal{O}(I_{I}\log_{2}(1/\epsilon)n(K(N_{t}+K+M)^{3}+nK(N_{t}+K+M)^{2}+n^{2}))$, where $I_{I}$ denotes the number of iterations of the inner loop, and $\epsilon$ denotes the accuracy of the SCA based method. In the outer loop, the main complexity depends on computing $\mu_{m,s}$ with (34), which is $\mathcal{O}(SN_{t}^{3})$. Finally, the overall computational complexity of the proposed Algorithm 2 can be represented as $\mathcal{O}(I_{O}(I_{I}\log_{2}(1/\epsilon)n(K(N_{t}+K+M)^{3}+nK(N_{t}+K+M)^{2}+n^{2})+SN_{t}^{3}))$, where $I_{O}$ is the number of iterations of the outer loop.