A balanced sensor scheduling for multitarget localization in a distributed multiple-input multiple-output radar network

Abstract

In this article, we consider the problem of optimally selecting a subset of transmitters from a transmitter set available to a multiple-input and multiple-output radar network. The aim is to minimize the location estimation error of underlying targets under a power constraint. We formulate it as a minimum-variance estimation problem and show that the underlying variance reduction function is submodular. From the properties of submodularity, we present a balanced selection policy which minimizes the worst-case error value using a minimax strategy. A greedy algorithm with guaranteed performance with respect to optimal solutions is given to efficiently implement the scheduling policy. The effectiveness and the efficiency of the proposed algorithm are demonstrated in simulated examples.

Introduction

The idea of multiple-input multiple-output (MIMO) radar was first introduced by Fishler et al.^1,2 Since then, it has been extensively investigated.^3–7 There are, broadly speaking, two categories of MIMO radars: collocated MIMO radars and distributed MIMO radars. We in this article focus on the latter, which has widely separated transmitters and receivers. It may perform much better than a conventional radar on information collection capability, ⁸ which comes from the synchronization of the network with a centralized processing capability and the utilization of the independency among multiple transmitter–receiver channels. ⁹ Comparing with a conventional radar, one pro of a distributed MIMO radar is its capability of attaining the spatial diversity gain. That is, the widely separated transmitters/receivers can provide rich information on the spatial scattering properties of a target since they are able to capture different cross sections of the target. To this end, target scattering scintillations might be avoided. Therefore, a distributed MIMO radar may outperform a conventional radar in the detection of low-observable targets. ¹⁰

In order to exploit the potential sensing capability of a distributed MIMO radar, it is of importance to develop a sophisticated signal processing and resource management system that can extract desired information from the MIMO radar as much as possible with various constraints. For example, which transmitters and receivers are involved in the operation have a great impact on the sensing capability of a distributed MIMO radar. For a multistatic radar network, its geometric topology significantly influences the performance of target localization. ⁹ The work of Godrich et al. ¹¹ showed that the performance of a distributed MIMO radar is related to signal-to-noise ratio (SNR), the number of transmitters and receivers involved, signal bandwidth, and so on. Accordingly, dynamic selection of transmitters/receivers of a distributed MIMO radar for maximizing its sensing performance is required. In practice, the network communication capability and power availability of a distributed MIMO radar restrict the number of transmitters/receivers that can be scheduled in a joint signal processing. The main challenge behind the selection of transmitters/receivers is to create a balance between the power that is consumed by transmitters/receivers and the amount of information acquired.

Much work has been done on the selection of transmitters/receivers for MIMO radars. Godrich et al. ⁶ used Cramér–Rao lower bound (CRLB) as an optimization metric, and considered two power allocation problems for a single target localization. One is to minimize the total power of transmitters subject to a predefined localization performance being met, while the power of each transmitter is restricted to an allowable range. The other is to maximize attainable localization performance with a given total power budget. Since both of the two problems are non-convex and non-linear, constraint relaxation and field decomposition were used to obtain approximate solutions. Godrich et al. ⁷ introduced two transmitter and receiver selection problems for a single target localization. Again CRLB was taken as a metric of optimization and two selection policies were considered. One is to minimize the number of transmitters and receivers subject to a given constraint on localization performance and the other is to select $k>0$ transmitters and receivers to minimize localization error. Both of the two problems were formulated as knapsack problems, which are typically NP-hard. Accordingly, greedy policy was proposed to achieve approximate solutions. Yan et al. ¹² considered the problem of joint beam selection and power allocation for multitarget tracking in an MIMO radar, where multiple beams transmit simultaneously. The Bayesian CRLB was used as the objective of the formulated mixed-integer non-linear optimization problem and a two-step solution based on relaxation and rounding was proposed to solve it. Yi et al. ¹³ considered the combinatorial problem of joint placement of transmitters and receivers for an MIMO radar and proposed a convex relaxation-based approximate algorithm. Xie et al. ¹⁴ considered the problem of node selection and power allocation for multitarget tracking in an MIMO radar, where the predicted conditional CRLB was used as the optimization criterion. A two-step semidefinite programming-based approach was proposed to derive an approximate solution. Xie et al. ¹⁵ considered the problem of receive-beam allocation when an MIMO radar works on “defocused transmit-focused receive” mode and proposed a convex relaxation-based approximate algorithm to solve the problem. Zhang et al. ¹⁶ considered the problem of joint subarray selection and power allocation for multitarget tracking with clutter in an MIMO radar, where the weighted sum of predicted conditional CRLB was used as the optimization criterion. A two-stage local search-based algorithm was proposed to solve the problem. Convex relaxation-based algorithms were proposed to solve the problem. Nosrati et al. ¹⁷ considered the problem of antenna selection for mitigating interference signals in an MIMO radar and proposed a relaxation-based approach to tackle the problem approximately. Yi et al. ¹⁸ proposed a reward-based iterative descending solution for joint beam and power scheduling in netted colocated MIMO radar systems. Yuan et al. ¹⁹ considered the power allocation problem in colocated MIMO radar systems based on the quality of service framework and proposed an iterative parallel search algorithm. While extensive work on selection of transmitters/receivers for an MIMO radar has been done,^6,7,12–19 the properties of the objective functions optimized were seldom considered in their work. The algorithms developed there were shown being effective as demonstrated by simulations but with no performance guarantee. This means that the performance of the algorithms could be very poor in some circumstances, which are not desired in safety-critical applications, such as aircraft surveillance.

Wang et al. ²⁰ considered the transmitter selection problem in an MIMO radar network for multitarget localization, which was formulated as an estimation variance minimization problem with a submodular variance reduction function. A weighted average policy (WAP) was used to optimize the selection of transmitters with a power budget constraint. Based on submodularity, a polynomial computational time algorithm with a guaranteed accuracy was proposed. Tohidi et al. ²¹ considered the problem of joint antenna and pulse placement for angle of arrival and velocity estimation in a colocated MIMO radar, where the CRLB for the angle and velocity estimator was used as the performance criterion and a submodularity-based algorithm with performance guarantee was proposed to solve the problem. Based on convex optimization, Yan et al. ²² developed two optimal resource allocation schemes for asynchronous multitarget tracking in heterogeneous radar networks. However, none of the abovementioned work considered the balanced selection policy, which is our focus in this article.

