# An Intelligent Prediction System for Mobile Source Localization Using Time
Delay Measurements

Hengnian Qi, Xiaoping Wu, and Naixue Xiong This work was supported by Zhejiang
Key R&D Plan 2017C03047 and Zhejiang Province Key Laboratory of Smart
Management & Application of Modern Agricultural Resources under Grant
2020E10017. _(Corresponding author: Xiaoping Wu)_ Hengnian Qi and Xiaoping Wu
are with the School of Information Engineering, Huzhou University, Huzhou
313000, China. (Email<EMAIL_ADDRESS>wxp@zafu.edu.cn).Naixue Xiong is with
the Department of Mathematics and Computer Science, Northeastern State
University, Tahlequah, OK, 74464, USA. (Email: xiongnaixue@gmail.com).

###### Abstract

In this paper, we introduce an intelligent prediction system for mobile source
localization in industrial Internet of things. The position and velocity of
mobile source are jointly predicted by using Time Delay (TD) measurements in
the intelligent system. To predict the position and velocity, the Relaxed
Semi-Definite Programming (RSDP) algorithm is firstly designed by dropping the
rank-one constraint. However, dropping the rank-one constraint leads to
produce a suboptimal solution. To improve the performance, we further put
forward a Penalty Function Semi-Definite Programming (PF-SDP) method to obtain
the rank-one solution of the optimization problem by introducing the penalty
terms. Then an Adaptive Penalty Function Semi-Definite Programming (APF-SDP)
algorithm is also proposed to avoid the excessive penalty by adaptively
choosing the penalty coefficient. We conduct experiments in both a simulation
environment and a real system to demonstrate the effectiveness of the proposed
method. The results have demonstrated that the proposed intelligent APF-SDP
algorithm outperforms the PF-SDP in terms of the position and velocity
estimation whether the noise level is large or not.

###### Index Terms:

mobile localization, semidefinite programming, time delay, rank-one solution,
intelligent prediction, industrial Internet of things

##  I Introduction

Nowadays, a popular trend in many application systems is the using of mobile
sources, such as Autonomous Vehicle (AV) and Unmanned Aerial Vehicle (UAV).
The using of mobile sources shows their great advantages for the flexible
mobility ability. Apparently, these mobile sources of AV and UAV can complete
some complicated tasks, such as location-based services, radar or sonar
navigation, target tracking, wireless network coverage and sensing
enhancement, data collection, etc [1, 2, 3, 4, 5, 6]. It is very crucial to
obtain the position parameter along with the velocity of mobile source when
the mobile source assists to complete these tasks. Due to the mobility of
mobile source, the position obtaining of mobile source is more complicated
compared with immobile source [7, 8, 9]. A common way to obtain the position
is to utilize the sensors with known positions. Besides the sensors, the
ranging information between source and sensor is also measured using some
technologies, such as Time Of Arrival (TOA) [10, 11], Time Difference Of
Arrival (TDOA) [12, 13], Received Signal Strength (RSS) [14, 15], Angle Of
Arrival (AOA) [16], and recently popular Time Delay (TD) [17, 18, 19]. Among
these technologies, TD is considered to be a simple and effective method to to
predict the position and velocity of the mobile source. [20].

Most of recent researches on the mobile source are focused on the obtaining of
motion parameters by using the motion sensor or Doppler shift measurements
[21, 22, 23]. In [24], the inertial navigation method is proposed to predict
the position of the mobile source using the motion sensors. When the motion
data is subjected to error, the filtering method is proposed to improve the
performance in [25, 26]. A linear algebraic method is put forward to predict
the motion parameters using the measurement of differential delay and Doppler
shift [27]. The direct acquisition of motion data depends on the hardware
devices, and the measurements may exist errors. Therefore, recent researches
are concentrated on how to utilize various ranging methods to directly predict
the motion parameters of mobile sources [15].

In this paper we propose an intelligent practice prediction method of motion
parameters including the position and velocity of mobile source by only using
the time delay (TD) measurements. The proposed method does not require any
motion sensors or Doppler shift measurements. To predict the position and
velocity of the mobile source, the optimization model is firstly built by
representing the measurement equation into a linear form. Then a Relaxed Semi-
Definite Programming (RSDP) form is designed by relaxing the non-convex model
into a convex problem. However, the solution to the RSDP problem is suboptimal
due to the relaxation of rank-one constraint. To obtain the rank-one solution,
we also propose the Penalty Function (PF) method by introducing the penalty
terms. Then the APF-SDP algorithm is designed to obtain the optimal
performance of the prediction problem by adaptively choosing the penalty
coefficient.

We also deploy a simulation environment and a real experiment system to fully
demonstrate the effectiveness of the proposed algorithms. The real system is
composed of a mobile car equipped with Ultrasonic module (UM) and nine
sensors. Besides the UM, motion sensors are also equipped to the mobile car
and used to measure the true position and velocity of the mobile source. The
experimental results show that APF-SDP significantly outperforms the PF-SDP in
terms of the position and velocity accuracy. 20% of the position error is
larger than 0.05 m for RSDP and 0.04 m for PF-SDP. However, it is reduced to
near 0.01 m for the APF-SDP. By adaptively choosing the penalty coefficient,
the performance of APF-SDP is still sufficiently close to the expected Cramér-
Rao Lower Bound (CRLB) of the prediction problem even if the number of sensors
is varied from 5 to 9. It confirms the advantage of APF-SDP by achieving the
rank-one constraint.

Our proposed TD-based localization model does not require any motion sensors
or Doppler shift measurements, so it differs from the existing mobile source
localization problem in [28, 24]. Moreover, to obtain the well performance of
the estimator, we propose the rank-one solution method using the penalty
function. Our rank-one solution of AFP-SDP method is globally convergent by
adaptively choosing the penalty coefficient and essentially different from the
existing rank-one solution methods proposed in [29, 30, 31]. The contributions
of this work are summarized as follows,

  1. 1.

To predict the position along with the velocity of the mobile source, dropping
the rank-one constraint produces the RSDP problem. To design a tighter convex
model, we propose the penalty function method to obtain a rank-one solution by
introducing the penalty terms.

  2. 2.

To avoid the excessive penalty, we also put forward the APF-SDP algorithm by
adaptively choosing the penalty coefficient. We have theoretically proven that
the APF-SDP can provide optimal performance of the prediction problem by
achieving the rank-one constraint.

  3. 3.

We have also developed an intelligent practice prediction system to
demonstrate our proposed algorithms using TD measurements. The results from
both the simulated and real experiments show that the APF-SDP algorithm
provides better performance than the PF-SDP whether the noise level is large
or not.

The rest of this paper is structured as follows. Related works are introduced
in Section II. Section III presents the system model and definitions. Section
IV in detail describes our proposed intelligent prediction system. In Section
V, the computational complexity of these proposed algorithms is derived. The
numerical simulations and real experiments are analyzed in Section VI. The
conclusions and future work are presented in Section VII. This paper contains
a number of symbols. Following the convention, the matrix is represented by
bold case letter. If the matrix is denoted by $(\ast)$, $(\ast)^{-1}$ and
$(\ast)^{T}$ indicate the matrix inverse and transpose operator, respectively.
$\lVert\ast\rVert$ denotes $\mathit{\ell}_{2}$ norm. $\bm{A}_{i,j}$ is the
$(i,j)$th element of matrix $\bm{A}$. $\bm{0}_{m}$ represents all-zero vector
with the length $m$, and $\bm{I}_{p}$ and $\bm{0}_{m\times n}$ are $m\times m$
identity and $m\times n$ zero matrices. Tr($\bm{A}$) and rank($\bm{A}$) stand
for the trace and rank of $\bm{A}$, respectively. $\bm{A}\succeq 0$ indicates
that $\bm{A}$ is positive semidefinite.

##  II Related Works

Based on the ranging information, various algorithms are proposed to predict
the position and the velocity of mobile source. The popular algorithms include
Maximum Likelihood Estimator (MLE) [32, 33], alternating direction method of
multiplier (ADMM) method [34, 35], linear estimator [36, 37], and convex
method [38, 39]. The numerical solution of MLE requires a very good initial
solution to guarantee its global convergence. The ADMM method provides the
optimal estimate of source position by converting the unconstrained nonlinear
problem into an equivalent constrained form. Due to the nonlinear nature of
the optimization model, the MLE or ADMM method will be trapped in local
optimum. To avoid the problem, the linear estimator directly represents the
unknown variables as an algebraic analytic form by converting the nonlinear
model into a linear problem. However, the constrained relationship among the
variables is difficult to be exploited in the linear estimator, so the
performance is not well enough. Convex method does not depend on the
initialization for its global convergence and gradually becomes a popular
method for the source position prediction problem. The convex methods can be
achieved by Semi-Definite Programming (SDP) [40, 41, 42] and Second Order Cone
Programming (SOCP) [43] which can be efficiently solved by using existing
algorithms such as interior point methods [44].

A common method to obtain an SDP problem is to relax the non-convex
optimization model into a convex form by dropping the rank-one constraint. The
rank-one relaxation may lead to produce a suboptimal solution that is not the
optimal solution of the original optimization problem. To obtain the rank-one
solution of the convex SDP problem, many mathematical methods are proposed to
deal with the troublesome problem [45, 46, 47, 48]. To solve the low-rank SDP
problems, the factorization method is introduced by obtaining a reformulation
of the original SDP problem in [49]. A modified interior point method is
proposed to solve low-rank SDP problem in [50]. Although these mathematical
tools deal with the low-rank SDP problems, the solutions are also not
guaranteed to be global convergence for their nonlinear or non-convex nature.

To obtain the rank-one solution with global convergence, recently the two-
stage method is proposed by refining the initial solution of the relaxed SDP
problem. In [29], a rank-reduction method is designed by solving an
incremental matrix of the solution when an initial rank-maximum solution has
been obtained with the relaxed SDP problem. The relaxed SDP problem provides a
suboptimal solution, so a feasible method to solve the rank-one solution is to
continuously tighten the relaxed problem. In [30], a set of SOCP constraints
is added to tighten the convex model and find a rank-one solution of the
original SDP problem. Above rank-one solutions strongly depend on the initial
solutions of the relaxed SDP problem. Therefore, if the initial solutions are
not accurate enough, the rank-one solution may fail. In [31], the Penalty
Function (PF) method is firstly proposed to obtain the rank-one solution of
the SDP problem by introducing the penalty term. The rank-one solution of PF
method provides a global optimum solution by choosing an appropriate penalty
coefficient. However, too large penalty coefficient will lead to the
occurrence of excessive penalty, which badly affects the performance of the
solution.

In this paper, we mainly propose an intelligent practice system, in which the
motion parameters including the position and the velocity of the mobile source
are predicted by only using the TD measurements. The proposed method does not
require any motion sensors or Doppler shift measurements, and is thus
applicable to track the mobile source, such as the AV or UAV.

##  III System Model and Definitions

Assuming that $M$ sensors with known positions are placed in a $p$-dimensional
($p=2$ or $p=3$) scenario. The known positions of the sensors are denoted by
$\bm{s}_{i}\in\mathbb{R}^{p}$, $i=1,2,\ldots,M$. In the same scenario, a
mobile source starts from an initial position $\bm{u}\in\mathbb{R}^{p}$ at a
constant velocity $\bm{v}\in\mathbb{R}^{p}$. Fig. 1 shows the diagram of TD
measurement between mobile source and sensor. A mobile source transmits a
signal at the initial position. Then the signal is received by the sensors and
immediately reflected to the mobile source, thus forming a time delay that is
formulated as:

$t_{i}^{o}=\frac{1}{c}(\|\bm{u}-\bm{s}_{i}\|+\|\bm{u}_{i}-\bm{s}_{i}\|),$ (1)

where $t_{i}^{o}$ is the true time delay, $c$ is the propagation speed of the
signal, $\bm{u}_{i}$ is the instantaneous position of the mobile source when
the signal from the $i$th sensor is received by the mobile source. Apparently,
we have

$\bm{u}_{i}=\bm{u}+\bm{v}t_{i}^{o}.$ (2)

Substituting (2) into (1) and multiplying $c$ on the both sides produce

$ct_{i}^{o}=\|\bm{u}-\bm{s}_{i}\|+\|\bm{u}+\bm{v}t_{i}^{o}-\bm{s}_{i}\|.$ (3)

Moving the term $\|\bm{u}-\bm{s}_{i}\|$ to the left side of (3) and squaring
on the both sides give

$2(\bm{u}-\bm{s}_{i})^{T}\bm{v}+2cd_{i}+(\bm{v}^{T}\bm{v}-c^{2})t_{i}^{o}=0,$
(4)

where $\bm{d}_{i}=\|\bm{u}-\bm{s}_{i}\|$. Therefore, the true time delay is
written as:

$t_{i}^{o}=\frac{2(\bm{u}-\bm{s}_{i})^{T}\bm{v}+2cd_{i}}{c^{2}-\bm{v}^{T}\bm{v}}.$
(5)

The measured time delay $t_{i}$ is always subject to error and given by

$t_{i}=t_{i}^{o}+n_{i},$ (6)

where $i=1,2,\ldots,M$, $n_{i}$ is the noise that constructs a vector form
$\bm{n}=[n_{1},n_{2},\ldots,n_{M}]^{T}$. The noise vector $\bm{n}$ is always
zero mean Gaussian distribution with covariance matrix $\bm{Q}_{n}$.

Figure 1: A diagram of TD measurements between sensor and mobile source.

The goal of our system model is to predict the initial position and velocity
of the mobile source using the noisy measurement $t_{i}$. In the following,
the relaxed SDP (RSDP) method is firstly proposed to predict the position and
velocity of the mobile source by dropping the rank-one constraint. The drop of
the rank-one constraint is considered as a relaxation method that leads to
produce an SDP solution with the rank higher than 1. So the performance of
RSDP is suboptimal. To improve the performance, the penalty function method is
proposed to obtain the rank-one solution by introducing the penalty terms.

##  IV Intelligent Prediction System

In this section, we in detail introduce the intelligent prediction system for
predicting the position and velocity of mobile source using TD measurement.
The convex SDP problem shows the advantages for its global convergence and can
be solved efficiently using interior-point algorithms. In our proposed system,
we mainly concentrate on the convex SDP solution to the mobile source
localization problem. To obtain a convex SDP form, a general approach is to
relax the non-convex optimization problem into a convex form, then the
position and the velocity of mobile source are predicted by extracting from
the SDP solution.

###  IV-A RSDP Algorithm

Substituting (6) into (4) yields the following expression:

$-2\bm{s}_{i}^{T}\bm{v}+2\bm{u}^{T}\bm{v}+t_{i}\bm{v}^{T}\bm{v}+2cd_{i}-c^{2}t_{i}=\varepsilon_{i},$
(7)

where $\varepsilon_{i}=(\bm{v}^{T}\bm{v}-c^{2})n_{i}$, $i=1,2,\ldots,M$. To
establish the optimization model, a new unknown vector $\bm{x}$ is defined by
$\bm{x}=[\bm{u}^{T},\bm{v}^{T},\bm{u}^{T}\bm{v},\bm{v}^{T}\bm{v},\bm{d}^{T},1]^{T}\in\mathbb{R}^{M+2p+3}$,
$\bm{d}=[d_{1},d_{2},\ldots,d_{M}]^{T}$. Then (7) is also represented by a
linear matrix form:

$\bm{Gx}=\bm{\varepsilon},$ (8)

where
$\bm{\varepsilon}=[\varepsilon_{1},\varepsilon_{2},\ldots,\varepsilon_{M}]^{T}$.
According to the expression of (7), $\bm{\varepsilon}\in\mathbb{R}^{M}$ and
$\bm{G}\in\mathbb{R}^{M\times(M+2p+3)}$ are also defined by

$\displaystyle\bm{\varepsilon}=\bm{B}\bm{n},$ (9a)
$\displaystyle\bm{B}=(\bm{v}^{T}\bm{v}-c^{2})\bm{I}_{M},$ (9b)
$\displaystyle\bm{G}=[\bm{g}_{1}^{T},\bm{g}_{2}^{T},\ldots,\bm{g}_{M}^{T}]^{T},$
(9c)
$\displaystyle\bm{g}_{i}=[\bm{0}_{p}^{T},-2\bm{s}_{i}^{T},2,t_{i},\bm{0}_{i-1}^{T},2c,\bm{0}_{M-i}^{T},-t_{i}c^{2}]^{T},$
(9d)

Therefore, the weighted least square (WLS) solution to the prediction problem
is formulated as

$\displaystyle\min_{\bm{x}}\quad(\bm{Gx})^{T}\bm{Q}^{-1}(\bm{Gx})$
$\displaystyle\mathrm{s.t.}\quad\|\bm{x}_{p+1:2p}\|^{2}=\bm{x}_{2p+2},$ (10a)
$\displaystyle\quad\quad\|\bm{x}_{1:p}-\bm{s}_{i}\|=\bm{x}_{2p+2+i},i=1,2\ldots,M,$
(10b) $\displaystyle\quad\quad\bm{x}_{1:p}^{T}\bm{x}_{p+1:2p}=\bm{x}_{2p+1},$
(10c) $\displaystyle\quad\quad\bm{x}_{M+2p+3}=1,$ (10d)

where $\bm{Q}\in\mathbb{R}^{M\times M}$ is the covariance matrix with respect
to $\bm{\varepsilon}$, and $\bm{Q}$ is further obtained by

$\displaystyle\bm{Q}=\bm{B}\bm{Q}_{n}\bm{B}.$ (11)

The equality condition (10a) is also rewritten as

$\|\bm{D}_{0}\bm{x}\|^{2}=\bm{x}_{2p+2},$ (12)

where $\bm{D}_{0}=[\bm{0}_{p},\bm{I}_{p},\bm{0}_{p\times(M+3)}]$. The equality
condition (10b) is also given by

$\|\bm{D}_{i}\bm{x}\|=\bm{x}_{2p+2+i},$ (13)

where $\bm{D}_{i}=[\bm{I}_{p},\bm{0}_{p\times(M+p+2)},\bm{s}_{i}]$,
$i=1,2,\ldots,M$. When the equality conditions (10a) and (10b) are equivalent
to the conditions (12) and (13), (10) is rewritten as

$\displaystyle\min_{\bm{x}}\quad(\bm{Gx})^{T}\bm{Q}^{-1}(\bm{Gx})$
$\displaystyle\mathrm{s.t.}\quad\|\bm{D}_{0}\bm{x}\|^{2}=\bm{x}_{2p+2},$
$\displaystyle\quad\quad\|\bm{D}_{i}\bm{x}\|=\bm{u}_{2p+2+i},i=1,2\ldots,M,$
$\displaystyle\quad\quad\bm{x}_{1:p}^{T}\bm{x}_{p+1:2p}=\bm{x}_{2p+1},$
$\displaystyle\quad\quad\bm{x}_{M+2p+3}=1.$ (14)

To further express (IV-A) as a convex form, a new unknown matrix is defined by
$\bm{X}=\bm{x}\bm{x}^{T}$. It is obviously shown that
$\mathrm{rank(\bm{X})}=1$. Thus, (IV-A) is given by

$\displaystyle\min_{\bm{X}}\quad\mathrm{Tr}(\bm{C}\bm{X})$
$\displaystyle\mathrm{s.t.}\quad\mathrm{Tr}(\bm{D}_{0}^{T}\bm{D}_{0}\bm{X})=\bm{X}_{2p+2,2p+2},$
(15a)
$\displaystyle\quad\quad\mathrm{Tr}(\bm{D}_{i}^{T}\bm{D}_{i}\bm{X})=\bm{X}_{2p+2+i,2p+2+i},i=1,2\ldots,M,$
(15b)
$\displaystyle\quad\quad\sum\nolimits_{i=1}^{p}\bm{X}_{i,p+i}=\bm{X}_{2p+1,2p+1},$
(15c) $\displaystyle\quad\quad\bm{X}_{M+2p+3,M+2p+3}=1,\bm{X}\succeq 0,$ (15d)
$\displaystyle\quad\quad\mathrm{rank}(\bm{X})=1,$ (15e)

where $\bm{C}=\bm{G}^{T}\bm{Q}^{-1}\bm{G}$. Due to the rank-one constraint of
(15e), the expression of (15) is non-convex. However, dropping the rank-one
condition of $\bm{X}$ yields a convex SDP form,

$\displaystyle\min_{\bm{X}}\quad\mathrm{Tr}(\bm{C}\bm{X})$
$\displaystyle\mathrm{s.t.}\quad\bm{A}_{i}\bm{X}=b_{i},i=1,\ldots,M+3,$ (16)

where $b_{i}=[\bm{0}_{1\times(M+2)},1]^{T}$, the sparse matrices $\bm{A}_{i}$
can be easily obtained by the expression of (15), $i=1,\ldots,M+3$. It is
obviously shown that the problem depicted by (IV-A) is a typical SDP problem
that can be effectively solved using a convex optimization package such as
SEDUMI and SDPT3. Unfortunately, dropping the rank-one constraint relaxes the
problem depicted by (15) and leads to produce a suboptimal SDP solution with a
rank higher than 1, which implies that the solution of (IV-A) may not be the
optimal solution of the original problem depicted by (15). To obtain the
optimal performance, we propose the penalty function method to meet the rank-
one constraint.

###  IV-B PF-SDP Algorithm

The convex problem depicted by (16) provides a relaxed SDP solution to the
mobile source localization problem. However, the relaxed SDP solution of (16)
is suboptimal due to the drop of rank-one constraint. To obtain the optimal
performance, we attempt to find a rank-one solution by introducing the penalty
terms. To design the PF-SDP method, we firstly prove the conclusion of Lemma
1.

######  Lemma 1.

For a given solution $\bm{X}\in\mathbb{R}^{N\times N}$ to the problem depicted
by (16), if $\bm{X}_{i,i}=\bm{X}_{i,N}^{2}$ $(N=M+2p+3$, $i=1,2,\ldots,N-1)$,
then $\mathrm{rank}(\bm{X})=1$.

###### Proof.