We consider the problem of selecting transmitters with balance between multiple targets under a power budget constraint for a distributed MIMO radar network consisting of multiple transmitters and receivers. We illustrate our problem in Figure 1, where two targets are localized by a distributed MIMO radar network consisting of five transmitters (For the sake of clarity, the receivers are not plotted.). Suppose that, because of a power budget constraint, we are allowed to select at most three transmitters to illuminate the two targets. The variance of the estimated target location is used as the performance metric. The right two plots of Figure 1 show the possible results of two different selection policies: unbalanced policy and balanced policy (BP). Although the total variance of the two targets obtained by unbalanced policy (top right plot) is less than that obtained by BP (bottom right plot), the variance of Target 1 under the former is higher than that under the latter, which is not desired. We seek to perform equally well with respect to the two targets. Formally speaking, the objective function of BP tries to maximize the least variance reduction (or, equivalently, minimize the greatest variance) using a maximin (minimax) strategy. After showing that the defined objective function under consideration is submodular, we then adapt the Saturation algorithm in Krause et al. ²³ to solve the problem. As WAP, the introduced greedy-based algorithm has a performance guarantee, that is, an upper bound on the least variance reduction is provided for the obtained solution of the algorithm. In summary, the innovation of this article is as follows: (1) a proof is provided that the variance reduction function derived based on CRLB is a non-decreasing submodular function and (2) a fast and balanced solution with performance guarantee is proposed for the problem of selecting transmitters under a power budget constraint in the presence of multiple targets for a distributed MIMO radar.

Figure 1.

Illustrations of balanced selection policy for a distributed MIMO radar. Left: a distributed MIMO radar for target localization. Top right: the result of unbalanced transmitter selection. Bottom right: the result of balanced transmitter selection.

The remainder of this article is organized as follows. The system model for target localization in a distributed MIMO radar network is described in section “Modeling.” There we derive the Fisher information matrix (FIM) of the underlying network parameterized by a target’s location. The transmitter selection problem in terms of BP objectives and the corresponding solutions are proposed in section “Problem formulation and solution.” We present the simulation results in section “Simulated scenarios” and conclude this article in section “Conclusion.”

Modeling

In this section, we describe the models we introduce for the distributed MIMO radar network and the measurements obtained by the MIMO radar network. We then provide a rigorous derivation of CRLB on the estimation of target localization. Based on CRLB, we define the variance reduction function as the performance metric as the transmitter selection and prove that it is submodular.

The MIMO radar network

Like the assumptions made in Godrich et al.,^6,7 we consider a distributed MIMO radar that consists of M transmitters and N receivers. We assume that all transmitters and receivers are time synchronized, widely separated, and located in a plane. Let ${\mathcal{S}}_{\mathrm{Tx}}:=\{(x_{m_{\mathrm{Tx}}},y_{m_{\mathrm{Tx}}}){\}}_{m=1}^M$ be the set of transmitter position, and ${\mathcal{S}}_{\mathrm{Rx}}:=\{(x_{n_{\mathrm{Rx}}},y_{n_{\mathrm{Rx}}}){\}}_{n=1}^N$ be the set of receiver position. Both ${\mathcal{S}}_{\mathrm{Tx}}$ and ${\mathcal{S}}_{\mathrm{Rx}}$ are assumed to be known. For the purpose of transmitter (receiver) scheduling and measurement processing, a centralized processing and control architecture is assumed to be used. That is, radar signals from all receivers are sent to the center and processed simultaneously. Denote $I(I\ge 2)$ by the number of targets in the region of interest, and $(x_i,y_i),i=1,2,\dots ,I$ by their locations. The variation of these locations, viewed by the transmitters and the receivers, is assumed to be small with respect to the system resolution capability. The latter is deemed to be far more better than a conventional radar alone in the sense that every receiver can independently process the returns from each of the individual transmitters which practically transmit waveforms almost orthogonal to each other. Denote $s_m(t)$ by a low-pass equivalent of each receiver, ${\mathcal{T}}_m$ by the duration of the mth transmitted signal, and ${\beta }_m$ by the waveform effective bandwidth. We assume that ${\int }_{{\mathcal{T}}_m}|s_m(t){|}^2\mathrm{dt}=1$ . As discussed in Lehmann et al. ²⁴ and Sun et al., ²⁵ an MIMO radar is able to localize targets and estimate unknown parameters, such as radar cross section (RCS) and propagation attenuation, of all targets in the region of interest through processing measurements collected by receivers.

Measurement model

The underlying MIMO radar network is designed as such that each of the N receivers is able to separate the mixture of transmitted signals by exploiting the orthogonality between the transmitted waveforms. Denote $R_{m{T}_x}^i$ and $R_{n{R}_x}^i$ by the ranges between transmitter m and target i, and between target i and receiver n, respectively. That is

$$\begin{array}{c}{R}_{m{T}_x}^i=\sqrt{{(x_{m_{\mathrm{Tx}}}-x_i)}^2+{(y_{m_{\mathrm{Tx}}}-y_i)}^2}\hfill \\ R_{n{R}_x}^i=\sqrt{{(x_i-x_{n_{\mathrm{Rx}}})}^2+{(y_i-y_{n_{\mathrm{Rx}}})}^2}\hfill \end{array}$$

(1)

Let ${\tau }_{m,i,n}(x_i,y_i)$ be the propagation time of a signal transmitted by transmitter m, reflected from target i, and received by receiver n. We have

$${\tau }_{m,i,n}(x_i,y_i)=\frac{R_{m{T}_x}^i+R_{n{R}_x}^i}{c}$$

(2)

where c is the light speed. The signal transmitted by transmitter m, and received by receiver n has the following baseband representation

$$\begin{array}{c}{z}_{m,n}\left(t\right)\hfill \\ =\sum _{i=1}^I\sqrt{{\alpha }_{m,i,n}{p}_{m_{\mathrm{Tx}}}}{h}_{m,i,n}{s}_m(t-{\tau }_{m,i,n})+v_{m,n}\left(t\right)\hfill \end{array}$$

(3)

with $p_{m_{\mathrm{Tx}}}$ being the power of the signal transmitted by transmitter m and ${\alpha }_{m,i,n}\propto (1/R_{m_{\mathrm{Tx}}}^i{R}_{n_{\mathrm{Rx}}}^i{)}^2$ being the path attenuation in the signal magnitude for target i. The term $h_{m,i,n}$ in equation (3) incorporates the effect of phase offsets between transmitter m, receiver n, and target i’s RCS impact on phase and amplitude. We assume that $h_{m,i,n}$ is deterministic and unknown. In addition, $h_{m,i,n}(1\le m\le M,1\le i\le I,1\le n\le N)$ are assumed to be mutually independent. $v_{m,n}(t)$ is circularly symmetric complex Gaussian noise with mean being zero. Meanwhile, it is spatially and temporally white and its autocorrelation function is ${\sigma }_w^2\delta (\tau )$ . The propagation path among transmitter m, target i, and receiver n is denoted as channel $(m,i,n)$ , of which the SNR may be expressed as^6,7

$$\mathrm{SN}{\mathrm{R}}_{m,i,n}\stackrel{\mathrm{\Delta }}{=}\frac{{\alpha }_{m,i,n}{\Vert h_{m,i,n}\Vert }^2{p}_{m_{\mathrm{Tx}}}}{{\sigma }_w^2}$$

(4)

Let $\mathit{u}=[{\mathit{u}}_1^T,\dots ,{\mathit{u}}_I^T{]}^T$ denote the state of all targets, where ${\mathit{u}}_i=[x_i,y_i,{\mathit{h}}_i^T{]}^T$ , and ${\mathit{h}}_i^T=[h_{1,i,1},h_{1,i,2},\dots ,h_{M,i,N}{]}^T$ . Upon receiving the signal represented by equation (3), the conditional probability density function (PDF) for measurement $z_{m,n}$ corresponding to transmitter m and receiver n can be written as^6,7

$$\begin{array}{c}f(z_{m,n}|\mathit{u})\hfill \\ =\frac{1}{{\left(2\pi {\sigma }_w^2\right)}^{\frac{1}{2}}}\exp \{-\frac{1}{{\sigma }_w^2}{\int }_{{\mathcal{T}}_m}|z_{m,n}\left(t\right)\hfill \\ -\sum _{i=1}^I\sqrt{{\alpha }_{m,i,n}{p}_{m_{\mathrm{Tx}}}}{h}_{m,i,n}{s}_m(t-{\tau }_{m,i,n}){|}^2\mathrm{dt}\}\hfill \\ =\frac{1}{{\left(2\pi {\sigma }_w^2\right)}^{\frac{1}{2}}}\exp \{-\frac{1}{{\sigma }_w^2}q(z_{m,n}|\mathbf{u})\}\hfill \end{array}$$

(5)

where

$$\begin{array}{c}q(z_{m,n}|\mathit{u})\hfill \\ ={\int }_{{\mathcal{T}}_m}|z_{m,n}\left(t\right)-\sum _{i=1}^I\sqrt{{\alpha }_{m,i,n}{p}_{m_{\mathrm{Tx}}}}{h}_{m,i,n}{s}_m(t-{\tau }_{m,i,n}){|}^2\mathrm{dt}\hfill \end{array}$$

(6)

Under a centralized measurement processing consideration as aforementioned, we denote the measurement vector at a specific time t as $\mathit{z}(t)=[z_{1,1}^T,z_{1,2}^T,\dots ,z_{M,N}^T{]}^T$ . Similarly, the following vector forms of the parameters aforementioned are used. $\mathit{\beta }=[{\beta }_1,{\beta }_2,\dots ,{\beta }_M{]}^T$ , ${\mathbf{\alpha }}_i=[{\alpha }_{1,i,1},{\alpha }_{1,i,2},\dots ,{\alpha }_{M,i,N}{]}^T$ , ${\mathit{p}}_{\mathrm{Tx}}=[p_{1_{\mathrm{Tx}}},\dots ,p_{M_{\mathrm{Tx}}}{]}^T$ , and ${\mathbf{\tau }}_i=[{\tau }_{1,i,1},{\tau }_{1,i,2},\dots ,{\tau }_{M,i,N}{]}^T$ .

As in Wang et al., ²⁰ our aim is to estimate the location of all targets in the region of interest based on the measurements using a portion of transmitters/receivers. We address this problem via minimizing the variance of localization error subject to the given constraint on power. Next, we derive a variance reduction function based on CRLB and show that it is submodular. Accordingly, the optimization problem can be efficiently addressed using the properties of submodularity.

The variance reduction function based on CRLB

It is understood from information theory that CRLB provides a lower bound for the mean square error (MSE) on parameter estimation via an unbiased estimator. The MSE of a maximum likelihood estimator (MLE) may approach the CRLB as SNR is large (say, over 10 dB).^11,26 Therefore, we derive the CRLB of target localization using measurements from the MIMO radar as a function of transmitter power allocation and it is served as the cost function for our optimization problem.

Given a vector parameter $\mathit{u}$ and an measurement vector $\mathit{z}(t)$ , any unbiased estimate $\hat{\mathit{u}}$ of $\mathit{u}$ satisfies the following inequality

$${\mathbf{E}}_{\mathit{u}|\mathit{z}\left(t\right)}\left\{(\hat{\mathit{u}}-\mathit{u}){(\hat{\mathit{u}}-\mathit{u})}^T\right\}\ge J_{\mathit{z}}^{-1}\left(\mathit{u}\right)$$

(7)

where $J_{\mathit{z}}^{-1}(\mathit{u})$ is the CRLB, and $J_{\mathbf{z}}(\mathbf{u})$ is the FIM which is defined as

$$J_{\mathbf{z}}\left(\mathbf{u}\right)={\mathbf{E}}_{\mathbf{z}}\left\{{\nabla }_{\mathbf{u}}\log f(\mathbf{z}|\mathbf{u}){\nabla }_{\mathbf{u}}^T\log f(\mathbf{z}|\mathbf{u})\right\}$$

(8)

where $f(\mathbf{z}|\mathbf{u})$ is the likelihood function of $\mathbf{z}$ conditional on $\mathbf{u}$ .