If $\bm{X}$ is a positive semidefinite solution of (16), its principal
$2\times 2$ submatrix is also a positive semidefinite matrix. It is noted that
$\bm{X}_{N,N}=1$. By selecting the $(i,i)$ and $(N,N)$ entry of $\bm{X}$ as
the principal diagonal elements, a new positive semidefinite submatrix is
constructed and given by

$\left[\begin{array}[]{ccc}\bm{X}_{i,i}&\bm{X}_{i,N}\\\ \bm{X}_{N,i}&1\\\
\end{array}\right]\succeq 0,$ (17)

where $i=1,2,\ldots,N-1$. Since $\bm{X}_{i,N}=\bm{X}_{N,i}$, we have
$\bm{X}_{i,i}\geq\bm{X}_{i,N}^{2}$. When the equality condition is satisfied,
it is also modified as

$\sqrt{\bm{X}_{i,i}}=\bm{X}_{i,N}.$ (18)

Every principal $3\times 3$ submatrix of positive semidefinite matrix $\bm{X}$
is also positive semidefinite. Similarly, selecting the $(i,i)$, $(j,j)$, and
$(N,N)$ entry of $\bm{X}$ as the principal diagonal elements and using the
equality (18), we also construct a $3\times 3$ submatrix that is given by

$\left[\begin{array}[]{cccc}\bm{X}_{i,i}&\bm{X}_{i,j}&\sqrt{\bm{X}_{i,i}}\\\
\bm{X}_{j,i}&\bm{X}_{j,j}&\sqrt{\bm{X}_{j,j}}\\\
\sqrt{\bm{X}_{i,i}}&\sqrt{\bm{X}_{j,j}}&1\end{array}\right]\succeq 0,$ (19)

where $i,j=1,2,\ldots,N-1$, $i<j$. (18) is also equivalent to the following
expression:

$\left[\begin{array}[]{ccc}0&\bm{X}_{i,j}-\sqrt{\bm{X}_{i,i}\bm{X}_{j,j}}\\\
\bm{X}_{i,j}-\sqrt{\bm{X}_{i,i}\bm{X}_{j,j}}&0\\\ \end{array}\right]\succeq
0,$ (20)

where $i,j=1,2,\ldots,N-1$, $i<j$. The equation (20) holds if and only if
$\bm{X}_{i,j}=\sqrt{\bm{X}_{i,i}\bm{X}_{j,j}}$ for all $i,j=1,2,\ldots,N-1$,
$i<j$. Therefore, the positive semidefinite matrix $\bm{X}$ is also rewritten
as

$\bm{X}=\left[\begin{array}[]{cccccccc}\bm{X}_{1,1}&\sqrt{\bm{X}_{1,1}\bm{X}_{2,2}}&\ldots&\sqrt{\bm{X}_{1,1}}\\\
\sqrt{\bm{X}_{1,1}\bm{X}_{2,2}}&\bm{X}_{2,2}&\ldots&\sqrt{\bm{X}_{2,2}}\\\
\vdots&\vdots&\vdots&\vdots\\\
\sqrt{\bm{X}_{1,1}}&\sqrt{\bm{X}_{2,2}}&\ldots&1\end{array}\right].\\\ $ (21)

It is obviously shown that $\mathrm{rank}(\bm{X})=1$. ∎

To achieve the equality (18), a feasible method is to restrict the magnitude
of $\bm{X}_{i,i}$ by introducing the penalty terms. Therefore, a new PF-based
SDP form is given by

$\displaystyle\min_{\bm{X}}\quad\mathrm{Tr}(\bm{C}\bm{X})+\eta\sum\nolimits_{i=1}^{N-1}\bm{X}_{i,i}$
$\displaystyle\mathrm{s.t.}\quad\bm{A}_{i}\bm{X}=b_{i},i=1,\ldots,M+3,$ (22)

where $\eta$ is a strictly positive constant and called as penalty
coefficient. When $\eta$ is gradually increased, $\bm{X}_{i,i}$ will
sufficiently close to $\bm{X}_{i,N}^{2}$. Eventually, large enough $\eta$ will
lead that $\mathrm{rank}(\bm{X})=1$.

When $\eta$ is chosen to be large enough, the SDP problem depicted by (IV-B)
provides a rank-one solution of the positive semidefinite $\bm{X}$. However, a
concerned issue is that whether the rank-one solution is the original solution
of (10) or not. It is directly shown that the problem depicted by (10) has
only $M+2p+3$ variables. Since there are $M+3$ constraints among the variables
in defined vector $\bm{x}$, so the problem has only $2p$ independent
variables. It is also observed that the dimension of $\bm{X}$ in the problem
depicted by (IV-B) is also $M+2p+3$ when $\mathrm{rank}(\bm{X})=1$. Due to the
similar constraints, the problem depicted by (IV-B) has also only $2p$
independent variables. Therefore, the rank-one solution of the SDP problem
depicted by (IV-B) is also the original solution of (10) since the convex SDP
form of (IV-B) is derived from the non-convex problem depicted by (10).

Figure 2: Localization error and optimal objective
$\mathrm{Tr}(\bm{C}\bm{X}^{*})$ with a test.

An appropriate penalty coefficient can ensure the effective working of penalty
terms. Although large enough $\eta$ will achieve the rank of the solved
$\bm{X}$ to be one, too large $\eta$ may lead to a risk of excessive penalty,
which will badly affect the predicted result. Hence, the choosing of $\eta$ is
crucial to obtain the well performance of the prediction problem. To choose an
appropriate penalty coefficient, we propose an adaptive penalty function-based
SDP (AFP-SDP) algorithm to obtain a rank-one solution for the SDP optimization
problem in the following.

###  IV-C APF-SDP Algorithm

Too large $\eta$ increases the proportion of penalty terms to the total cost
function of (22), so it will weaken the original objective
$\mathrm{Tr}(\bm{C}\bm{X})$. To avoid the chosen $\eta$ to be too large,
$\eta$ starts from a small value and increases gradually, until the rank-one
condition is satisfied. To clearly observe the effect of the penalty
coefficient in the problem depicted by (22), we conducted a random test and
the results are plotted in Fig.1, where $\bm{X}^{*}$ and
$\mathrm{Tr}(\bm{C}\bm{X}^{*})$ represent the optimal solution and objective,
respectively. In the test, $\mathrm{log}\eta$ is gradually increased from -6
to 3 (i.e. $\eta$ is gradually increased from $10^{-6}$ to $10^{3}$). As can
be seen that the penalty terms almost do not work when $\mathrm{log}\eta$ is
increased from -6 to -4. The localization error is gradually reduced for the
working of penalty terms when $\mathrm{log}\eta$ is increased from -4 to -2.
The optimal performance is achieved when $\mathrm{log}\eta$ is set to the
range of (-2, 0). However, if $\mathrm{log}\eta$ is continuously to be
increased, the performance will become worse. The localization error will
sharply increased due to the occurrence of excessive penalty when
$\mathrm{log}\eta$ is larger than 1. The occurrence of excessive penalty badly
affects the localization result, so it is crucial to detect whether the
penalty is excessive or not. It can be seen from Fig.2 that the optimal
objective $\mathrm{Tr}(\bm{C}\bm{X}^{*})$ is simultaneously increased with the
localization error when the excessive penalty occurs. Therefore, the optimal
objective $\mathrm{Tr}(\bm{C}\bm{X}^{*})$ can be used to detect the occurrence
of excessive penalty.

######  Lemma 2.

If $\bm{X}^{*}$ is an optimal rank-one solution of SDP problem depicted by
(16), then
$\mathrm{Tr}(\bm{C}\bm{X}^{*})\thicksim\sum_{i=1}^{M}\lambda_{i}\chi^{2}(1)$.

###### Proof.

If $\bm{X}^{*}$ is an optimal rank-one solution of SDP problem depicted by
(16), a new vector $\bm{x}^{*}$ is defined by
$\bm{X}^{*}=\bm{x}^{*}\bm{x}^{*T}$, where $\bm{x}^{*}\in\mathbb{R}^{M+2p+3}$.
Since the problem is derived from (10), $\bm{x}^{*}$ is also an optimal
solution of problem depicted by (10). Then it yields

$\displaystyle\mathrm{Tr}(\bm{C}\bm{X}^{*})=\bm{x}^{*T}\bm{C}\bm{x}^{*}.$ (23)

The prediction error of $\bm{x}^{*}$ is denoted by $\triangle\bm{x}^{*}$, then
$\bm{x}^{*}+\triangle\bm{x}^{*}=\bm{x}^{o}$, where $\bm{x}^{o}$ represents the
true value of the defined $\bm{x}$.

The constrained optimization problem depicted by (10) has only $2p$
independent variables that constructs a vector
$\bm{z}=[\bm{u}^{T},\bm{v}^{T}]^{T}$. If the prediction error of $\bm{z}$ is
denoted as $\triangle\bm{z}$, we can obtain

$\displaystyle\frac{\partial(\bm{G}\bm{x})}{\partial\bm{x}}\frac{\partial\bm{x}}{\partial\bm{z}}\triangle\bm{z}=\bm{\varepsilon},$
(24)

where $\frac{\partial(\bm{G}\bm{x})}{\partial\bm{x}}=\bm{G}$,
$\frac{\partial\bm{x}}{\partial\bm{z}}$ is denoted as $\bm{H}$ and given by

$\displaystyle\bm{H}=\left[\begin{array}[]{ccccccc}\bm{I}_{p}&\bm{0}_{p\times
p}&\bm{v}&\bm{0}_{p}&\frac{\partial\bm{d}}{\partial\bm{u}}&\bm{0}_{p}\\\
\bm{0}_{p\times p}&\bm{I}_{p}&\bm{u}&2\bm{v}&\bm{0}_{p\times M}&\bm{0}_{p}\\\
\end{array}\right]^{T},$ (27)

where
$\frac{\partial\bm{d}}{\partial\bm{u}}=[\frac{\bm{u}-\bm{s}_{1}}{\|\bm{u}\|},\frac{\bm{u}-\bm{s}_{2}}{\|\bm{u}\|},\ldots,\frac{\bm{u}-\bm{s}_{M}}{\|\bm{u}\|}]$.
Therefore, (24) is rewritten as

$\displaystyle\bm{P}\triangle\bm{z}=\bm{\varepsilon},$ (28)

where $\bm{P}=\bm{G}\bm{H}$. Since $\bm{\varepsilon}=\bm{B}\bm{n}$, the WLS
solution to (28) is

$\displaystyle\triangle\bm{z}=(\bm{P}^{T}\bm{W}\bm{P})^{-1}\bm{P}^{T}\bm{W}\bm{B}\bm{n},$
(29)

where $\bm{W}=\bm{Q}^{-1}$. Since $\triangle\bm{x}^{*}=\bm{H}\triangle\bm{z}$,
the prediction error $\triangle\bm{x}^{*}$ is obtained by

$\displaystyle\triangle\bm{x}^{*}=\bm{E}\bm{n},$ (30)

where $\bm{E}=\bm{H}(\bm{P}^{T}\bm{W}\bm{P})^{-1}\bm{P}^{T}\bm{W}\bm{B}$.
Since $\bm{x}^{*}=\bm{x}^{o}-\triangle\bm{x}^{*}$, the optimal objective is
approximately given by