By assuming that different transmitter–target–receiver channels are independent, the FIM for target i with respect to $x,y$ is derived as

$$J_i^{x,y}=\frac{8{\pi }^2}{c^2{\sigma }_w^2}\sum _{m=1}^M\sum _{n=1}^N{\beta }_m^2{p}_m{\alpha }_{m,i,n}{\left|h_{m,i,n}\right|}^2\left[\begin{array}{cc}{u}_a\hfill & u_c\hfill \\ \multicolumn{1}{c}{u_c}& u_b\hfill \end{array}\right]$$

(9)

where

$$\begin{array}{c}{u}_a={(\frac{x_{m{T}_x}-x_i}{R_{m{T}_x}^i}+\frac{x_{n{R}_x}-x_i}{R_{n{R}_x}^i})}^2\hfill \\ u_b={(\frac{y_{m{T}_x}-y_i}{R_{m{T}_x}^i}+\frac{y_{n{R}_x}-y_i}{R_{n{R}_x}^i})}^2\hfill \\ u_c=(\frac{x_{m{T}_x}-x_i}{R_{m{T}_x}^i}+\frac{x_{n{R}_x}-x_i}{R_{n{R}_x}^i})(\frac{y_{m{T}_x}-y_i}{R_{m{T}_x}^i}+\frac{y_{n{R}_x}-y_i}{R_{n{R}_x}^i})\hfill \end{array}$$

(10)

The derivation of equation (9) is given in Appendix 1.

Therefore, denoted by ${\mathbf{C}}_{x_i,y_i}$ , the CRLB with respect to each target i’s location is written as

$${\mathbf{C}}_{x_i,y_i}={\left\{J_i^{x,y}\right\}}^{-1}$$

(11)

Given the set of transmitters ${\mathbf{T}}_x:=\{1,2,\dots ,M\}$ and the set of receivers ${\mathbf{R}}_x:=\{1,2,\dots ,N\}$ , for $A_T\subseteq {\mathbf{T}}_x$ and $A_R\subseteq {\mathbf{R}}_x$ , the CRLB with respect to target i’s location by selecting $A_T$ and $A_R$ may be denoted by ${\mathbf{C}}_{x_i,y_i}(A_T,A_R)$ . In other words, ${\mathbf{C}}_{x_i,y_i}(A_T,A_R)=\{J_i^{x,y}(A_T,A_R){\}}^{-1}$ , where

$$\begin{array}{c}{J}_i^{x,y}(A_T,A_R)\hfill \\ =\frac{8{\pi }^2}{c^2{\sigma }_w^2}\sum _{m\in A_T}\sum _{n\in A_R}{\beta }_m^2{p}_m{\alpha }_{m,i,n}{\left|h_{m,i,n}\right|}^2\left[\begin{array}{cc}{u}_a\hfill & u_c\hfill \\ \multicolumn{1}{c}{u_c}& u_b\hfill \end{array}\right]\hfill \end{array}$$

(12)

We define the variance reduction function $f_i(A_T,A_R)$ with respect to $A_T,A_R$ as

$$f_i(A_T,A_R)\stackrel{\mathrm{\Delta }}{=}\mathrm{trace}\left({\mathbf{C}}_{x_i,y_i}(Ø,Ø)\right)-\mathrm{trace}\left({\mathbf{C}}_{x_i,y_i}(A_T,A_R)\right)$$

(13)

We next show that the objective functions equation (13) satisfy the following intuitive diminishing returns property: adding a pair of transmitter-receiver reduces variance more if we have assigned few pairs of transmitter-receiver so far, and less if we already assigned lots of pairs of transmitter-receiver. This intuition can be formalized using the combinatorial concept of submodularity.^27,28 More specifically, let E be a ground set. A set function $f:2^E\to \mathcal{R}$ is called submodular, if for all $A\subseteq B\subseteq E$ and $j\in E\backslash B$

$$f(A\cup \left\{j\right\})-f\left(A\right)\ge f(B\cup \left\{j\right\}))-f\left(B\right)$$

Theorem 1

The objective function defined by equation (13) is a non-decreasing submodular function. The proof is given in Appendix 1.

Problem formulation and solution

The balanced transmitter selection problem

In the final section, variance reduction of target localization estimation as a function of transmitting signal power is derived based on CRLB, which serves as a minimization cost function for selecting a subset of transmitters and receivers. Such a selection/scheduling is subject to a total transmitting power constraint. For simplicity, we assume that all receivers are activated for all time to receive signals, although the proposed algorithms can be extended to the case of selecting transmitters and receivers simultaneously.

In Wang et al., ²⁰ WAP, which optimizes an average localization performance for all targets through weighting, was proposed for selecting a subset of transmitters which maximizes the objective function under the total power constraint. As aforementioned, if we optimize for the average performance with WAP, it can happen that a few of the targets are poorly localized, that is, there may exist some targets whose localization error variance reduction is insignificant, which can be problematic for security-critical applications. Alternatively, we may optimize variance reduction of the targets based on a balanced selection. Specifically, we maximize the minimum of each target’s variance reduction subject to the total power constraint. That is

$$\begin{array}{c}\underset{A_T\subseteq {\mathbf{T}}_x}{max}\underset{i\in I}{min}{f}_i\left(A_T\right)\hfill \\ f_i\left(A_T\right)=\mathrm{trace}\left({\mathbf{C}}_{x_i,y_i}(Ø,Ø)\right)-\mathrm{trace}\left({\mathbf{C}}_{x_i,y_i}(A_T,{\mathbf{R}}_x)\right)\hfill \\ \mathrm{s}.\mathrm{t}.\sum _{m\in A_T}{p}_m\le p_{\mathrm{total}}\hfill \end{array}$$

(14)

In equation (14), the objective function of BP tries to maximize the least variance reduction using a maximin strategy. In this way, we seek to perform equally well with respect to each target. As the robust submodular observation selection (RSOS) problem in Krause et al. ²³ (we will next describe it in detail) which is proved to be NP-hard, equation (14) is an NP-hard problem. Therefore, approximate solutions are sought. In particular, approximate algorithms with polynomial-time computational complexity and performance guarantee are preferred.

A balanced subset selection policy

Although each $f_i$ in equation (14) is submodular, in general, the function $\underset{i}{min}{f}_i(S)$ is not submodular any more. Therefore, the theory of Nemhauser et al. ²⁹ is not applied, that is, the greedy algorithm of Nemhauser et al. ²⁹ cannot obtain an approximate solution with a performance guarantee. Following the idea of Krause et al., ²³ we relax the knapsack constraint in equation (14) and obtain a relaxed version of equation (14), that is

$$\begin{array}{c}{max}_{c,A_T\subseteq {\mathbf{T}}_x}c,\hfill \\ \mathrm{s}.\mathrm{t}.f_i\left(A_T\right)\ge c,i=1,\dots ,I,\sum _{m\in A_T}{p}_m\le \alpha p_{\mathrm{total}}\hfill \end{array}$$

(15)

where $\alpha \ge 1$ is the relaxation factor, and recall that $f_i=\mathrm{trace}({\mathbf{C}}_{x_i,y_i}(Ø,Ø))-\mathrm{trace}({\mathbf{C}}_{x_i,y_i}(A_T,{\mathbf{R}}_x))$ . If $\alpha =1$ , equation (15) is the same as equation (14).

Next, we adapt the Saturation algorithm in Krause et al. ²³ to solve equation (15), and propose a subset selection algorithm for transmitters as in Algorithm 1.

Algorithm 1. Algorithm pseudo code for equation (15).
Input: $f_i,p_{\mathrm{total}}>0,U={\mathbf{T}}_x$ ;
Output $A_{\mathrm{best}}\subseteq U$ ;
1: $c_{min}\leftarrow 0$ ; $c_{max}\leftarrow \underset{i}{min}{f}_i({\mathbf{T}}_x)$ ; $A_{\mathrm{best}}\leftarrow Ø$ ;
2: while $(c_{max}-c_{min})\ge \epsilon$ do
3:   $c\leftarrow (c_{min}+c_{max})/2$ ;      ${\overline{F}}_c(A)\leftarrow \frac{1}{I}\sum _{i}min\{f_i(A),c\}$ ;
4:   $\hat{A}\leftarrow Ø$ ;
5:  while ${\overline{F}}_c(\hat{A})<c$ do
6:    $g=\arg \underset{m\in {\mathbf{T}}_x\backslash \hat{A}}{max}\frac{{\overline{F}}_c(A\cup \{m\})-{\overline{F}}_c(A)}{p_m}$ ;
7:    $\hat{A}\leftarrow \hat{A}\cup \{g\}$ ;
8:  end while
9:  if $\sum _{m\in \hat{A}}{p}_m>\alpha p_{\mathrm{total}}$ then
10:    $c_{max}\leftarrow c$ ;
11:  else
12:    $c_{min}\leftarrow c;A_{\mathrm{best}}=\hat{A}$ ;
13:  end if
14: end while
15: return $A_{\mathrm{best}}$

In the vein of the Saturation algorithm, ²³ in Algorithm 1, c is an upper bound to equation (14). The lower bound $c_{min}$ and the upper bound $c_{max}$ are improved sequentially by bisection search to find “best”c. At the beginning, we initiate the upper bound to be $\underset{i}{min}{f}_i({\mathbf{T}}_x)$ and the lower bound to be $0$ . At each iteration, we check the feasibility of $c=(c_{min}+c_{max})/2$ . To do it, for each target i, we define a truncation function of the submodular function $f_i$ as $min\{f_i(A),c\}$ and then define ${\overline{F}}_c(A)=1/I\sum _{i}min\{f_i(A),c\}$ . According to the property of submodular functions, both the function $min\{f_i(A),c\}$ and the function ${\overline{F}}_c(A)$ are submodular. Therefore, a solution $\hat{A}$ can be obtained using a greedy algorithm with performance guarantee. If c is not feasible, that is, $\sum _{m\in \hat{A}}{p}_m>\alpha p_{\mathrm{total}}$ , then let $c_{max}=c$ ; otherwise, save $\hat{A}$ as the current best solution $A_{\mathrm{best}}$ and let $c_{min}=c$ . Finally, Algorithm 1 guarantees to find a solution $A_{\mathrm{best}}$ , which is a superset of the optimal solution of equation (14), achieves a value greater than or equal to the optimal value of equation (14) within a small approximation error. ²³ Following Krause et al., ²³ we have the following theorem for Algorithm 1. The correctness of Theorem 2 can be proved by applying Theorem 5 in Krause et al. ²³ to equation (14).

Theorem 2

For the objective function defined in equation (14) and $\alpha =1+\log (\underset{m\in {\mathbf{T}}_x}{max}\sum _i{f}_i(\{m\}))$ , Algorithm 1 returns a solution $A_{\mathrm{best}}$ , such that

$$\begin{array}{c}\multicolumn{1}{c}{\underset{i}{min}{f}_i\left(A_{\mathrm{best}}\right)\ge \underset{\sum _{m\in A}{p}_m\le p_{\mathrm{total}}}{max}\underset{i}{min}{f}_i\left(A\right)}\\ \sum _{m\in A_{\mathrm{best}}}{p}_m\le \alpha p_{\mathrm{total}}\hfill \end{array}$$

(16)

The computational complexity of Algorithm 1 is $\mathcal{O}(|{\mathbf{T}}_x{|}^{2}I\log (I\underset{i}{min}{f}_i({\mathbf{T}}_x)))$ .

More specifically, for any real number $\alpha \ge max\{1,1+\log (\underset{m\in {\mathbf{T}}_x}{max}\sum _i{f}_i(\{m\}))\}$ , the solution returned by Algorithm 1 is no less than the optimum value of equation (14). That is, the value obtained by Algorithm 1 provides an upper bound of the optimum value of equation (14). Clearly, the smaller the value of $\alpha$ , the tighter the upper bound is. Note that, $1+\log (\underset{m\in {\mathbf{T}}_x}{max}\sum _i{f}_i(\{m\}))$ is the theoretical lower bound that $\alpha$ takes to guarantee the correctness of Theorem 2. Based on our observation in the following simulation, a reasonable solution can be found by Algorithm 1 even when $\alpha =1$ . In addition, the algorithm runs less time than that calculated theoretically.