$\displaystyle\bm{x}^{*T}\bm{C}\bm{x}^{*}\approx\bm{x}^{*T}\bm{C}\bm{x}^{o}.$
(31)

Substituting $\bm{x}^{*}=\bm{x}^{o}-\triangle\bm{x}^{*}$ into the right side
of (31) yields

$\displaystyle\bm{x}^{*T}\bm{C}\bm{x}^{*}=\bm{x}^{oT}\bm{C}\bm{x}^{o}-\bm{x}^{oT}\bm{C}\triangle\bm{x}^{*},$
(32)

where $\bm{C}$ is interfered by the noise. Let
$\bm{C}=\bm{C}^{o}+\triangle\bm{C}$, $\bm{C}^{o}$ and $\triangle\bm{C}$ denote
the true value and the error term caused by noise, respectively. Similarly, we
also define $\bm{G}=\bm{G}^{o}+\triangle\bm{G}$ in which $\bm{G}^{o}$ and
$\triangle\bm{G}$ represent the true value and the error term,

$\displaystyle\triangle\bm{G}=[\triangle\bm{g}_{1}^{T},\triangle\bm{g}_{2}^{T},\dots,\triangle\bm{g}_{M}^{T}]^{T}$
(33a)
$\displaystyle\triangle\bm{g}_{i}=[\bm{0}_{2p+1}^{T},n_{i},\bm{0}_{M}^{T},-c^{2}n_{i}]^{T}.$
(33b)

According to the definition of $\bm{C}=\bm{G}^{T}\bm{W}\bm{G}$, we can obtain
the true value $\bm{C}^{o}$ and the error term $\triangle\bm{C}$,

$\displaystyle\bm{C}^{o}=\bm{G}^{oT}\bm{W}\bm{G}^{o}$ (34a)
$\displaystyle\triangle\bm{C}=\bm{G}^{oT}\bm{W}\triangle\bm{G}+\triangle\bm{G}^{T}\bm{W}\bm{G}^{o}+\triangle\bm{G}^{T}\bm{W}\triangle\bm{G}.$
(34b)

Since $\bm{G}^{o}\bm{x}^{o}=\bm{0}_{M}$, the first term of the right side in
(32) is also further written as

$\displaystyle\bm{x}^{oT}\bm{C}\bm{x}^{o}\approx(\triangle\bm{G}\bm{x}^{o})^{T}\bm{W}(\triangle\bm{G}\bm{x}^{o}),$
(35)

where $\triangle\bm{G}\bm{x}^{o}$ is also given by

$\displaystyle\triangle\bm{G}\bm{x}^{o}=\bm{F}\bm{n}$ (36a)
$\displaystyle\bm{F}=\mathrm{diag}(\bm{T}\bm{x}^{o})$ (36b)
$\displaystyle\bm{T}=[\bm{t}_{1}^{T},\bm{t}_{2}^{T},\ldots,\bm{t}_{M}^{T}]^{T}$
(36c)
$\displaystyle\bm{t}_{i}=[\bm{0}_{2p+1}^{T},1,\bm{0}_{M}^{T},-c^{2}]^{T}.$
(36d)

Therefore, (35) is also rewritten as

$\displaystyle\bm{x}^{oT}\bm{C}\bm{x}^{o}=\bm{n}^{T}\bm{F}^{T}\bm{W}\bm{F}\bm{n}.$
(37)

Similarly, the second term of the right side in (32) can also rewritten as

$\displaystyle\bm{x}^{oT}\bm{C}\triangle\bm{x}^{*}=\bm{x}^{oT}(\triangle\bm{G}^{T}\bm{W}\bm{G}^{o}+\triangle\bm{G}^{T}\bm{W}\triangle\bm{G})\triangle\bm{x}^{*}.$
(38)

Neglecting the high order term of (38) also yields

$\displaystyle\bm{x}^{oT}\bm{C}\triangle\bm{x}^{*}\approx(\triangle\bm{G}\bm{x}^{o})^{T}\bm{W}\bm{G}^{o}\triangle\bm{x}^{*}.$
(39)

Using the expression of (28) and (34a), we can rewrite (39) as

$\displaystyle\bm{x}^{oT}\bm{C}\triangle\bm{x}^{*}=\bm{n}^{T}\bm{F}^{T}\bm{W}\bm{G}^{o}\bm{E}\bm{n}.$
(40)

According to the expressions (37) and (40), (32) is also rewritten as

$\displaystyle\bm{x}^{*T}\bm{C}\bm{x}^{*}\approx\bm{n}^{T}\bm{F}^{T}\bm{W}(\bm{F}-\bm{G}\bm{E})\bm{n},$
(41)

where $\bm{G}^{o}$ is approximately equal to $\bm{G}$. Let
$\bm{L}=\bm{F}^{T}\bm{W}(\bm{F}-\bm{G}\bm{E})$. Therefore,
$\bm{x}^{*T}\bm{C}\bm{x}^{*}$ conforms to the chi-square distribution
$\sum_{i=1}^{M}\lambda_{i}\chi^{2}(1)$, which is also denoted by

$\displaystyle\bm{x}^{*T}\bm{C}\bm{x}^{*}\thicksim\sum_{i=1}^{M}\lambda_{i}\chi^{2}(1),$
(42)

where $\lambda_{i}$ is the eigenvalue of the matrix $\bm{L}\bm{Q}_{n}$,
$i=1,2,\ldots,M$. ∎

When $\eta$ is chosen to be too small, the rank of $\bm{X}^{*}$ may be larger
than 1. The increasing of $\eta$ can ensure that the rank of $\bm{X}^{*}$ is
gradually close to be one. To evaluate the rank of $\bm{X}^{*}$, we firstly
propose an eigenvalue method according to the following threshold value,

$\displaystyle\tau=\frac{\widetilde{\lambda}_{N-1}}{\widetilde{\lambda}_{N}},$
(43)

where $\widetilde{\lambda}_{N-1}$ and $\widetilde{\lambda}_{N}$ are the
$N-1$th and the $N$th eigenvalue of $\bm{X}^{*}$, $N=M+2p+3$. It is obviously
shown that the rank of $\bm{X}^{*}$ tends to be one when $\tau$ is
sufficiently close to zero.

To ensure the well performance of the PF-SDP, it is crucial to choose an
optimal penalty coefficient which depends on various factors including the
positions of sensors and mobile source, the number of sensors, and the noise
level. Therefore, the optimal $\eta$ is not invariable under the different
situations. In Algorithm 1, we propose an Adaptive Penalty Function-based SDP
(APF-SDP) algorithm by adaptively choosing an optimal penalty coefficient.

In Algorithm 1, two conditions is required to be judged. One is the condition
$\tau<\delta$ ($\delta$ is a constant that is sufficiently close to zero),
which is considered as rank-one condition and used to judge whether the rank
of $\bm{X}^{*}$ is one or not. The other is the condition
$\bm{C}\bm{X}<\epsilon$ ($\epsilon$ is also a threshold value determined by
the chi-square distribution of $\bm{C}\bm{X}^{*}$), which is used to detect
whether the penalty is excessive or not. $\eta$ starts from a small value and
increases gradually with $\eta=\gamma\eta$ ($\gamma>1$), where $\gamma$ is
called as step length. When $\tau<\delta$ is satisfied, it is considered to
meet the rank-one condition. Then using the condition $\bm{C}\bm{X}<\epsilon$,
we further detect whether the penalty is excessive or not. When meeting the
condition $\bm{C}\bm{X}<\epsilon$, we accept the solution that is considered
as a final solution $\bm{X}^{*}$. If not, we relax the rank-one condition with
$\delta=\alpha\delta$ ($\alpha>1$) and resolve the problem depicted by (22)
until two conditions hold.

Input: $\bm{A}$, $\bm{b}$, $\bm{C}$, $\eta$, $\alpha>1$, $\gamma>1$,
$\delta>0$, $\epsilon>0$

Output: $\bm{X}^{*}$

1 while  _1_ do

2 Using (17) to calculate the optimal $\bm{X}$;

3 if  _$\tau <\delta$_ then

4 if  _$\mathrm{Tr}(\bm{C}\bm{X}) <\epsilon$_ then

5 $\bm{X}^{*}=\bm{X}$;

6 break;

7 else

8 $\delta=\alpha\delta$;

9 goto 2;

10 end if

11

12 else

13 $\eta=\gamma\eta$;

14 end if

15

16 end while

return $\bm{X}^{*}$

Algorithm 1 APF-SDP Algorithm TABLE I: Positions of ten sensors in simulations(km) $x$ | 0.145 | -0.020 | -0.207 | 0.400 | -0.358 | 0.187 | -0.604 | 0.621 | -0.099 | 0.205
---|---|---|---|---|---|---|---|---|---|---
$y$ | -0.385 | 0.199 | -0.163 | -0.464 | 0.765 | -0.167 | -0.473 | 0.409 | 0.736 | -0.456

##  V Computational Complexity

The worst-case complexity of solving an SDP problem in each iteration is
$O(m^{2}n^{2})$ [51, 52], where $m$ is the number of the equality constraints,
$n$ is the dimension of the SDP cone. The number of iteration is bounded by
$O(\sqrt{n}\mathrm{ln}(1/\varsigma))$, where $\varsigma$ is the solution
precision, $\sqrt{n}$ is the number of iterations caused by barrier parameter.
For the proposed PF-SDP, we have $m=M+3$ and $n=M+2p+3$. Therefore, the
computational complexity of the PF-SDP in (22) is given by

$\text{PF-SDP Complexity}\simeq
O\big{(}(M+3)^{2}N^{2.5}\mathrm{ln}(1/\varsigma))\big{)},$ (44)

where $N=(M+2p+3)$. To choose an appropriate penalty coefficient, the APF-SDP
requires $\kappa$ PF-SDP solutions, where $\kappa$ depends on $\eta_{0}$, the
initial penalty coefficient, $\eta^{*}$, the optimal penalty coefficient, and
the step length $\gamma$. Therefore, it needs to be satisfied that
$\eta_{0}\gamma^{\kappa}\geq\eta^{*}$. We can further obtain

$\kappa\geq\frac{\mathrm{ln}\eta^{*}-\mathrm{ln}{\eta^{0}}}{\mathrm{ln}{\gamma}}.$
(45)

The complexity of APF-SDP in Algorithm 1 is $\kappa$ times of that of the PF-
SDP and given by

$\text{PF-SDP Complexity}\simeq
O\big{(}\kappa(M+3)^{2}N^{2.5}\mathrm{ln}(1/\varsigma))\big{)}.$ (46)

##  VI Performance Analysis

The RSDP, PF-SDP and APF-SDP algorithms are proposed to predict the initial
position and velocity of the mobile source using the time delay measurements.
To evaluate the performance of these proposed algorithms, the simulations were
firstly conducted in a 2-D scenario. We have also done the simulations for the
3-D case and the observations were similar. The positions of ten sensors are
randomly generated according to the uniform distribution $\mathcal{U}(-1,1)$
and listed in Tab. I. In the same region, a mobile source starts from the
initial position (0.2, -0.4) km with a constant velocity (-1, 1) m/s. The
propagation speed $c$ is randomly drawn from the range $(0.3,0.4)$ km/s. TD
measurements are generated based on the model equation (1), and the noise in
(6) is modeled as zero mean Gaussian distribution with covariance
$\bm{Q}_{n}=\sigma^{2}\bm{I}_{M}$. The performance is evaluated using Mean
Square Error (MSE) defined by