Numerical simulation and analysis

In this section, we demonstrate the performance of Algorithm 1 considering two simulated scenarios on multitarget localization of a distributed MIMO radar. Under a small scale of an MIMO radar, we are able to compare the performance of Algorithm 1 with that of the optimal solutions. We also compare the performance differences between Algorithm 1, WAP, ²⁰ and the particle swarm optimization (PSO) algorithm.^30–32

PSO, which conducts search using a population of particles, is a kind of heuristic algorithms that can be applied to minimax optimization problems. ³⁰ To ensure an improved convergence of particles, we adopt the binary PSO algorithm in Nezamabadi-pour et al. ³¹ to solve equation (14). In the initialization of the binary PSO algorithm, a population of particles is generated. Each particle represents a potential solution, that is, a binary code, of which the length is the total number of transmitters M. The mth bit of the code is randomly set as 0 or 1. The transmitter m is selected if the mth bit of the code is 1; otherwise, the transmitter m is not selected. For each iteration, the following steps are performed until a termination condition, such as the maximum number of iterations is reached, is met. (1) The objective function $\underset{i\in I}{min}{f}_i(A_T)$ is used to evaluate each particle and the fitness value of each particle is derived. (2) For each particle, update its best fitness value so far, namely, the local optimal value, and the corresponding code. (3) Find the best particle of the population so far, that is, the particle with the maximum fitness value over all iterations so far, namely, the global optimal particle. (4) Update the velocity of each particle based on the current code of this particle, the corresponding local optimal code, the global optimal particle, and the parameters (inertia weight, cognitive parameter, and social parameter). The velocity of each particle is limited in the range $(v_{min},v_{max})$ , which is then converted to a real number $a\in [0,1]$ using the function $2\times |1/1+e^{-v}-0.5|$ . Accordingly, the code of the particle is changed with probability a. The parameter settings of the PSO algorithm are as follows. The number of particles is 50. The maximum number of iterations is 100. The inertia weight w gradually decreased from $0.6$ to $0.1$ . Both the cognitive parameter and the social parameter are set as 2. Note that there are many variants of PSO, for example, the adaptive PSO algorithm with supervised learning and control that adaptively chooses parameters, ³³ the comprehensive learning PSO algorithm with local search that enhances the diversity of particles and global search ability, ³⁴ and the surrogate-assisted PSO algorithm with Pareto active learning that deals with computationally expensive objective functions. ³⁵ In addition, a hyper-heuristic framework can save computational time and find high-quality solutions more efficiently. ³⁶ The ideas of the abovementioned work^33–36 might be used to improve the performance of the binary PSO algorithm in this article, and we leave it as our future work.

Performance of the three algorithms is measured in terms of localization accuracy and computational efficiency. Note that in the simulations, once we derive the results of variance reduction, we then calculate the corresponding variance (i.e. $\mathrm{trace}({\mathbf{C}}_{x_i,y_i}),i=1,2$ ) for easy of exposition.

Simulated scenarios

The performance of Algorithm 1 will be checked via Scenario 1, by which we also compare the performance of Algorithm 1 with WAP and PSO. Statistical results of Algorithm 1 and PSO are generated via Scenario 2.

• Scenario 1: the MIMO radar consists of $M=12$ transmitters and $N=10$ receivers. The power of each transmitter is $2\mathrm{kW}$ . The maximum total transmitting power is constrained not to exceed $12\mathrm{kW}$ . As illustrated in Figure 2, the locations of each transmitters $T_1~T_{12}$ are $(-3\times {10}^5,2\times {10}^5)$ , $(-4.73\times {10}^5,1\times {10}^5)$ , $(-4.73\times {10}^5,-1\times {10}^5)$ , $(-3\times {10}^5,-2\times {10}^5)$ , $(-1\times {10}^5,2\times {10}^5)$ , $(1\times {10}^5,2\times {10}^5)$ , $(-1\times {10}^5,-2\times {10}^5)$ , $(1\times {10}^5,-2\times {10}^5)$ , $(3\times {10}^5,2\times {10}^5)$ , $(4.73\times {10}^5,1\times {10}^5)$ , $(4.73\times {10}^5,-1\times {10}^5)$ , and $(3\times {10}^5,-2\times {10}^5)$ , respectively. The coordinates of each receivers $R_1~R_{10}$ are $(-4\times {10}^5,1.73\times {10}^5)$ , $(-5\times {10}^5,0)$ , $(-4\times {10}^5,-1.73\times {10}^5)$ , $(-2\times {10}^5,2\times {10}^5)$ , $(-2\times {10}^5,-2\times {10}^5)$ , $(2\times {10}^5,2\times {10}^5)$ , $(2\times {10}^5,-2\times {10}^5)$ , $(4\times {10}^5,1.73\times {10}^5)$ , $(5\times {10}^5,0)$ , and $(4\times {10}^5,-1.73\times {10}^5)$ , respectively. Two targets appear at locations $(-3\times {10}^5,0)$ and $(3\times {10}^5,0)$ , respectively. A lossless path and uniform scattering are assumed for the channels of transmitter–Target 1–receiver, that is, $|h_{m,1,n}|=1$ , $1\le m\le M$ , and $1\le n\le N$ . For the channels of transmitter-Target 2-receiver, we assume that the path loss exists and is measurable, that is, $|h_{m,2,n}|=l·|h_{m,1,n}|$ , $1\le m\le M$ , and $1\le n\le N$ , where $0\le l\le 1$ is a decay factor.

• Scenario 2: $100$ Monte Carlo runs are carried out in each of the following three cases:

• Case 1: assume there are $M=18$ transmitters and $N=9$ receivers and their distribution is shown in Figure 3. The power of each transmitter is $2\mathrm{kW}$ . There are five targets which are uniformly distributed in the area. The characteristic of each transmitter–target–receiver channel is the same. The maximum total transmitting power $p_{\mathrm{total}}$ is constrained to $12,140,160,\mathrm{and}180\mathrm{kW}$ .

• Case 2: we assume that the number of transmitters is $16,17,18,19,\mathrm{and}20$ . The transmitters are distributed evenly on a circle with radius $r=2\times {10}^5\mathrm{m}$ and their azimuths are $2\pi /\mathrm{Mm},1\le m\le M$ , respectively. The number and distribution of the receivers are the same as those in Case $1$ . The power of each transmitter is $2\mathrm{kW}$ . There are five targets, and they are uniformly distributed in the surveillance area. The characteristic of each transmitter–target–receiver channel is the same. The total transmitting power $p_{\mathrm{total}}$ is constrained to $12\mathrm{kW}$ .

• Case 3: assume that the configuration of the MIMO radar network is the same as that in Case $1$ , except that $p_{\mathrm{total}}\le 12\mathrm{kW}$ , and there are $6,7,8,9,\mathrm{and}10$ targets which are uniformly distributed in the area.

Figure 2.

Illustrations of locations of transmitters, receivers, and targets in Scenario $1$ . The red solid circles represent the targets, the blue rectangles represent the transmitters, and the black diamonds represent the receivers.

Figure 3.

Illustrations of locations of transmitters and receivers in Case $1$ of Scenario $2$ . The blue rectangles represent the transmitters, and the black diamonds represent the receivers.

Simulation results and analysis

The observations and analysis of the simulations are summarized as follows.

BP

In Scenario 1, the signal decay factor l is assigned to $1,0.8,0.6,0.4,0.3,\mathrm{and}0.2$ , respectively. In other words, the detection of Target $2$ is getting harder and harder as l is getting smaller. We run Algorithm 1 for BP. Let the relaxation factor $\alpha =1$ , that is, the constraint in equation (15) does not change. The selected transmitters of Algorithm 1 are shown in Figure 4. The derived variances (instead of variance reductions) of the two targets (i.e. $\mathrm{trace}({\mathbf{C}}_{x_i,y_i}),i=1,2$ ) under each value of l are presented in Table 1.

Figure 4.

Selected transmitters of optimal solution (OPT), Algorithm 1 (Alg 1), WAP, and PSO in Scenario 1: (a) $l=1$ , (b) $l=0.8$ , (c) $l=0.6$ , (d) $l=0.4$ , (e) $l=0.3$ , and (f) $l=0.2$ .

Table 1.

The variances of the two targets under each value of l obtained by OPT, Alg 1, WAP, and PSO in Scenario 1.

l	OPT	Alg 1	WAP	PSO
1.0	(5.18, 5.18)	(6.24, 6.24)	(5.18, 5.18)	(6.27, 6.27)
0.8	(7.25, 7.25)	(9.09, 8.86)	(5.18, 8.09)	(7.25, 7.25)
0.6	(10.2, 11.1)	(11.7, 12.9)	(7.59, 11.4)	(10.2, 11.1)
0.4	(20.4, 23.1)	(14.1, 23.6)	(7.59, 25.6)	(20.4, 23.1)
0.3	(29.8, 39.1)	(20.4, 41.1)	(7.59, 45.4)	(20.8, 39.1)
0.2	(29.8, 87.9)	(29.8, 87.9)	(13.8, 94.4)	(29.8, 87.9)

OPT: optimal solution; Alg 1: Algorithm 1; WAP: weighted average policy; PSO: particle swarm optimization.

From Table 1 and Figure 4, we observe that, when $l=1$ , transmitters $T_2,T_4$ , and $T_5$ which gain a better observation on Target $1$ and transmitters $T_6,T_{10}$ , and $T_{12}$ which gain a better observation on Target $2$ are selected by Algorithm 1. Such a selection yields a balanced observation of the two targets when the characterization of all channels is the same. Subsequently, the variances for the two targets are the same.

In the case that $l=0.8$ , which signifies that the network observability of Target $2$ is getting worse, transmitters $T_3$ and $T_7$ which have a better observation on Target $1$ and transmitters $T_6,T_8,T_{10},\mathrm{and}{T}_{12}$ which may yield a better observation on Target $2$ are selected by Algorithm 1. More transmitting power resources are assigned to Target $2$ to keep a balanced observation of the two targets as a result of signal decay. As l is decreasing, which means that the network observability of Target $2$ becomes worse and worse, more and more transmitting power resources are assigned to Target $2$ by Algorithm 1. For example, when $l=0.2$ , the selected transmitters are $T_6,T_8,T_9,T_{10},T_{11},\mathrm{and}{T}_{12}$ , all of which are close to Target $2$ . These results validate Algorithm 1 for BP.

It is worth mentioning that since $\alpha =1$ , the performance guarantee of Theorem 2 is no longer valid, that is, the variance reduction value obtained by Algorithm 1 is not necessarily to be greater than the optimum, which is evidenced in Table 1. Nevertheless, compared with the optimum, Algorithm 1 is still able to yield very good approximations. On the other side, simulation results show that if the relaxation factor $\alpha$ satisfies the assumption of Theorem 2, Algorithm 1 indeed obtains the approximation guarantee provided by Theorem 2. For saving of space, the relevant simulation results are not presented in this article.

Comparison of Algorithm 1 with WAP and PSO

We run WAP and PSO on Scenario $1$ . For WAP, the weights of the two targets are assumed to be equal to unity. The results are also shown as in Figure 4 and Table 1. Clearly, the selection of transmitter sets for each l is different comparing with those selected by Algorithm 1. For example, transmitters $T_1,T_2,\mathrm{and}{T}_3$ close to Target $1$ and transmitters $T_{10},T_{11},\mathrm{and}{T}_{12}$ close to Target $2$ are selected by WAP under $l=1$ . In the same scenario, $T_2,T_4,T_5,T_6,T_{10},\mathrm{and}{T}_{12}$ are selected by Algorithm 1. From Table 1, it is seen that Algorithm 1 does not perform better than WAP when $l=1$ . This is because that in this case, the two targets are totally symmetric. WAP, which aims to maximize total variance reduction of the two targets, also achieves balanced variance reduction by chance. As the decreasing of l, the observability of Target 2 becomes worse and worse, more transmitters close to Target 2 may be selected by WAP. However, Table 1 shows that the results obtained by Algorithm 1 are more balanced comparing with WAP when $l=0.4,0.3,\mathrm{and}0.2$ . PSO performs surprisingly well on Scenario $1$ except that when $l=1$ . However, it is noted that PSO does not have performance guarantee and its performance depends on the initialization of particles. Next, we will give a statistical performance comparison of Algorithm 1 and PSO.