$\left\\{\begin{array}[]{ll}\mathrm{MSE}(\bm{u})=\frac{1}{K}\sum_{k=1}^{K}\|\bm{u}_{k}-\bm{u}^{o}\|^{2}\\\
\mathrm{MSE}(\bm{v})=\frac{1}{K}\sum_{k=1}^{K}\|\bm{v}_{k}-\bm{v}^{o}\|^{2},\end{array}\right.$
(47)

where $\bm{u}^{o}$ and $\bm{v}^{o}$ represent the true initial position and
velocity, $\bm{u}_{k}$ and $\bm{v}_{k}$ are the predicted results in the $k$th
Monte Carlo (MC) run, $K$ is the number of MC runs ($K$ is set to 1000 in our
simulations). The proposed algorithms are also compared with the Cramér-Rao
Lower Bound (CRLB) of the prediction problem.

###  VI-A Penalty Coefficient

The increasing of $\eta$ ensures that the rank of the SDP solution gradually
tends to be one. However, the performance of the rank-one solution may be not
optimal due to the occurrence of excessive penalty. Therefore, it is required
to ensure that the rank-one solution has no the excessive penalty. To evaluate
the improved performance of PF-SDP compared with the RSDP, the Logarithmic MSE
difference (LOGMSED) is defined by

$\mathrm{LOGMSED}=\mathrm{10log(MSE})_{\mathrm{PF}}-\mathrm{10log(MSE})_{\mathrm{R}},$
(48)

where $\mathrm{10log(MSE})_{\mathrm{PF}}$ and
$\mathrm{10log(MSE})_{\mathrm{R}}$ represent the MSE in log-scale of PF-SDP
and RSDP, respectively.

The first eight sensors listed in Tab. I are used to determine the initial
position and velocity of the mobile source. The noise level
$10\mathrm{log}\sigma^{2}$ is set to -40, -20, 0, respectively. Fig. 3(a)
illustrates the LOGMSED($\bm{u}$) performance when the penalty coefficient
$\eta$ is increased from $10^{-6}$ to $10^{3}$ (i.e. log$\eta$ is increased
from -6 to 3). The penalty terms do not basically work when $\eta$ is too
small. As can be seen that LOGMSED($\bm{u}$) is near zero when log$\eta$ is
set to -6. When $\eta$ gradually increases, LOGMSED($\bm{u}$) becomes a
negative value, which illustrates that the MSE of PF-SDP is smaller than that
of the RSDP due to the working of penalty terms. When log$\eta$ is set to (0,
2) for $10\mathrm{log}\sigma^{2}=-40$, (-1, 1) for
$10\mathrm{log}\sigma^{2}=-20$, and (-2, 0) for $10\mathrm{log}\sigma^{2}=0$,
LOGMSED($\bm{u}$) reaches the least value, indicating the optimal performance
of PF-SDP. Therefore, the optimal range of $\eta$ is variable due to the
different noise level. When log$\eta$ is larger than the optimal range,
LOGMSED($\bm{u}$) will be sharply increased, which confirms the occurrence of
excessive penalty. For instance, when log$\eta$ is set to (2,3),
LOGMSED($\bm{u}$) becomes a positive value for $10\mathrm{log}\sigma^{2}=0$,
illustrating that the performance of PF-SDP is worse than that of the RSDP.

When the parameter setup is the same with that of Fig. 3(a), Fig. 3(b) shows
the LOGMSED($\bm{v}$) with the increasing of $\eta$. It is also observed that
the LOGMSED($\bm{v}$) becomes a negative value when log$\eta$ is gradually
increased from -6. Then the LOGMSED($\bm{v}$) reaches a least value that
manifests the optimal performance of PF-SDP. For instance, when
$10\mathrm{log}\sigma^{2}$ is set to -40, the optimal log$\eta$ is set to the
range of (0, 2), where the PF-SDP provides better performance and has almost
13 dB bias compared with the RSDP. When log$\eta$ is larger than the optimal
range, the LOGMSED($\bm{v}$) will be sharply increased due to the occurrence
of excessive penalty.

(a) logarithmic MSE difference of position.

(b) logarithmic MSE difference of velocity.

Figure 3: Performance comparison as the $\eta$ varies. TABLE II: Number of rank-one solutions when $\eta$ varies(1000 MC runs for each value of $\eta$, $10\mathrm{log}\sigma^{2}=0$) log$\eta$ | -6 | -5 | -4 | -3 | -2 | -1 | 0 | 1 | 2
---|---|---|---|---|---|---|---|---|---
Number | 0 | 3 | 4 | 12 | 28 | 306 | 825 | 1000 | 1000
TABLE III: Number of rank-one solutions when $\sigma^{2}$ varies (1000 MC runs for each value of $10\mathrm{log}\sigma^{2}$) 10log$\sigma^{2}$ | -60 | -50 | -40 | -30 | -20 | -10 | 0 | 10 | 20
---|---|---|---|---|---|---|---|---|---
Number ($\mathrm{log}{\eta}=-1$) | 0 | 1 | 9 | 37 | 69 | 169 | 293 | 525 | 739
Number ($\mathrm{log}{\eta}=0$) | 2 | 34 | 112 | 354 | 531 | 619 | 834 | 999 | 1000
Number ($\mathrm{log}{\eta}=1$) | 205 | 542 | 675 | 769 | 872 | 999 | 1000 | 1000 | 1000

(a) Performance of predicted position.

(b) Performance of predicted velocity.

Figure 4: Performance comparison as $\sigma^{2}$ varies.

The rank of PF-SDP solution will gradually tend to be one as $\eta$ increases.
If meeting $\tau<\delta$, it is considered as a rank-one solution. When
10log$\sigma^{2}$ is set to 0 dB, the number of rank-one solutions is shown in
Tab. II, where $\delta$ is set to $10^{-5}$. As can be seen that the number of
rank-one solutions is increased with log$\eta$, since larger $\eta$ provides
tighter SDP cone constraints. When log$\eta$ is larger than 1, all 1000 MC
runs of PF-SDP solutions meet the rank-one condition. When 10log$\sigma^{2}$
is increased from -60 dB to 20 dB, the number of rank-one solutions is also
illustrated in Tab. III, where log$\eta$ is set to -1, 0, and 1, respectively.
The number of rank-one solutions is also dramatically increased with
10log$\sigma^{2}$. When 10log$\sigma^{2}$ is set to -60 dB, the number of
rank-one solutions is 2 for log$\eta=0$, 205 for log$\eta=1$. However, all
1000 MC runs meet the rank-one condition for log$\eta=0$ and log$\eta=1$ when
10log$\sigma^{2}$ is increased to 20 dB.

###  VI-B Performance with Simulations

(a) Performance of predicted position.

(b) Performance of predicted velocity.

Figure 5: Performance comparison as $c$ varies.

(a) Performance of predicted position.

(b) Performance of predicted velocity.

Figure 6: Performance comparison as number of sensors varies.

The parameters in APF-SDP are listed as follows, $\alpha=5$, $\gamma=10$, and
$\delta=10^{-5}$. The first eight sensors listed in Tab. I are used to locate
the mobile source which starts from the initial position (0.2, -0.4) km with a
constant velocity (-1, 1) m/s. When the noise level 10log$\sigma^{2}$ is
increased from -60 dB to 20 dB, Fig. 4(a) and Fig. 4(b) illustrates the
logarithmic MSEs of predicted initial position and velocity, respectively. It
is obviously shown that the logarithmic MSE performance becomes worse as
10log$\sigma^{2}$ increases. When log$\eta$ is set to -6 and -4, the
logarithmic MSEs of PF-SDP are almost same with that of the RSDP, indicating
the no working of penalty terms. However, the performance of PF-SDP becomes
better when log$\eta$ is set to -2. Especially, the logarithmic MSE of PF-SDP
is closer to the CRLB at large noise level, since the rank of PF-SDP solution
is easier to be one. Compared with the RSDP, the AFP-SDP provides better
performance in reaching the CRLB accuracy. Hence, it confirms the advantage of
the APF-SDP in the position and velocity prediction.

(a) The scene in real experiments.

(b) CDF of position error

(c) CDF of velocity error.

Figure 7: Performance comparison in real experiments.

The noise level 10log$\sigma^{2}$ is set to -40 dB and the other parameters
are the same with those in Fig. 4. Fig. 5(a) and Fig. 5(b) also shows the
logarithmic MSEs of predicted initial position and velocity when the
propagation speed $c$ is increased from 0.1 km/s to 1 km/s. It is shown that
the performance of predicted position or velocity becomes worse as the
propagation speed increases. The performance of PF-SDP is always better that
of the RSDP when log$\eta$ is set to -2, indicating the working of penalty
terms. The APF-SDP always provides the almost CRLB performance in the
prediction of position and velocity when $c$ is varied from 0.1 km/s to 1km/s.
It confirms the advantages of AFP-SDP which can adaptively choose the penalty
coefficient and provides optimal performance at the entire range of $c$.

Finally, we also examine the impact of sensors on the performance of these
proposed algorithms when the sensors are selected from Tab. I in the order of
list. The other parameter setup is also the same with that in Fig.4. Fig. 6(a)
and Fig. 6(b) show the logarithmic MSEs of predicted initial position and
velocity as the number of sensors varies. As can be seen that the logarithmic
MSEs of PF-SDP (log$\eta$ = -2) are the same with those of RSDP at when there
are only five or six sensors. However, the PF-SDP (log$\eta$ = -2) provides
better performance compared with the RSDP when the number of sensors is larger
than seven. The proposed APF-SDP always reach the CRLB performance. It is also
shown that adding more sensors can reduce the logarithmic MSE. However, how
much reduction in logarithmic MSE by using additional sensors depends on the
positions of the new sensors included.

###  VI-C Real Experiments

To verify the performance of our proposed algorithms, we also conducted the
real experiments using a mobile car equipped with Ultrasonic module (UM).
Besides the UM, motion sensors are also equipped to the mobile car and used to
measure the parameters including the velocity and moving direction. Nine
sensors are placed at the positions listed in Tab. IV. The UM transmits the
ultrasonic signals to the sensors at the initial position, then the signals
are received by the sensors and reflected to the UM of mobile car. Extensive
tests show that the measurement noise in the range-difference is sufficiently
close to be zero-mean Gaussian distribution with an predicted noise variance
$6.3\times 10^{-3}$ ms2. The initial position of the mobile car is set to (1,
2) m, and the velocity of the mobile car is in the range of (1, 3) m/s.

TABLE IV: Positions of nine sensors in real experiment(m) $x$ | 5 | 5 | -5 | -5 | 5 | -5 | 0 | 0 | 0
---|---|---|---|---|---|---|---|---|---
$y$ | 5 | -5 | 5 | -5 | 0 | 0 | 5 | -5 | 0

(a) RMSE performance of predicted position.

(b) RMSE performance of predicted velocity.

Figure 8: Performance comparison using experimental data as number of sensors
varies.

We collect 200 sampling data that is used to predict the initial position and
and velocity of the mobile car. Fig. 7(b) illustrates the cumulative
distribution function (CDF) of position error with 200 runs with these
different algorithms. As can be seen that 20% of the position error is larger
than 0.05 m for RSDP. However, it is reduced to near 0.01 m for the APF-SDP,
indicating the advantage of our APF-SDP in the performance of position
prediction. Fig. 7(c) shows the CDF of velocity error using 200 sampling data.
20% of the velocity error is larger than 0.08 m/s for RSDP, but it is reduced
to 0.06 m/s for our proposed APF-SDP. Therefore, the APF-SDP performs better
than the RSDP in the prediction of velocity, confirming the advantage of the
APF-SDP by achieving rank-one constraint.

Selecting the sensors from Tab. IV in the order of list, we also verify the
performance of these proposed algorithm when the number of sensors is
increased from 5 to 9. To clearly illustrate the effect of these different
algorithms, we use the root mean square error (RMSE) to evaluate the
performance with 200 sampling data. The position RMSE performance of these
algorithms is plotted in Fig. 8(a), where “CRLB expected” represents the best
achievable performance expected using the predicted noise variance. The RSDP
performs not well even if the number of sensors increases. However, the
proposed APF-SDP almost reaches the “CRLB expected” performance, which is
consistent with the simulation results. Fig. 8(b) illustrates the velocity
RMSE as the number of sensors increases. The performance of PF-SDP is better
than that of the RSDP when log$\eta$ is set to -2. However, the PF-SDP can not
provide the “CRLB expected” performance since it can not adaptively choose the
penalty coefficient.

###  VI-D Industrial Applications

There are many position and velocity prediction problems in industrial
Internet of things, such as the position obtaining of mobile AV and the
tracking of UAV, which are illustrated in Fig. 9(a) and Fig. 9(b),
respectively. Most of the current research on these problems is focused on the
position obtaining of mobile source using some ranging information or the
velocity prediction using motion sensors or Dropper shift measurements. The
motion sensor is an extra hardware device which will increase the system cost.
Since the velocity obtaining of Dropper shift measurement is subject to the
error of (2, 4) m/s, it is considered to be unacceptable in some application
systems. Our proposed method does not require any motion sensors or Dropper
shift measurements and realize the prediction of position along with the
velocity of mobile source.

Our experimental results show that the mean position RMSE of APF-SDP is
smaller than 0.01 m when the UM equipped to the mobile source is used to
collect the TD information. Apparently, the signal can be acoustic or
electromagnetic, and the signal can travels in the underwater or underground
scenarios. Due to the different propagation speed of the signal, our proposed
algorithms provides distinct accuracy performance in the prediction of
position and velocity, which is also demonstrated in Fig. 5. Therefore, a
feasible method to improve the accuracy performance is to reduce the noise
level at large propagation speed. Moreover, the precise timing is very
important to ensure the performance especially when the propagation speed of
electromagnetic signal reaches $3\times 10^{8}$ m/s.

(a) 2-D scenario of mobile AV.

(b) 3-D scenario of mobile UAV.

Figure 9: Different scenarios of our proposed system in industrial
applications.

##  VII Conclusions and Future Work

We introduce an intelligent practice prediction system for mobile source using
the TD measurements. To predict the position and velocity of mobile source,
the RSDP algorithm is firstly proposed to obtain a convex SDP problem by
dropping the rank-one constraint. The performance of RSDP is very poor due to
the drop of rank-one constraint. Then the PF-SDP algorithm is put forward to
obtain better performance by introducing the penalty terms. In the PF-SDP, the
penalty coefficient is crucial to ensure the well performance of the PF-SDP
and its optimal value is variable. Therefore, the AFP-SDP algorithm is
proposed by adaptively choosing the penalty coefficient. Compared with the
RSDP and PF-SDP, the AFP-SDP shows its performance advantages in the
prediction of position and velocity. We have also done the simulations and the
real experiments for the 2-D case. For the future work, our proposed methods
are not applied in the 3-D scenario. Therefore, we will focus the future work
on the position prediction of UAV and extend the system to the 3-D case.

## References

  * [1] Y. Wang, A. V. Vasilakos, J. Ma, and N. Xiong, “On Studying the Impact of Uncertainty on Behavior Diffusion in Social Networks,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , vol. 45, no. 2, pp. 185–197, 2015.
  * [2] Y. Zhou, D. Zhang, and N. Xiong, “Post-cloud Computing Paradigms: A Survey and Comparison,” _Tsinghua Science and Technology_ , vol. 12, no. 12, pp. 714–732, 2017.
  * [3] H. Zheng, W. Guo, and N. Xiong, “A Kernel-Based Compressive Sensing Approach for Mobile Data Gathering in Wireless Sensor Network Systems,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , vol. 48, no. 12, pp. 2315–2327, 2018.
  * [4] X. Wu, S. Wang, H. Feng, J. Hu, and G. Wang, “Motion Parameter Capturing of Multiple Mobile Targets in Robotic Sensor Networks,” _IEEE Access_ , vol. 6, pp. 24 375–24 390, 2018.
  * [5] M. Wu, N. Xiong, and L. Tan, “Adaptive Range-Based Target Localization Using Diffusion Gauss-Newton Method in Industrial Environments,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 11, pp. 5919–5930, 2019\.
  * [6] H. Huang, A. V. Savkin, M. Ding, and C. Huang, “Mobile Robots in Wireless Sensor Networks in Wireless Snesor Networks: Survey on Tasks,” _Computer Networks_ , vol. 148, no. 2, pp. 1–19, 2019.
  * [7] L. Shu, Y. Zhang, Z. Yu, L. T. Yang, M. Hauswirth, and N. Xiong, “Context-aware Cross-layer Optimized Video Streaming in Wireless Multimedia Sensor Networks,” _The Journal of Supercomputing_ , vol. 54, no. 1, pp. 94–121, 2010.
  * [8] W. Fang, Y. Li, H. Zhang, N. Xiong, J. Lai, and A. V. Vasilakos, “On the Throughput-Energy Tradeoff for Data Transmission Between Cloud and Mobile Devices,” _Information Sciences_ , vol. 283, pp. 79–93, 2014.
  * [9] Y. Liu, M. Ma, X. Liu, N. N. Xiong, A. Liu, and Y. Zhu, “Design and Analysis of Probing Route to Defense Sink-Hole Attacks for Internet of Things Security,” _IEEE Transactions on Network Science and Engineering_ , vol. 7, no. 1, pp. 356–372, 2020.
  * [10] J.-M. Pak, C. K. Ahn, P. Shi, Y. S. Shmaliy, and M. T. Lim, “Distributed Hybrid Particle/FIR Filtering for Mitigating NLOS Effects in TOA-Based Localization Using Wireless Sensor Networks,” _IEEE Transactions on Industrial Electronics_ , vol. 64, no. 6, pp. 5182–5191, 2017.
  * [11] Y. Kang, Q. Wang, J. Wang, and R. Chen, “A High-Accuracy TOA-Based Localization Method Without Time Synchronization in a Three-Dimensional Space,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 1, pp. 173–182, 2019.
  * [12] B. Huang, L. Xie, and Z. Yang, “TDOA-Based Source Localization With Distance-Dependent Noises,” _IEEE Transactions on Wireless Communications_ , vol. 14, no. 1, pp. 468–480, 2015.
  * [13] Y. Sun, K. C. Ho, and Q. Wan, “Solution and Analysis of TDOA Localization of a Near or Distant Source in Closed Form,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 2, pp. 320–335, 2019.
  * [14] Y. Zhang, S. Xing, Y. Zhu, F. Yan, and L. Shen, “RSS-Based Localization in WSNs Using Gaussian Mixture Model via Semidefinite Relaxation,” _IEEE Communication Letters_ , vol. 21, no. 6, pp. 1329–1332, 2017.
  * [15] X. Wu, X. Zhu, S. Wang, and L. Mo, “Cooperative Motion Parameter Estimation Using RSS Measurements in Robotic Sensor Networks,” _Journal of Network and Computer Applications_ , vol. 136, pp. 57–70, 2019.
  * [16] H.-J. Shao, X.-P. Zhang, and Z. Wang, “Efficient Closed-Form Algorithms for AOA Based Self-Localization of Sensor Nodes Using Auxiliary Variables,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 10, pp. 2580–2594, 2014.
  * [17] A. Amar, G. Leus, and B. Friedlander, “Emitter Localization Given Time Delay and Frequency Shift Measurements,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 48, no. 2, pp. 1826–1837, 2012.
  * [18] C.-H. Park and J.-H. Chang, “Closed-Form Localization for Distributed MIMO Radar Systems Using Time Delay Measurements,” _IEEE Transactions on Wireless Communications_ , vol. 15, no. 2, pp. 1480–1490, 2016.
  * [19] T. Jia, K. C. Ho, H. Wang, and X. Shen, “Effect of Sensor Motion on Time Delay and Doppler Shift Localization: Analysis and Solution,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 22, pp. 5881–5895, 2019.
  * [20] ——, “Accurate Localization of AUV in Motion by Explicit Solution Using Time Delays,” in _2020 IEEE International Conference on Acoustics, Speech and Signal Processing_ , 05 2020, pp. 4871–4875.
  * [21] N. H. Nguyen and K. Dogancay, “Closed-Form Algebraic Solutions for 3-D Doppler-Only Source Localization,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 10, pp. 6822–6836, 2018.
  * [22] S. Yu, X. Su, and J. Zeng, “Cooperative Localization with Velocity-Assisted for Mobile Wireless Sensor Networks,” in _18th International Symposium on Communications and Information Technologies_ , 09 2018, pp. 131–136.
  * [23] L. Carlino, F. Bandiera, A. Coluccia, and G. Ricci, “Improving Localization by Testing Mobility,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 13, pp. 3412–3423, 2019.
  * [24] H. Chen, F. Li, and Y. Wang, “SoundMark: Accurate Indoor Localization via Peer-Assisted Dead Reckoning,” _IEEE Internet of Things Journal_ , vol. 5, no. 6, pp. 4803–4815, 2018.
  * [25] Y. Yu, “Consensus-Based Distributed Mixture Kalman Filter for Maneuvering Target Tracking in Wireless Sensor Networks,” _IEEE Transactions on Vehicular Technology_ , vol. 65, no. 10, pp. 8669–8681, 2016.
  * [26] L. Li, M. Yang, C. Wang, and B. Wang, “Hybrid Filtering Framework Based Robust Localization for Industrial Vehicles,” _IEEE Transactions on Industrial Informatics_ , vol. 14, no. 3, pp. 941–950, 2018.
  * [27] L. Yang, L. Yang, and K. C. Ho, “Moving Target Localization in Multistatic Sonar by Differential Delays and Doppler Shifts,” _IEEE Signal Processing Letters_ , vol. 23, no. 9, pp. 1160–1164, 2016.
  * [28] T. Tirer and A. J. Weiss, “High Resolution Localization of Narrowband Radio Emitters Based on Doppler Frequency Shifts,” _Signal Processing_ , vol. 141, pp. 288–298, 2017.
  * [29] A. Lemon, A. M.-C. So, and Y. Ye, “Low-Rank Semidefinite Programming: Theory and Applications,” _Foundations and Trends in Optimization_ , vol. 106, pp. 51–60, 2016.
  * [30] G. Wang and K. C. Ho, “Convex Relaxation Methods for Unified Near-Field and Far-Field TDOA-Based Localization,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 4, pp. 2346–2360, 2019.
  * [31] H. Qi, X. Wu, and L. Jia, “Semidefinite Programming for Unified TDOA-based Localization Under Unknown Propagation Speed,” _IEEE Communication Letters_ , 2020.
  * [32] A. Simonetto and G. Leus, “Distributed Maximum Likelihood Sensor Network Localization,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 6, pp. 1424–1437, 2014.
  * [33] E. Dranka and R. F. Coelho, “Robust Maximum Likelihood Acoustic Energy Based Source Localization in Correlated Noisy Sensing Environments,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 9, no. 2, pp. 259–267, 2015.
  * [34] J. Liang, Y. Chen, H. C. So, and Y. Jing, “Circular/hyperbolic/elliptic Localization via Euclidean Norm Elimination,” _Signal Processing_ , vol. 148, pp. 102–113, 2018.
  * [35] C. Zhang and Y. Wang, “Sensor Network Event Localization via Nonconvex Nonsmooth ADMM and Augmented Lagrangian Methods,” _IEEE Transactions on Control of Network Systems_ , vol. 6, no. 4, pp. 1473–1485, 2019.
  * [36] T.-K. Le and N. Ono, “Closed-Form and Near Closed-Form Solutions for TDOA-Based Joint Source and Sensor Localization,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 5, pp. 1207–1221, 2017.
  * [37] S. Tomic, M. Beko, and M. Tuba, “A Linear Estimator for Network Localization Using Integrated RSS and AOA Measurements,” _IEEE Signal Processing Letter_ , vol. 26, no. 3, pp. 405–409, 2019.
  * [38] G. Wang, A. M.-C. So, and Y. Li, “Robust Convex Approximation Methods for TDOA-Based Localization Under NLOS Conditions,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 13, pp. 3281–3296, 2016.
  * [39] J. Zheng and X. Wu, “Convex Optimization Algorithms for Cooperative RSS-based Sensor Localization,” _Pervasive and Mobile Computing_ , vol. 37, pp. 78–93, 2017.
  * [40] Z. Su, G. Shao, and H. Liu, “Semidefinite Programming for NLOS Error Mitigation in TDOA Localization,” _IEEE Communications Letters_ , vol. 22, no. 7, pp. 1430–1433, 2018.
  * [41] Y. Zhang, S. Xing, Y. Zhu, F. Yan, and L. Shen, “RSS-Based Localization in WSNs Using Gaussian Mixture Model via Semidefinite Relaxation,” _IEEE Communication Letters_ , vol. 21, no. 6, pp. 1329–1332, 2017.
  * [42] Z. Wang, H. Zhang, T. Lu, and T. A. Gulliver, “Cooperative RSS-Based Localization in Wireless Sensor Networks Using Relative Error Estimation and Semidefinite Programming,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 1, pp. 483–497, 2019.
  * [43] W. Wang, G. Wang, F. Zhang, and Y. Li, “Second-Order Cone Relaxation for TDOA-Based Localization Under Mixed LOS/NLOS Conditions,” _IEEE Signal Processing Letters_ , vol. 23, no. 12, pp. 1872–1876, 2016.
  * [44] C. Helmberg, F. Rendl, R. J. Vanderbei, and H. Wolkowicz, “An Interior-Point Method for Semidefinite Programming,” _SIAM Journal on Optimization_ , vol. 6, no. 2, pp. 342–361, 1996.
  * [45] S. Burer and R. D. C. Monteiro, “Local Minima and Convergence in Low-Rank Semidefinite Programming,” _Mathematical Programming_ , vol. 103, pp. 427–444, 2005.
  * [46] S. Burer and C. Choi, “Computational Enhancements in Low-rank Semidefinite Programming,” _Optimization Methods and Software_ , vol. 21, no. 3, pp. 493–512, 2006.
  * [47] G. Yuan, Z. Zhang, B. Ghanem, and Z. Hao, “Low-Rank Semidefinite Programming: Theory and Applications,” _Neurocomputing_ , pp. 1–156, 2013.
  * [48] M. Wu, L. Zhong, and N. Xiong, “The Distributed Gauss-Newton Methods for Solving the Inverse of Approximated Hessian with Application to Target Localization,” in _Proceedings of the 2020 3rd International Conference on Computer Science and Software_ , 2020, pp. 22–26.
  * [49] M. Journee, F. R. Bach, P.-A. Absil, and R. Sepulchre, “Low-Rank Optimization on the Cone of Positive Semidefinite Matrices,” _SIAM Journal on Optimization_ , vol. 20, no. 5, pp. 2327–2351, 2010.
  * [50] R. Y. Zhang and J. Lavaei, “Modified Interior-point Method for Large-and-sparse Low-rank Semidefinite Programs,” in _2017 IEEE 56th Annual Conference on Decision and Control_ , 2017, pp. 5640–5647.
  * [51] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge, U.K.: Cambridge University Press, 2004.
  * [52] S. Arora and B. Barak, _Computational Complexity: A Modern Approach_. Cambridge, U.K.: Cambridge University Press, 2009.


0\end{pmatrix},\qquad\Psi_{t}(\xi_{1})=\int_{0}^{t}\mathrm{d}\tau\,\psi\left(g\varkappa
Y(\tau)\right),$
$\Psi_{\Delta}(\xi_{2})\simeq\Delta\int_{0}^{1}\mathrm{d}x\,\psi\left(\frac{\varkappa}{\Delta}\,S_{x\Delta}P_{0}\begin{pmatrix}0\\\
1\end{pmatrix}\right)\simeq\Delta\int_{0}^{1}\mathrm{d}x\,\psi\left(\varkappa\left(\frac{1}{\Delta}\,+xZ\right)P_{0}\begin{pmatrix}0\\\
1\end{pmatrix}\right)\\\
=\Delta\int_{0}^{1}\mathrm{d}x\,\psi\left(\frac{\varkappa}{\Delta}\,P_{0}\begin{pmatrix}0\\\
1\end{pmatrix}+\varkappa gxP_{1}\begin{pmatrix}1\\\
0\end{pmatrix}\right)=-\frac{A_{44}}{2\Delta}\,\varkappa^{2}-\frac{A_{11}}{6}\,\Delta
g\varkappa^{2}.$

We also change the initial time from 0 to $t_{0}$, which means
$\hat{\pi}_{0}\to\hat{\pi}_{t_{0}}$. Then (99) becomes

$\mathbb{E}\left[\exp\left\\{\mathrm{i}\varkappa\,\frac{X_{2}(t_{0}+t+\Delta)-X_{2}(t_{0}+t)}{\Delta}\right\\}\right]\\\
{}\simeq\int_{\mathbb{R}^{2}}\mathrm{d}x\operatorname{Tr}\left\\{\hat{\pi}_{t_{0}}(x)W_{1}(g\varkappa
P_{1}Y(t))\right\\}\exp\left\\{\mathrm{i}x_{1}gY_{3}(t)\varkappa-\frac{A_{44}}{2\Delta}\,\varkappa^{2}+\int_{0}^{t}\mathrm{d}\tau\,\psi\big{(}g\varkappa
Y(\tau)\big{)}\right\\}.$ (114)

Note that there is no dependence on $x_{2}$ inside the exponential which means
that these probabilities do not depend on the distribution of $X_{2}(t_{0})$.
For large $t$, $Y(t)$ vanishes and any memory of the state $\hat{\pi}_{t_{0}}$
is lost. For $t=0$ we get $Y_{i}(0)=\delta_{i1}$ and (114) becomes

$\mathbb{E}\left[\exp\left\\{\mathrm{i}\varkappa\,\frac{X_{2}(t_{0}+\Delta)-X_{2}(t_{0})}{\Delta}\right\\}\right]\\\
{}\simeq\operatorname{Tr}\left\\{\rho_{t_{0}}\mathrm{e}^{\mathrm{i}Qg\varkappa}\right\\}\exp\left\\{-\frac{A_{44}}{2\Delta}\,\varkappa^{2}\right\\},\qquad\rho_{t_{0}}=\int_{\mathbb{R}^{2}}\mathrm{d}x\,\hat{\pi}_{t_{0}}(x).$
(115)

This means that the distribution of the increment per unit of time is the
quantum distribution of the position in the state $\rho_{t_{0}}$ smoothed by a
large Gaussian. The state $\rho_{t_{0}}$ depends on the whole previous
history, and it could be the conditional state with respect to previous
observations (an argument not developed here).

##  7 Conclusions

In this article we have afforded the problem of giving a unified treatment of
interacting classical and quantum systems. The formalism we adopted allows to
unify quantum master equations, Liouville equation, Kolmogorov-Fokker-Planck
equation. The treatment has been restricted to the case of a Markovian and
quasi-free dynamics. The notion of quasi-free dynamics, Sec. 2.3, generalizes
the well known Gaussian case by including contributions of “jump type”. The
problem of finding the most general Markovian quasi-free dynamics has been
solved, see Theor. 1; in the construction, connections with the classical
Lévy-Khintchine formula have been used, Sec. 2.4.3.

The restriction to the quasi-free case allows to give explicitly the dynamical
map, without the need of solving the evolution equations, involving unbounded
generators. However, also the formal structure of the generator is given (Sec.
3), because this gives hints to construct hybrid dynamics out of the quasi-
free case. As a byproduct, we have obtained the most general quasi-free
quantum dynamical semigroup, together with the structure of its generator, and
we have presented examples which underline the role of the “jump”
contributions (Sec. 4).

The results in the pure classical case (Sec. 5) are all essentially known.
However, this case allows to introduce the idea that not only probabilities at
one time are involved in the theory, but also transition probabilities, multi-
time probabilities, probabilities for stochastic processes in continuous time.
In the general case (Sec. 6), these ideas are translated in the introduction
of _transition instruments_ (Sec. 6.1), which allow for the construction of
multi-time probabilities for the classical component. By analyzing the
generator of the hybrid dynamics (Sec. 3.2), it is possible to identify the
interactions responsible of the action of the classical system on the quantum
one and the interactions responsible of the flux of information from the
quantum system to the classical one. The classical component can be observed,
in principle without disturbance, and information about the quantum system can
be collected. This last feature connects dynamical semigroups for hybrid
systems to the theory of quantum measurements in continuous time. Let us
stress that the whole construction respects the general idea that no
information can be extracted from a quantum system without dissipation;
indeed, by asking the vanishing of the dissipative terms in the quantum part
of the generator, we get that also the terms responsible of the flux of
information towards the classical component have to vanish. Finally, in Sec.
6.2, an example is elaborated, which shows how a classical system can inject
noise of Gaussian and jump type into the quantum one, and how it can extract
information from the quantum component.

Many points are left open for future developments. In the theory of open
systems, the notion of _stochastic dilation_ or _stochastic unraveling_ has
been introduced, to mean the case of quantum master equations written in form
of (classical) stochastic differential equations [42, 13, 31]. This can be
done also for the hybrid case, as shown in [41], but not all possibilities
have been explored. Another old idea for hybrid systems is to treat the
classical component as quantum, but with only commuting operators appearing as
observables. The connections of Lévy processes with Bose quantum fields and
_quantum stochastic calculus_ have been worked out in detail [32, 43]. By
using these means, in the case of the monitoring in time of a quantum system,
the observed signal has been dilated to a quantum system and represented by
quantum observables, commuting among themselves also at differen times in the
Heisenberg picture [44, 45]. The same construction for the dynamics of general
hybrid systems is to be done, as it is open the problem of not quasi-free
hybrid systems, possibly with unbounded generators [46, 25].

Inside the theory of quantum measurements in continuous time, also the notions
of _quantum filtering_ , _quantum trajectories_ , _conditional states_ have
been developed [42, 13, 45, 50, 49, 47, 48]; the introduction of these notions
for hybrid systems is clearly possible and we expect this to be fruitful. The
classical component of the hybrid can observed, or partially observed, and the
whole filtering theory for stochastic processes enters into play. As we wrote,
in the hybrid case there is transfer of information also from the classical to
the quantum component and it could be interesting to connect this to the idea
of _feedback_ on quantum systems [13]. Obviously, the final open problem is to
find physical application of the new ideas on the dynamics of hybrid systems
in quantum optics and quantum information, beyond the quasi-free case.

## References

  * [1] L. Diósi, N. Gisin, W.T. Strunz, Quantum approach to coupling classical and quantum dynamics, Phys. Rev. A 61 (2000) 022108.
  * [2] L. Diósi, Hybrid quantum-classical master equations, Phys. Scr. T163 (2014) 014004.
  * [3] A. Sergi, D. Lamberto, A. Migliore, A. Messina, Quantum–classical hybrid systems and Ehrenfest’s Theorem, Entropy 25 (2021) 602.
  * [4] L. Dammeier, R.F. Werner, Quantum-classical hybrid systems and their quasifree transformations, arxiv:2208.05020 (2022).
  * [5] L. Diósi, Hybrid completely positive Markovian quantum-classical dynamics, Phys. Rev. A 107 (2023) 062206.
  * [6] G. Manfredi, A. Rittaud, C. Tronci, Hybrid quantum-classical dynamics of pure-dephasing systems, J. Phys. A: Math. Theor. (2023).
  * [7] A.S. Holevo, Statistical Structure of Quantum Theory, Lecture Notes in Physics m 67 (Springer, Berlin, 2001).
  * [8] A. Barchielli, G. Lupieri, Instruments and mutual entropies in quantum information, Banach Center Publ. 73 (2006) 65–80.
  * [9] A. Barchielli, G. Lupieri, A quantum analogue of Hunt’s representation theorem for the generator of convolution semigroups on Lie groups, Probab. Theory Rel. Fields 88 (1991) 167–194.
  * [10] A. Barchielli, A. S. Holevo, G. Lupieri, An analogue of Hunt’s representation theorem in quantum probability, J. Theor. Probab. 6 (1993) 231–265.
  * [11] A. Barchielli, A. M. Paganoni, A note on a formula of Lévy–Khinchin type in quantum probability, Nagoya Math. J. 141 (1996) 29–43.
  * [12] A. Barchielli, A. M. Paganoni, F. Zucca, On stochastic differential equations and semigroups of probability operators in quantum probability, Stoch. Proc. Appl. 73 (1998) 69–86.
  * [13] H.M. Wiseman and G.J. Milburn, Quantum Measurement and Control (Cambridge University Press, Cambridge, 2010).
  * [14] K. Sato, Lévy processes and infinitely divisible distributions (Cambridge University Press, Cambridge, 1999)
  * [15] D. Applebaum, Lévy Processes and Stochastic Calculus, Second Edition (Cambridge University Press, 2009).
  * [16] A.S. Holevo, Conditionally positive definite functions in quantum probability, Proc. Int. Cong. Math. 1011-1020 (1986) Berkeley.
  * [17] A. Barchielli, Convolution semigroups in quantum probability, Semesterbericht Funktionalanalysis Tübingen, Band 12, Sommersemester 1987, pp. 29–42.
  * [18] A. S. Holevo, Limit theorems for repeated measurements and continuous measurement processes, in L. Accardi, W. von Waldenfels (eds.), Quantum Probability and Applications IV, Lecture Notes in Mathematics 1396 (Springer, Berlin, 1989) pp. 229–257.
  * [19] A. Barchielli, G. Lupieri, Convolution semigroups in quantum probability and quantum stochastic calculus. In L. Accardi, W. von Waldenfels (eds.), Quantum Probability and Applications IV. Lecture Notes in Mathematics, Vol. 1396 (Springer, Berlin, 1989) pp. 107–127.
  * [20] A. Barchielli, Probability operators and convolution semigroups of instruments in quantum probability, Probab. Theory Rel. Fields 82 (1989) 1–8.
  * [21] A. Barchielli, Some Markov semigroups in quantum probability. In L. Accardi, W. von Waldenfels (eds.), Quantum Probability and Applications V, Lecture Notes in Mathematics, vol. 1442 (Springer, Berlin, 1990) pp. 86–98.
  * [22] A. Barchielli, G. Lupieri, Convolution semigroups of instruments in quantum probability. In Probability Theory and Mathematical Statistics Vol. I, ed. by B. Grigelionis et al. (Mokslas, Vilnius, and VSP, Utrecht, 1990) pp. 78–90.
  * [23] A. Barchielli, On the quantum theory of measurements continuous in time, Rep. Math. Phys. 33 (1993) 21–34.
  * [24] A.M. Chebotarev, F. Fagnola, Sufficient conditions for conservativity of minimal quantum dynamical semigroups, J. Funct. Anal. 153 (1998) 382–404.
  * [25] I. Siemon, A.S. Holevo, and R.F. Werner, Unbounded generators of dynamical semigroups, Open Systems & Information Dynamics 24 1740015 (2017).
  * [26] P. Vanheuverzwijn, Generators for quasi-free completely positive semi-groups, Ann. I. H. P. A 29 (1978) 123–138; Errata: Generators for quasi-free completely positive semi-groups Ann. I. H. P. A 30 (1979) 83.
  * [27] B. Demoen, P. Vanheuverzwijn, A. Verbeure, Completely positive quasi-free maps of the CCR-algebra, Rep. Math. Phys. 15 (1979) 27–39.
  * [28] M. Hellmich, Quasi-free semigroups on the CCR algebra, Rep. Math. Phys. 66 (2010) 277–298.
  * [29] T. Heinosaari, A.S. Holevo, M.M. Wolf, The semigroup structure of Gaussian channels, J. Quantum Inf. Comp., 10:0619-0635 (2010).
  * [30] A.S. Holevo, On translation covariant quantum Markov equations, Izv. Math. 59 (1995) 427–443.
  * [31] A.S. Holevo, Covariant quantum Markovian evolutions, J. Math. Phys. 37 (1996) 1812–1832.
  * [32] K.R. Parthasarathy, An Introduction to Quantum Stochastic Calculus (Birkhäuser, Basel, 1992).
  * [33] A.S. Holevo, A noncommutative generalization of conditionally positive definite functions, in L. Accardi, W. von Waldenfels (eds.), Quantum Probability and Applications III. Lecture Notes in Mathematics, Vol. 1303 (Springer, Berlin, 1988), pp. 128–148.
  * [34] A.V. Skorohod, Random Processes with Independent Increments (Kluwer, Dordrecht, 1991).
  * [35] K. Sato, M. Yamazato, Operator-selfdecomposable distributions as limit distributions of processes of Ornstein-Uhlenbeck type, Stoc. Proc. Appl. 17 (1984) 73–100.
  * [36] B. Vacchini, Description of decoherence by means of translation-covariant master equations and Lévy processes, in J.C. García, R. Quezada, S.B. Sontz, Quantum Probability and Related Topics, QP-PQ: Quantum Probability and White Noise Analysis Vol. 23 (World Scientific, Singapore, 2008) pp. 254–265.
  * [37] B. Vacchini, K. Hornberger, Quantum linear Boltzmann equation, Phys. Rep. 478 (2009) 71–120.
  * [38] B. Vacchini, Quantum optical versus quantum Brownian motion master equation in terms of covariance and equilibrium properties, J. Math. Phys. 43 (2002) 5446–5458.
  * [39] A. Barchielli, Quantum stochastic equations for an opto-mechanical oscillator with radiation pressure interaction and non-Markovian effects, Rep. Math. Phys. 77 (2016) 315–333.
  * [40] A. Barchielli, M. Gregoratti, Quantum optomechanical system in a Mach-Zehnder interferometer, Phys. Rev. A 104 (2021) 013713.
  * [41] A. Barchielli, A. S. Holevo, Constructing quantum measurement processes via classical stochastic calculus, Stoch. Proc. Appl. 58 (1995) 293–317. F
  * [42] A. Barchielli, M. Gregoratti, Quantum Trajectories and Measurements in Continuous Time: The Diffusive Case, Lect. Notes Phys. 782 (Springer, Berlin & Heidelberg, 2009).
  * [43] K.R. Parthasarathy, From CCR to Lévy processes: An excursion in quantum probability, Expo. Math. 36 (2018) 302–324.
  * [44] A. Barchielli, Measurement theory and stochastic differential equations in quantum mechanics, Phys. Rev. A 34 (1986) 1642–1649.
  * [45] A. Barchielli, Continual Measurements in Quantum Mechanics and Quantum Stochastic Calculus. In S. Attal, A. Joye, C.-A. Pillet (eds.), Open Quantum Systems III, Lect. Notes Math. 1882 (Springer, Berlin, 2006), pp. 207–291.
  * [46] R. Castro Santis, A. Barchielli, Quantum stochastic differential equations and continuous measurements: unbounded coefficients, Rep. Math. Phys. 67 (2011) 229–254.
  * [47] V.P. Belavkin, Nondemolition measurements, nonlinear filtering and dynamic programming of quantum stochastic processes. In A. Blaquière (ed.), Modelling and Control of Systems, Lecture Notes in Control and Information Sciences, vol. 121 (Springer, Berlin, 1988) pp. 245–265.
  * [48] V.P. Belavkin, A new wave equation for a continuous nondemolition measurement, Phys. Lett. A 140 (1989) 355-358.
  * [49] A. Barchielli, V. P. Belavkin, Measurements continuous in time and a posteriori states in quantum mechanics, J. Phys. A: Math. Gen. 24 (1991) 1495–1514.
  * [50] P. Zoller and C.W. Gardiner, Quantum noise in quantum optics: the stochastic Schrödinger equation. In S. Reynaud, E. Giacobino & J. Zinn-Justin eds., Fluctuations quantiques, (Les Houches 1995) (North-Holland, Amsterdam, 1997) pp. 79–136.