Statistical performance

As a way to analyze the robustness of Algorithm 1, we perform Monte Carlo multiple runs on the three cases in Scenario $2$ , where at each case multitarget locations are randomly generated from a uniform distribution within the region of interest. We also run PSO on Scenario $2$ . Tables 2–4 show the statistical results, where $A_{\mathrm{mean}}^{\mathrm{Alg}1}$ , $A_{\mathrm{mean}}^{\mathrm{PSO}}$ , $T_{\mathrm{mean}}^{\mathrm{Alg}1}$ , $T_{\mathrm{mean}}^{\mathrm{PSO}}$ , and $T_{\mathrm{OPT}}$ represent the average approximation factor and the average run-time achieved by Algorithm 1 and PSO, and the run-time of the exhaustive enumeration algorithm (to provide the optimal solution), respectively. Specifically, we define

$$\begin{array}{c}{A}_{\mathrm{mean}}^{\mathrm{Alg}1}=\frac{1}{K}\sum _{k=1}^K\frac{\underset{i=1,2,\dots ,I}{max}\mathrm{trace}\left(C_{x_i,y_i}\left(S_{\mathrm{Alg}1}^k\right)\right)}{\underset{i=1,2,\dots ,I}{max}\mathrm{trace}\left(C_{x_i,y_i}\left(S_{\mathrm{OPT}}^k\right)\right)}\hfill \\ A_{\mathrm{mean}}^{\mathrm{PSO}}=\frac{1}{K}\sum _{k=1}^K\frac{\underset{i=1,2,\dots ,I}{max}\mathrm{trace}\left(C_{x_i,y_i}\left(S_{\mathrm{PSO}}^k\right)\right)}{\underset{i=1,2,\dots ,I}{max}\mathrm{trace}\left(C_{x_i,y_i}\left(S_{\mathrm{OPT}}^k\right)\right)}\hfill \end{array}$$

where $S_{\mathrm{Alg}1}^k$ , $S_{\mathrm{PSO}}^k$ , and $S_{\mathrm{OPT}}^k$ are the solutions obtained by Algorithm 1, PSO, and exhaustive enumeration algorithm of kth Monte Carlo simulation, respectively. From Table 2, it is seen that with the increase in power budget (from $120$ to $180$ ), Algorithm 1 achieves good approximation to its optimum value. The run-time of Algorithm 1 increases slowly. The run-time of the optimal algorithm increases very fast. From Table 3, it is seen that with the increase in the number of transmitters (from $16$ to $20$ ), Algorithm 1 achieves good approximation to its optimum value again. The run-time increase of Algorithm 1 is not obvious. Again, the run-time of the optimal algorithm increases fast. Table 4 shows that with the increase in targets (from $6$ to $10$ ), Algorithm 1 achieves good approximation to its optimum value once again. These results clearly demonstrate that Algorithm 1 is significantly more computationally efficient than the enumeration algorithm of which the computational cost increases exponentially with the increase in the number of network elements. Regarding the performance comparison between Algorithm 1 and PSO, it is seen that Algorithm 1 is a little bit worse (better) than PSO in terms of approximation factor (run-time). On the whole, Algorithm 1 and PSO have comparable performance. However, recall that the former has the advantage of theoretical performance guarantee, but the latter does not.

Table 2.

Statistical results of Algorithm 1 and PSO (Case $1$ , Scenario $2$ ).

$p_{\mathrm{total}}$	$A_{\mathrm{mean}}^{\mathrm{Alg}1}\backslash A_{\mathrm{mean}}^{\mathrm{PSO}}$	$T_{\mathrm{mean}}^{\mathrm{Alg}1}(s)\backslash T_{\mathrm{mean}}^{\mathrm{PSO}}(s)$	$T_{\mathrm{OPT}}(s)$
120	1.05\1.02	1.13\1.76	13.6
140	1.04\1.01	1.32\2.01	29.0
160	1.03\1.01	1.58\2.35	54.0
180	1.02\1.01	1.80\2.73	85.5

PSO: particle swarm optimization.

Table 3.

Statistical results of Algorithm 1 and PSO (Case 2, Scenario 2).

M	$A_{\mathrm{mean}}^{\mathrm{Alg}1}\backslash A_{\mathrm{mean}}^{\mathrm{PSO}}$	$T_{\mathrm{mean}}^{\mathrm{Alg}1}(s)\backslash T_{\mathrm{mean}}^{\mathrm{PSO}}(s)$	$T_{\mathrm{OPT}}(\mathrm{s})$
16	1.06\1.02	0.92\1.74	5.88
17	1.08\1.03	1.04\1.75	9.07
18	1.07\1.02	1.07\1.66	12.9
19	1.05\1.03	1.17\1.65	19.2
20	1.07\1.03	1.26\1.63	27.9

PSO: particle swarm optimization.

Table 4.

Statistical results of Algorithm 1 and PSO (Case 3, Scenario 2).

I	$A_{\mathrm{mean}}^{\mathrm{Alg}1}\backslash A_{\mathrm{mean}}^{\mathrm{PSO}}$	$T_{\mathrm{mean}}^{\mathrm{Alg}1}(s)\backslash T_{\mathrm{mean}}^{\mathrm{PSO}}(s)$	$T_{\mathrm{OPT}}(s)$
6	1.05\1.02	1.31\2.00	15.5
7	1.06\1.02	1.58\2.38	18.5
8	1.07\1.02	1.88\2.79	21.7
9	1.07\1.03	2.12\3.13	24.2
10	1.08\1.03	2.27\3.36	25.7

PSO: particle swarm optimization.

Conclusion

We have addressed the problem of transmitter subset selection subject to a power budget constraint in a distributed MIMO radar for minimizing the variance of target localization error. For a given transmitter/receiver configuration, we explicitly derived the CRLB, which is a function of target location, and depicts the minimum achievable localization error. Based on the CRLB, for a subset of transmitters, we defined a variance reduction function and used it as the optimization objective. We proved that the variance reduction function is a non-decreasing submodular function. By exploiting submodularity of the objective function, a balanced selection policy, which has performance guarantee and polynomial-time computational complexity, was adopted to solve the underlying problem.