Sensor Control in Anti-Submarine Warfare—A Digital Twin and Random Finite Sets Based Approach

Abstract

Since the submarine has become the major threat to maritime security, there is an urgent need to find a more efficient method of anti-submarine warfare (ASW). The digital twin theory is one of the most outstanding information technologies, and has been quite popular in recent years. The most influential change produced by digital twin is the ability to enable real-time dynamic interactions between the simulation world and the real world. Digital twin can be regarded as a paradigm by means of which selected online measurements are dynamically assimilated into the simulation world, with the running simulation model guiding the real world adaptively in reverse. By combining digital twin theory and random finite sets (RFSs) closely, a new framework of sensor control in ASW is proposed. Two key algorithms are proposed for supporting the digital twin-based framework. First, the RFS-based data-assimilation algorithm is proposed for online assimilating the sequence of real-time measurements with detection uncertainty, data association uncertainty, noise, and clutters. Second, the computation of the reward function by using the results of the proposed data-assimilation algorithm is introduced to find the optimal control action. The results of three groups of experiments successfully verify the feasibility and effectiveness of the proposed approach.

1. Introduction

Submarines are the main combat forces of modern maritime warfare, and the major threats to maritime security. Anti-submarine warfare (ASW) is a type of warfare that depends on surface warships, aircraft, or submarines to fight against enemy submarines. The key of ASW is to quickly identify and localize as many enemy submarines as possible. Sensor control is the key technology for the victory in ASW, so we focus on the innovation of the online sensor control method. Many works have been done to apply simulation-based approaches in naval warfare research, but there are quite a few effective methods for combining simulation technologies with the real ASW in real time. In this paper, we study how to control the sensor of anti-submarine ships in ASW by employing simulation theory, random finite set (RFS) theory, and digital twin theory.

The sensor control problem is also known as sensor management problem. The sensor equipped on the anti-submarine ship can perform many different actions including moving to designated areas, searching in certain directions, etc. In real ASW, the anti-submarine ship usually takes tactical actions to estimate the accurate distance and guarantee the observability [1,2,3]. Different actions have different effectiveness; some actions will be effective for sensing the submarine, while some others not. Here the goal of sensor control is to ensure maximum efficiency of sensors and provide more accurate measurements to the simulation system. Finding the optimal control action is an urgent need in practical application. In this paper, the objective of sensor control in ASW is choosing control actions online so that the utility of sensors is maximized. Here, sensor control means sequential decision making, where each control decision improves the utility of sensors and helps to win the victory in ASW.

In [4], a practical implementation of sensor control based on Kullback–Leibler divergence is proposed for the first time. In [5], the Rényi divergence is used for deriving the reward function. In [6], Kullback–Leibler discrimination is used as the reward function for sensor control, but the authors did not provide the implementation method. In [7], the posterior expected number of targets (PENT)-based reward function is given. PENT has been successfully applied in some applications [8,9]. In [10], the Cauchy-Schwarz divergence is used, and the analytical solution is also derived. These works all depend on the expected divergence or information gain between the predicted and posterior densities for choosing the optional control action. Their major limitations are the problems of significant computational cost, and too little use of simulation methods. The simulation methods are the key for taking more complex dynamic application scenarios into consideration. It is worth mentioning that some divergences usually fail to be analytically computed.

RFS-based Bayesian inference has an advantage in dealing with situations where the number of targets and received detections are random, and the targets’ states and detections are also random. The RFS-based Bayesian inference can successfully overcome the limitations of the conventional vector-based Bayesian inference, and is gradually applied to the sensor control problem. In [11], the Rényi divergence and RFS-based reward function is derived, and the probability hypothesis density-based filter is adopted for implementing the sensor control method. In [12], the RFS and Rényi divergence-based reward function is derived, and the CB-MemBer filter-based implementation of the sensor control method is proposed. In [13], the RFS-based analytical solution for Cauchy-Schwarz divergence is derived. In [14], the Cauchy-Schwarz divergence and generalized labeled multi-Bernoulli filter-based sensor control method is derived. In [15], a labeled multi-Bernoulli filter and Cauchy-Schwarz divergence-based constrained sensor control method is proposed. In [16], the authors control the sensors by using the expected risk reduction between the multi-target predicted and updated densities.

The framework of partially observed Markov decision process (POMDP) is used to study the problem of making the next decision by using the past decisions and past observations [17]. The current information state, a set of admissible sensor actions and the reward function associated with each action are the main elements of POMDP. The above-mentioned works all adopt the framework of POMDP. POMDP is widely used in robot control, moving target searching, object identification, and so on. POMDP depends on the assumption that the states of the studied system are determined by the Markov process. POMDP is lack of use of the prediction and evaluation ability of simulation theory, so it fails to take many complex simulation scenarios into consideration. Since it focuses on the basic mathematical theory, its ability to apply new information technologies is also weak. POMDP depends on the connection of the virtual space and the physical space, but it fails to tell us how to connect these different spaces. In this paper, we propose to use a more promising framework to solve the problem of sensor control in ASW.

Recently, digital twin is seen as the best way to enable the interoperability and integration between the real and simulated worlds [18]. At present, there is no formal and consistent definition of digital twin, and there are many different versions of definitions of digital twin. Different industries and application fields have different definition perspectives and methods. However, the core ideas of all definitions are similar. Grieves defines the digital twin as “a set of virtual information constructs that fully describes a potential or actual physical manufactured product from the micro atomic level to the macro geometrical level” [19,20,21]. Shafto et al. define the digital twin as “an integrated multiphysics, multiscale simulation of a vehicle or system that uses the best available physical models, sensor updates, fleet history, etc., to mirror the life of its corresponding flying twin” [22,23]. These definitions are the most cited, but for the research of sensor control in ASW, they seem to be unsuitable. In this paper, digital twin depends on creating the simulation models for related real entities in the digital way to simulate their behaviors, serving as a bridge between the real and the simulated systems [24,25]. The adaptability of the digital twin-based framework to ASW is reflected in the following aspects:

• Digital twin breaks the barriers between the real ASW and the simulated ASW, and enables real-time communication in both directions.
• Digital twin can make full use of the prediction, evaluation, and analysis ability of simulation system to evaluate available courses of actions for sensor control.
• Digital twin paves a way for the cyber-physical integration of ASW, which is an important bottleneck to enable intelligent and adaptive decision making.
• Digital twin can enable the integrated application of new-generation information technologies, such as internet of things (IoT), 5G, AI, cloud, edge computing, and so on.
• Digital twin can take full advantage of POMDP and the simulation-based approach to support more complex application scenarios.

The backbone technology of digital twin is simulation for prediction, evaluation and analysis [24]. The simulation-based prediction with high confidence is the fundamental function of digital twin. The vision of the digital twin itself refers to a comprehensive functional description, evaluation, and prediction together with all available operational data of an entity, target, or system, which includes more or less all information which could be useful in all the current and subsequent phases [26]. The digital twin in ASW is not only used to describe and predict the behaviors of the submarine in the real ASW, but also to derive and evaluate solutions and course of actions (COAs) relevant to the real ASW.

The digital twin theory has been quite popular in recent years, but most of the research related to it only focuses on the theoretical research and requirement analysis. The practical application and implementation of digital twin has rarely been mentioned. This paper is the first attempt to apply and implement the digital twin theory. Simulation models for ASW are very complicated now, but still fail to describe the real ASW accurately. The main reason is the separation of simulation system and the real system, and the second reason is the failure to simultaneously integrate online measurements into the running simulation models.

In the proposed digital twin-based framework of sensor control in ASW, the online measurements are dynamically assimilated into the running simulation models, and the running simulation models guide the real ASW process in reverse. The intuitive application of digital twin for ASW is to obtain estimated states or parameters of the real ASW system by combining the real-time measurements with a simulation model [27]. Since real-time measurements can indicate the latest updated states of the real system, we focus on the problem of effectively assimilating continuous streams of data into running simulation models. At the same time, we have also studied the computation of reward function by using the results of the proposed data-assimilation algorithm.

The rest of the paper is structured as follows. We give the digital twin and RFS-based framework of sensor control in ASW in Section 2. The proposed RFS-based models for digital twin are introduced in Section 3. The RFS-based data-assimilation algorithm is detailed in Section 4. Section 5 describes the computation of the reward function by using the results of the proposed data-assimilation algorithm. Experimental results are detailed in Section 6, and the conclusions are given in Section 7.

2. Digital Twin and RFS-Based Framework of Online Sensor Control

As it is shown in Figure 1, the digital twin in ASW is mainly used for decision making. The digital twin can be regarded as a virtual equivalent or dynamic digital representation of the real ASW [28]. The simulated ASW in this framework is used to predict the emergent behaviors in the real ASW, evaluate the COAs and choose the best one for the operator. The simulated ASW evolves with the real ASW along the whole life cycle and integrates the currently available and commonly required data and knowledge. We can get the prediction, evaluation, and analysis of an enemy submarine by means of precise simulations. Digital twin can assist in ensuring information continuity throughout the whole operation, sensor control, and system behavior predictions in ASW based on simulations.

To improve the coordination between the simulation system and the real system for sensor control in ASW, we propose the digital twin-based framework of online sensor control by incorporating RFS theory. The technical view of the proposed framework is given in Figure 2. Since there are two constituent objects (simulated space and physical space) in digital twin, in this paper, we propose two corresponding technologies to support the implementation of the digital twin-based framework: one is the RFS-based data-assimilation algorithm for assimilating real-time measurements into the running simulation model, and the other one is the computation of the reward function by using the results of the proposed data-assimilation algorithm for finding the optimal control action.

Just as it is shown in Figure 2, the proposed feed-back control loop incorporates real-time measurements into the running simulation model while dynamically managing the physical sensors to refine measurements. The physical ASW space provides the simulation inputs to the simulation model, and it also provides the real states to the physical sensors.

The RFS-based simulation model provides the predicted states to the RFS-based measurement model. The simulation model helps to analyze potential alternative solutions for the anti-submarine ship, and evaluate the impact of possible control actions for the online sensor control method. The RFS-based measurement model characterizes the behaviors of the physical sensor on the anti-submarine ship. It uses the predicted states outputted by the RFS-based simulation model to generate the predicted measurements and provides them to the data-assimilation process. The RFS-based simulation model and measurement model are the main elements of the simulated ASW space.

The digital twin is supported by data-assimilation algorithm which can incorporate real-time measurements into the running simulation models for more accurate prediction of the physical ASW system, and it can also evaluate the COAs [29]. RFS-based data-assimilation process is the foundation of digital twin and is in charge of fusing real-time measurements to estimate the states of the physical ASW space [30]. There are two tasks for the data-assimilation process, the first is to dynamically update the current simulation states of the physical ASW space and provide the updated states to the simulation model for subsequent simulation running; the second is to provide the updated states and corresponding weights to the reward function computation module for computing the reward function.

3. RFS-Based Modeling of the Simulated ASW

The simulated ASW in digital twin depends on two kinds of models to describe the physical ASW: one is the Markov transition density-based simulation model; the other one is the measurement likelihood-based measurement model. Here we use the RFS-based simulation model to model the state transition of the enemy submarine, and the RFS-based measurement model to model the physical sensor of anti-submarine ship.

3.1. RFS-Based Data Model

Conventional estimation techniques fail to support digital twin for ASW, because many sophisticated simulation models of ASW cannot provide the analytical mathematical structures for deriving the functional forms of probability distribution. The sequential Monte Carlo (SMC) method which is also named as particle filter (PF), is the most widely used data-assimilation algorithm in recent years [29]. PF-based data-assimilation algorithm in traffic simulation is presented by Xu in [31] and Wu in [32]. PF-based data-assimilation algorithm in wildfire simulation is presented by Hu in [33,34]. They all use the same non-parametric statistic inference method based on PF, because PF has no assumptions on the distribution and linearity of the simulation model.

The conventional PF-based data-assimilation algorithm depends on the vector-based representation of data including states and measurements. The vector-based representation makes the vector-based data-assimilation algorithm have the following essential disadvantages:

• It is based on the assumption that the studied system is a single dynamic system that is permanently active. It cannot be used for the dynamic system that switches on and off randomly. Switching is quite common for submarine activity, for example, a submarine may enter and leave a battle area at random instance.
• It is based on the assumption that the detection is perfect with no false detections and no missed detections, and it also needs the number and ordering of measurements to be previously designated. Furthermore, it cannot jointly estimate the number of submarines and the states of each submarine.

Being different from the conventional vector-based representation, the RFS-based representation can take the more complex situations into consideration. The RFS-based representation of states can support the transition of submarines from one mode to another. RFS-based representation of states enables the submarine number to be constantly varying. For example, completely new submarines can enter a scene randomly. Submarines can likewise leave a scene, as when disappearing behind some other occlusion, or they can be damaged or destroyed. The RFS-based representation of measurements enables us to take the imperfection of the sensors into consideration. The sensors on the anti-submarine ship can fail to generate a measurement of the submarine state, or pick up false measurements. Since there is no information on which state generates the measurement, and the number of measurements is a random variable, RFS rather than vector-based representation can be more useful.

3.2. RFS-Based Measurement Model

The measurement model is in charge of mapping the simulation states to the measurements collected by the sensor. The sensor on the anti-submarine ship provides measurements with imperfect enemy submarine detection including noises, clutters, and missed detections. Detection uncertainty and clutters in the measurements can be described by the union of Bernoulli RFS and Poisson RFS as introduced in [35,36]. The sensor sequentially gives an unordered finite set of measurements ${\mathbf{Z}}_k\subset \mathcal{Z}$, and never indicates which of these measurements is generated from the enemy submarine.

The measurements produced by the sensor can be mathematically modeled as an RFS $\mathbf{Z}=\{{\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_m\}$. Its major advantage is that both the measurements’ number $m=\left|\mathbf{Z}\right|$ and the value of the constituent vector $\mathbf{z}\in \mathbf{Z}$ in the measurement space $\mathcal{Z}\subseteq {\mathbb{R}}^{n_z}$ are random. It makes no assumptions on the order of detections in the RFS $\mathbf{Z}$. The measurements for ASW can be represented as $\mathbf{Z}=\mathbf{C}\cup \mathbf{W}$, here $\mathbf{C}$ is the Poisson random finite subset of false detections and $\mathbf{W}$ is the Bernoulli random finite subset corresponding to the enemy submarine.

The behavior of the sensor can be described by the conventional likelihood function $g\left(z\right|\mathbf{x})$ which characterizes the probability that measurement z is generated by the submarine state $\mathbf{x}$. Here we can have $g\left(z\right|{\mathbf{x}}_k)=\mathcal{N}(z;h\left({\mathbf{x}}_k\right),{\omega }_k)$. If a measurement $z\in {\mathbf{Z}}_k$ is generated from an enemy submarine, then the relationship between $z\in {\mathbf{Z}}_k$ and the submarine state ${\mathbf{x}}_k$ can be described by following equation:

$$z=h\left({\mathbf{x}}_k\right)+{\omega }_k$$

(1)

where ${\omega }_k$ is a zero-mean independent Gaussian noise with variance ${\sigma }_{\theta }^2$, and $h\left({\mathbf{x}}_k\right)$ is in charge of modeling the true submarine bearing or range at time k. In the following part of this section, we will derive the specific mathematical form for the measurement model denoted by $\phi \left(\mathbf{Z}\right|\mathbf{X})$, where $\mathbf{X}$ is the Bernoulli RFS of simulation states. The measurement model $\phi \left(\mathbf{Z}\right|\mathbf{X})$ can be written in two specific forms, one for $\mathbf{X}=\varnothing$ and the other for $\mathbf{X}=\left\{\mathbf{x}\right\}$.

Firstly, if the enemy submarine does not exist in the real operational area, the measurement set will only consist of clutters. It means that $\mathbf{X}=\varnothing$ and $\mathbf{Z}=\mathbf{C}\cup \varnothing =\mathbf{C}$. Here we model the number of clutters in the measurement set by the following Poisson distribution:

$$P\left\{\right|\mathbf{C}|=s\}=\frac{e^{-\lambda }{\lambda }^s}{s!},s=0,1,2,\dots$$

(2)

here $\lambda$ is the average number of clutters. Clutters are modeled as independent identically distributed random vectors conditioned on $\left|\mathbf{C}\right|$. The values of these random vectors are taken from the measurement space $\mathcal{Z}$ with probability density function (PDF) $c\left(z\right)$. Here $c\left(z\right)={\left(2\pi \right)}^{-1}$ is the time invariant spatial distribution of clutters. The measurement model can be modeled as follows:

$$\phi \left(\mathbf{Z}\right|\varnothing )=\kappa \left(\mathbf{Z}\right)=e^{-\lambda }\prod _{z\in \mathbf{Z}}\lambda c\left(z\right)$$

(3)

Secondly, if the enemy submarine exists in the real operational area with state $\mathbf{X}=\left\{\mathbf{x}\right\}$. The measurement RFS $\mathbf{W}$ corresponding to the enemy submarine can be modeled as a Bernoulli RFS conditioned on $\mathbf{X}=\left\{\mathbf{x}\right\}$. If the submarine fails to be detected, we have $\mathbf{W}=\varnothing$. If the submarine is detected and causes a measurement z, we have $\mathbf{W}=\left\{z\right\}$. To derive the measurement model $\phi \left(\mathbf{Z}\right|\left\{\mathbf{x}\right\})$, we should specify the PDF of the RFS $\mathbf{W}$ conditioned on $\left\{\mathbf{x}\right\}$. This PDF can be modeled as follows:

$$\eta \left(\mathbf{W}\right|\left\{\mathbf{x}\right\})=\left\{\begin{array}{cc}1-p_D\left(\mathbf{x}\right)\hfill & \mathrm{if}\mathbf{W}=\varnothing \hfill \\ p_D\left(\mathbf{x}\right)·g\left(z\right|\mathbf{x}),\hfill & \mathrm{if}\mathbf{W}=\left\{z\right\}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(4)

where $p_D\left(\mathbf{x}\right)$ is the probability of detecting the submarine state $\mathbf{x}$.

The measurement model $\phi \left(\mathbf{Z}\right|\left\{\mathbf{x}\right\})$ can be represented as $\phi \left(\mathbf{Z}\right|\left\{\mathbf{x}\right\})={\sum }_{\mathbf{W}\subseteq \mathbf{Z}}\eta \left(\mathbf{W}\right|\left\{\mathbf{x}\right\})\kappa (\mathbf{Z}\setminus \mathbf{W})$, here ∖ is the set-difference operation and $\kappa$ was defined by (3). Since the RFS $\mathbf{W}$ is a Bernoulli RFS, and the summation is computable, we can get the measurement model for the case $\mathbf{X}=\left\{\mathbf{x}\right\}$ as follows:

$$\begin{array}{cc}\hfill \phi \left(\mathbf{Z}\right|\left\{\mathbf{x}\right\})& =\eta (\varnothing |\left\{\mathbf{x}\right\})·\kappa \left(\mathbf{Z}\right)+\sum _{z\in \mathbf{Z}}\eta \left(\left\{z\right\}\right|\left\{\mathbf{x}\right\})·\kappa (\mathbf{Z}\setminus \left\{z\right\})\hfill \end{array}$$

(5)

$$\begin{array}{c}\hfill =\kappa \left(\mathbf{Z}\right)\left[1-p_D\left(\mathbf{x}\right)+p_D\left(\mathbf{x}\right)\sum _{z\in \mathbf{Z}}g\left(z\right|\mathbf{x})\frac{\kappa (\mathbf{Z}\setminus \{z\left\}\right)}{\kappa \left(\mathbf{Z}\right)}\right]\end{array}$$

(6)

In summary, the RFS-based measurement model for ASW can be represented as follows:

$$\phi \left(\mathbf{Z}\right|\mathbf{X})=\left\{\begin{array}{cc}\kappa \left(\mathbf{Z}\right),\hfill & \mathrm{if}\mathbf{X}=\varnothing \hfill \\ \kappa \left(\mathbf{Z}\right)\left[1-p_D\left(\mathbf{x}\right)+p_D\left(\mathbf{x}\right)\sum _{z\in \mathbf{Z}}g\left(z\right|\mathbf{x})\frac{\kappa (\mathbf{Z}\setminus \{z\left\}\right)}{\kappa \left(\mathbf{Z}\right)}\right],\hfill & \mathrm{if}\mathbf{X}=\left\{\mathbf{x}\right\}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(7)

3.3. RFS-Based Simulation Model

Digital twin requires a simulation model that can predict possible future states. To give a unified description of the enemy submarine presence/absence in the operational area and its kinematic state, we employ the Bernoulli RFS-based simulation model to describe the dynamics of enemy submarine at discrete time k. The state space is $\varnothing \cup \mathcal{S}\left(\mathcal{X}\right)$, where $\mathcal{S}\left(\mathcal{X}\right)$ is a set of all singleton $\left\{\mathbf{x}\right\}$ and $\mathbf{x}\in \mathcal{X}$. Here we use PDF to denote the uncertainty of the submarine’s states which are evolving according to the state space model in discrete time as follows:

$${\mathbf{X}}_k=f_{k|k-1}\left({\mathbf{X}}_{k-1}\right)+{\mathbf{v}}_k$$

(8)

here $f_{k|k-1}$ denotes the deterministic part of the true evolution equation which is in charge of mapping the state to the next time step, and $\mathbf{v}$ is the stochastic part of the true evolution that we fail to capture deterministically, and it makes assimilating online measurements necessary.

If the enemy submarine is existing in the real operational area, the state ${\mathbf{X}}_k$ will be a singleton, and can be modeled by the Markov process whose Markov transition density is denoted by ${\pi }_{k|k-1}\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$ during the simulation interval ${\mathbf{T}}_k=t_k-t_{k-1}$. The submarine state vector is adopted as follows:

$${\mathbf{x}}_k^t={\left[x_k^t{\dot{x}}_k^t{y}_k^t{\dot{y}}_k^t\right]}^T$$

(9)

where $(x_k^t,y_k^t)$ and $({\dot{x}}_k^t,{\dot{y}}_k^t)$ is respectively the enemy submarine position and velocity in Cartesian coordinate. The state vector of the anti-submarine ship ${\mathbf{x}}_k^o$ is similarly represented by ${\mathbf{x}}_k^o={\left[x_k^o{\dot{x}}_k^o{y}_k^o{\dot{y}}_k^o\right]}^T$. Here ${\pi }_{k|k-1}\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$ is the nearly constant velocity (NCV) model and relies on the relative state vector which is defined as follows:

$${\mathbf{x}}_k={\left[x_k{\dot{x}}_k{y}_k{\dot{y}}_k\right]}^T={\left[x_k^t-x_k^o{\dot{x}}_k^t-{\dot{x}}_k^o{y}_k^t-y_k^o{\dot{y}}_k^t-{\dot{y}}_k^o\right]}^T.$$

(10)

The specific mathematical form of ${\pi }_{k|k-1}\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$ is described by the state space model as follows:

$${\mathbf{x}}_{k+1}={\mathbf{F}}_k{\mathbf{x}}_k-{\mathbf{S}}_{k+1,k}+{\epsilon }_k$$

(11)

here ${\mathbf{F}}_k$ is the deterministic part of the true evolution equation, ${\mathbf{S}}_{k+1,k}$ is a deterministic matrix related to the effect of anti-submarine ship acceleration, and ${\epsilon }_k\sim \mathcal{N}(0,{\mathbf{Q}}_k)$ is stochastic part of the true evolution equation. Here we have ${\pi }_{k|k-1}\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})=\mathcal{N}(\mathbf{x};{\mathbf{F}}_k{\mathbf{x}}_{k-1}-{\mathbf{S}}_{k+1,k},{\mathbf{Q}}_k)$. Their specific forms are as follows:

$${\mathbf{F}}_k={\mathbf{I}}_2\otimes \left[\begin{array}{cc}1& T_k\\ 0& 1\end{array}\right],{\mathbf{Q}}_k={\mathbf{I}}_2\otimes \varpi \left[\begin{array}{cc}\frac{T_k^3}{3}& \frac{T_k^2}{2}\\ \frac{T_k^2}{2}& T_k\end{array}\right],{\mathbf{S}}_{k+1,k}=\left[\begin{array}{c}{x}_{k+1}^o-x_k^o-T_k{\dot{x}}_k^o\\ {\dot{x}}_{k+1}^o-{\dot{x}}_k^o\\ y_{k+1}^o-y_k^o-T_k{\dot{y}}_k^o\\ {\dot{y}}_{k+1}^o-{\dot{y}}_k^o\end{array}\right]$$

(12)

where ⊗ is the Kronecker product, and $\varpi$ is the intensity of process noise. Here we assume $T_k=T=const$, so we can get ${\mathbf{F}}_k=\mathbf{F}$ and ${\mathbf{Q}}_k=\mathbf{Q}$.

To characterize the enemy submarine’s appearance and disappearance during the simulation interval, we use a binary random variable ${\xi }_k\in \{0,1\}$ to denote its existence. Here ${\xi }_k=1$ indicates that the enemy submarine exists at time k, and ${\xi }_k=0$ indicates that it does not exist at time k. The dynamics of ${\xi }_k$ is described by the first-order 2-state Markov chain with a transitional probability matrix $\mathsf{\Xi }$ as it is shown in Figure 3. The elements of $\mathsf{\Xi }$ are defined as ${\left[\mathsf{\Xi }\right]}_{ij}=P\{{\xi }_k=j-1|{\xi }_{k-1}=i-1\}$ for $i,j\in \{1,2\}$. $\mathsf{\Xi }$ is defined as follows:

$$\mathsf{\Xi }=\left[\begin{array}{cc}(1-p_b)& p_b\\ (1-p_s)& p_s\end{array}\right]$$

(13)

where $p_b=P\{{\xi }_{k+1}=1|{\xi }_k=0\}$ is the probability that the submarine appears in the operational area during the simulation interval, and $p_s=P\{{\xi }_{k+1}=1|{\xi }_k=1\}$ is the probability that the submarine is still in the operational area during the simulation interval. Here $p_b$ and $p_s$ are assumed to be known. If the submarine appears during the simulation interval $T_k$, PDF $b_{k|k-1}\left(\mathbf{x}\right)$ can be used to denote its PDF.

Now we derive the RFS-based simulation model of enemy submarine. Since we represent the submarine’s simulation states by using Bernoulli RFS, we consider two different situations for the RFS-based simulation model.

Firstly, if the enemy submarine does not exist in the real operational area at time k, we can get ${\mathbf{X}}_k=\varnothing$. In addition, the submarine can appear in the operational area with probability $p_{b,k}$ and have kinematic state ${\mathbf{x}}_k$ with PDF $b_{k|k-1}\left({\mathbf{x}}_k\right)$, or remain absent from the operational area with probability $1-p_{b,k}$. The simulation model $f_{k|k-1}\left({\mathbf{X}}_k\right|\varnothing )$ for state ${\mathbf{X}}_k$ is specified as follows:

$$f_{k|k-1}\left({\mathbf{X}}_k\right|\varnothing )=\left\{\begin{array}{cc}1-p_{b,k},\hfill & \mathrm{if}{\mathbf{X}}_k=\varnothing \hfill \\ p_{b,k}·b_{k|k-1}\left({\mathbf{x}}_k\right),\hfill & \mathrm{if}{\mathbf{X}}_k=\left\{{\mathbf{x}}_k\right\}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(14)

Secondly, if the enemy submarine is existing in the real operational area at time $k-1$, which means ${\mathbf{X}}_{k-1}=\left\{{\mathbf{x}}_{k-1}\right\}$, it can come through to the next time step with probability $p_{s,k}\left({\mathbf{x}}_{k-1}\right)$ and transit to ${\mathbf{x}}_k$ with PDF ${\pi }_{k|k-1}\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$, or disappear with probability $1-p_{s,k}\left({\mathbf{x}}_{k-1}\right)$. Thus, the simulation model for state ${\mathbf{X}}_k$ at time step k is characterized by

$$f_{k|k-1}\left({\mathbf{X}}_k\right|\left\{{\mathbf{x}}_{k-1}\right\})=\left\{\begin{array}{cc}1-p_{s,k}\left({\mathbf{x}}_{k-1}\right),\hfill & \mathrm{if}{\mathbf{X}}_k=\varnothing \hfill \\ p_{s,k}\left({\mathbf{x}}_{k-1}\right)·{\pi }_{k|k-1}\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1}),\hfill & \mathrm{if}{\mathbf{X}}_k=\left\{{\mathbf{x}}_k\right\}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(15)

From (14) and (15), we can know that the Bernoulli RFS can either be an empty set or a nonempty set with one element only. In addition, the probabilities of these two cases can be modeled respectively as $1-q$ and q. If it has one element, the element will be spatially distributed over $\mathcal{S}\left(\mathcal{X}\right)\subseteq {\mathbb{R}}^{n_x}$ according to the standard PDF $s\left(\mathbf{x}\right)$. So, the simulation model for the simulation state denoted by the Bernoulli RFS ${\mathbf{X}}_k$ at time k can be given by:

$$f\left({\mathbf{X}}_k\right)=\left\{\begin{array}{cc}1-q_{k|k},\hfill & \mathrm{if}{\mathbf{X}}_k=\varnothing \hfill \\ q_{k|k}·s_{k|k}\left({\mathbf{x}}_k\right),\hfill & \mathrm{if}{\mathbf{X}}_k=\left\{{\mathbf{x}}_k\right\}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(16)

where $q_{k|k}=p\left\{\right|{\mathbf{X}}_k|=1|{\mathbf{Z}}_{1:k}\}$ is the probability of submarine’s existence in the real operational area, $s_{k|k}\left({\mathbf{x}}_k\right)=p\left({\mathbf{x}}_k\right|{\mathbf{Z}}_{1:k})$ is the spatial PDF of the submarine.

4. RFS-Based Data-Assimilation Algorithm

Unforeseen operational entities can enter the designated real operational area and critical assets can be destroyed, thereby invalidating the current simulation setting [37]. It is possible that the simulation execution may not proceed as expected because of the dynamic and changing environment. Since assimilating real-time measurements into the running simulation models can significantly improve the accuracy of the simulation results, the implementation of digital twin theory depends on the simulation system’s use of the information from the real world.

Here an RFS-based data-assimilation algorithm is proposed for online assimilating the sequence of measurement sets in the present of noise, false alarms, data association uncertainty, and detection uncertainty. The proposed RFS-based data-assimilation algorithm can overcome the limitations of the standard vector-based data-assimilation algorithms very well.

4.1. Data Assimilation with RFS-Based Models

Since data assimilation usually updates the posterior distribution of simulation states by using the simulation model and measurement model, it can be regarded as a Bayesian inference procedure from a probabilistic point of view [38]. Here we also use the RFS-based models and Bayesian inference for data assimilation.

The RFS-based data assimilation estimates recursively the posterior PDF of submarine’s states by using the RFS-based simulation model, measurement model, and online measurements. It usually consists of two stages, the prediction stage and the update stage. Since we have got a Bernoulli RFS-based simulation model (16), the posterior PDF at time k, denoted by $f_{k|k}\left({\mathbf{X}}_k\right|{\mathbf{Z}}_{1:k})$ is completely specified by two terms:

• the posterior submarine’s existence probability $q_{k|k}=p\left\{\right|{\mathbf{X}}_k|=1|{\mathbf{Z}}_{1:k}\}$;
• the posterior spatial PDF of ${\mathbf{X}}_k=\left\{\mathbf{x}\right\}$ denoted by $s_{k|k}\left({\mathbf{x}}_k\right)=p\left({\mathbf{x}}_k\right|{\mathbf{Z}}_{1:k})$.

Denote the posterior PDF of simulation states at time $k-1$ as $f_{k-1|k-1}\left({\mathbf{X}}_{k-1}\right|{\mathbf{Z}}_{1:k-1})$. The RFS-based prediction equation of the data-assimilation process is:

$$\begin{array}{cc}\hfill f_{k|k-1}\left({\mathbf{X}}_k\right|{\mathbf{Z}}_{1:k-1})& =\int f_{k|k-1}\left({\mathbf{X}}_k\right|{\mathbf{X}}^{{}^{\prime }})f_{k-1|k-1}\left({\mathbf{X}}^{{}^{\prime }}\right|{\mathbf{Z}}_{1:k-1})\delta {\mathbf{X}}^{{}^{\prime }}\hfill \end{array}$$

(17)

$$\begin{array}{c}=f_{k|k-1}\left({\mathbf{X}}_k\right|\varnothing )f_{k-1|k-1}(\varnothing |{\mathbf{Z}}_{1:k-1})+\int f_{k-1|k-1}\left(\left\{{\mathbf{x}}^{{}^{\prime }}\right\}\right|{\mathbf{Z}}_{1:k-1})d{\mathbf{x}}^{{}^{\prime }}.\hfill \end{array}$$

(18)

At time $k-1$, the posterior PDF is characterized by the pair $(q_{k-1|k-1},s_{k-1|k-1}\left(\mathbf{x}\right))$, the prediction equation of data assimilation can be derived from (18) as:

$$q_{k|k-1}=p_b(1-q_{k-1|k-1})+p_s{q}_{k-1|k-1}$$

(19)

$$s_{k|k-1}\left(\mathbf{x}\right)=\frac{p_b(1-q_{k-1|k-1})b_{k|k-1}\left(\mathbf{x}\right)+p_s{q}_{k-1|k-1}\int {\pi }_{k|k-1}\left(\mathbf{x}\right|{\mathbf{x}}^{{}^{\prime }})s_{k-1|k-1}\left({\mathbf{x}}^{{}^{\prime }}\right)d{\mathbf{x}}^{{}^{\prime }}}{p_b(1-q_{k-1|k-1})+p_s{q}_{k-1|k-1}}$$

(20)

In Bayesian theory, the updated PDF is calculated by:

$$f_{k|k}\left({\mathbf{X}}_k\right|{\mathbf{Z}}_{1:k})=\frac{{\phi }_k\left({\mathbf{Z}}_k\right|{\mathbf{X}}_k)f_{k|k-1}\left({\mathbf{X}}_k\right|{\mathbf{Z}}_{1:k-1})}{\int {\phi }_k\left({\mathbf{Z}}_k\right|\mathbf{X})f_{k|k-1}\left(\mathbf{X}\right|{\mathbf{Z}}_{1:k-1})\delta \mathbf{X}}=\frac{{\phi }_k\left({\mathbf{Z}}_k\right|{\mathbf{X}}_k)f_{k|k-1}\left({\mathbf{X}}_k\right|{\mathbf{Z}}_{1:k-1})}{f_k\left({\mathbf{Z}}_k\right|{\mathbf{Z}}_{1:k-1})}$$

(21)

Given the measurement model (7) and the prediction Equations (19) and (20), the update equation can be derived from (21) as:

$$q_{k|k}=\frac{1-{\delta }_k}{1-{\delta }_k·q_{k|k-1}}·q_{k|k-1}$$

(22)

$$s_{k|k}\left(\mathbf{x}\right)=\frac{1-p_D+p_D{\sum }_{z\in {\mathbf{Z}}_k}\frac{g_k\left(z\right|\mathbf{x})}{\lambda c\left(z\right)}}{1-{\delta }_k}·s_{k|k-1}\left(\mathbf{x}\right)$$

(23)

where

$${\delta }_k=p_D\left(1-\sum _{z\in {\mathbf{Z}}_k}\frac{\int g_k\left(z\right|\mathbf{x})s_{k|k-1}\left(\mathbf{x}\right)d\mathbf{x}}{\lambda c\left(z\right)}\right).$$

(24)

4.2. SMC-Based Calculation

In the general cases, the RFS-based prediction equation and update equation cannot be solved analytically [39]. Here we propose a SMC-based calculation. We use the particle system ${\{w_k^{\left(i\right)},{\mathbf{x}}_k^{\left(i\right)}\}}_{i=1}^N$ to approximate the spatial PDF $s_{k|k}\left(\mathbf{x}\right)$. Here ${\mathbf{x}}_k^{\left(i\right)}$ is the state of particle i and $w_k^{\left(i\right)}$ is its weight. As $s_{k|k}\left(\mathbf{x}\right)$ is a conventional PDF, the corresponding weights of the particles should be normalized, i.e., ${\sum }_{i=1}^N{w}_k^{\left(i\right)}=1$.

As it is shown in Figure 4, the RFS-based simulation model firstly runs to the next time point and generates the predicted states of ASW. Weight updating relies on the difference between the online measurements and the predicted measurements generated from the predicted states. Suppose at time $k-1$, the submarine’s existence probability is $q_{k-1|k-1}$, and the spatial PDF is approximated by

$${\widehat{s}}_{k-1|k-1}\left(\mathbf{x}\right)=\sum _{i=1}^N{w}_{k-1}^{\left(i\right)}{\delta }_{{\mathbf{x}}_{k-1}^{\left(i\right)}}\left(\mathbf{x}\right).$$

(25)

here ${\delta }_{\mathbf{c}}\left(\mathbf{x}\right)$ is the Dirac delta function concentrated at state $\mathbf{c}$. The predicted submarine’s existence probability $q_{k|k-1}$ can be computed by (19). According to (20) and (25), the predicted spatial PDF depends on the sum of two terms. Consequently $s_{k|k-1}\left(\mathbf{x}\right)$ can be approximated by the particle system as a weighted point mass representation:

$${\widehat{s}}_{k|k-1}\left(\mathbf{x}\right)=\sum _{i=1}^{N+B_k}{w}_{k|k-1}^{\left(i\right)}{\delta }_{{\mathbf{x}}_{k|k-1}^{\left(i\right)}}\left(\mathbf{x}\right).$$

(26)

here the particles are drawn from two importance densities, ${\rho }_k$ for persisting particles and ${\beta }_k$ for birth particles:

$${\mathbf{x}}_{k|k-1}^{\left(i\right)}\sim \left\{\begin{array}{cc}{\rho }_k\left(\mathbf{x}\right|{\mathbf{x}}_{k-1}^{\left(i\right)},{\mathbf{Z}}_k)\hfill & i=1,\dots ,N\hfill \\ {\beta }_k\left(\mathbf{x}\right|{\mathbf{Z}}_k)\hfill & i=N+1,\dots ,N+B_k\hfill \end{array}\right.$$

(27)

with weights

$$w_{k|k-1}^{\left(i\right)}=\left\{\begin{array}{cc}\frac{p_s·q_{k-1|k-1}}{q_{k|k-1}}\frac{{\pi }_{k|k-1}\left({\mathbf{x}}_k^{\left(i\right)}\right|{\mathbf{x}}_{k-1}^{\left(i\right)})}{{\rho }_k\left({\mathbf{x}}_{k|k-1}^{\left(i\right)}|{\mathbf{x}}_{k-1}^{\left(i\right)},{\mathbf{Z}}_k\right)}{w}_{k-1}^{\left(i\right)},\hfill & i=1,\dots ,N\hfill \\ \frac{p_b·(1-q_{k-1|k-1})}{q_{k|k-1}}\frac{b_{k|k-1}\left({\mathbf{x}}_{k|k-1}^{\left(i\right)}\right)}{{\beta }_k\left({\mathbf{x}}_{k|k-1}^{\left(i\right)}\right|{\mathbf{Z}}_k)}\frac{1}{B_k},\hfill & i=N+1,\dots ,N+B_k.\hfill \end{array}\right.$$

(28)

where $B_k$ is the number of submarine-birth particles drawn from the importance density ${\beta }_k$.

The simplest choice of importance density ${\rho }_k\left(\mathbf{x}\right|{\mathbf{x}}_{k-1},{\mathbf{Z}}_k)$ is the transitional density ${\pi }_{k|k-1}\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1})$. If there is little prior knowledge of the action plan of enemy submarine, we should assume that the enemy submarine can appear anywhere in the state space $\mathcal{S}\left(\mathcal{X}\right)$. So, we model $b_{k|k-1}\left(\mathbf{x}\right)$ by using the uniform distribution over $\mathcal{S}\left(\mathcal{X}\right)$. The birth importance density ${\beta }_k$ in (27) needs to have the same support as $b_{k|k-1}\left(\mathbf{x}\right)$ (i.e., the entire $\mathcal{S}\left(\mathcal{X}\right)$) [40].

The update step of the data-assimilation algorithm is implemented according to (22)–(24). First, for every $z\in {\mathbf{Z}}_k$, the integral $I_k\left(z\right)=\int g_k\left(z\right|\mathbf{x})s_{k|k-1}\left(\mathbf{x}\right)d\mathbf{x}$, which appears in (24), is approximately calculated as follows:

$$I_k\left(z\right)\approx \sum _{i=1}^{N+B_k}{w}_{k|k-1}^{\left(i\right)}·g_k\left(z\right|{\mathbf{x}}_{k|k-1}^{\left(i\right)})$$

(29)

Then based on (24), ${\delta }_k$ can be calculated by

$${\delta }_k\approx p_D\left(1-\sum _{z\in {\mathbf{Z}}_k}\frac{I_k\left(z\right)}{\lambda c\left(z\right)}\right).$$

(30)

The submarine’s existence probability is updated by using (22), and the corresponding weights are updated according to (23):

$${\tilde{w}}_{k|k}^{\left(i\right)}\approx \left[1-p_D+p_D\sum _{z\in {\mathbf{Z}}_k}\frac{g_k\left(z\right|{\mathbf{x}}_{k|k-1}^{\left(i\right)})}{\lambda c\left(z\right)}\right]·w_{k|k-1}^{\left(i\right)}$$

(31)

These weights should be normalized to get the normalized importance weights:

$$w_{k|k}^{\left(i\right)}=\frac{{\tilde{w}}_{k|k}^{\left(i\right)}}{{\sum }_{j=1}^{N+B_k}{\tilde{w}}_{k|k}^{\left(j\right)}}$$

(32)

for $i=1,\dots ,N+B_k$. At last, we resample N times from ${\{w_{k|k}^{\left(i\right)},{\mathbf{x}}_{k|k-1}^{\left(i\right)}\}}_{i=1}^{N+B_k}$ to avoid sample degeneracy.

5. Computation of Reward Function

We propose an online sensor control method by using the predicted states, updated states, and their corresponding weights generated by the proposed data-assimilation algorithm. Online sensor control is mainly used for finding the optimal control action from a set of admissible control actions. Thus, it means sequential decision making, where each decision generates measurements that provide an additional information for data assimilation. In the digital twin for ASW, the control action is determined in the present of uncertainty both in the measurement space and the state space. Here an information theoretic approach is proposed for online sensor control. In this approach, the posterior PDF is used to represent the uncertain states, and the reward function is regarded as a measure of the information gain related to each action.

5.1. Derivation of Reward Function

In this paper, the online sensor control means the online selection of headings for individual anti-submarine ship, to maximize the use efficiency of its measurement system. Here control actions are ranked by using the quantity of information predicted to be gained from their execution. The data-assimilation enhanced simulation model is used to rapidly predict the possible output information of alternative control actions.

Reward function is used to measure the reduction in the information gain, in comparison with the current information state. The information gain can be characterized by using various information measures [41]. The Fisher information is typically used as a criterion for optimization in the absence of detection uncertainty [42,43,44]. Here we use the Rényi divergence-based reward function. The Rényi information divergence provides a way to measure the dissimilarity between two probability densities [45]. The Rényi divergence between any two probability densities $p_0\left(\mathbf{x}\right)$ and $p_1\left(\mathbf{x}\right)$ is described as:

$$I_{\alpha }(p_1,p_0)=\frac{1}{\alpha -1}log\int p_1^{\alpha }\left(\mathbf{x}\right)p_0^{1-\alpha }\left(\mathbf{x}\right)d\mathbf{x}$$

(33)

where $\alpha \ge 0$ is the factor that reflects how much we emphasize the tails of two probability distributions.

Let ${\mathbf{u}}_k\in {\mathbb{U}}_k$ denote the control action chosen for controlling the sensor at time $t_k$ in order to collect the future measurements at $t_{k+1}$. Here ${\mathbb{U}}_k$ denotes the set of admissible control actions at time $t_k$. In general, both the simulation model $f_{k+1|k}$ and the measurement model ${\phi }_{k+1}$ depend on the control action ${\mathbf{u}}_k\in {\mathbb{U}}_k$. Then the prediction Equation (18) and update Equation (21) for RFS-based data assimilation can be rewritten as follows:

$$f_{k+1|k}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k},{\mathbf{u}}_{0:k})=\int f_{k+1|k}\left(\mathbf{X}\right|{\mathbf{X}}^{{}^{\prime }},{\mathbf{u}}_k)f_{k|k}\left({\mathbf{X}}^{{}^{\prime }}\right|{\mathbf{Z}}_{1:k},{\mathbf{u}}_{0:k-1})\delta {\mathbf{X}}^{{}^{\prime }}$$

(34)

$$f_{k+1|k+1}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k+1},{\mathbf{u}}_{0:k})=\frac{{\phi }_{k+1}\left({\mathbf{Z}}_{k+1}\right|{\mathbf{X}}_{k+1},{\mathbf{u}}_k)f_{k+1|k}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k},{\mathbf{u}}_{0:k})}{\int {\phi }_k\left({\mathbf{Z}}_{k+1}\right|\mathbf{X},{\mathbf{u}}_k)f_{k+1|k}\left(\mathbf{X}\right|{\mathbf{Z}}_{1:k},{\mathbf{u}}_{0:k})\delta \mathbf{X}}$$

(35)

The optimal control action to be applied at time k is defined by maximizing the expected Rényi information divergence according to equation

$${\mathbf{u}}_k=arg\underset{\mathbf{v}\in {\mathbb{U}}_k}{max}\mathbb{E}\left\{\varphi (\mathbf{v},p\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k},{\mathbf{u}}_{1:k-1}),{\mathbf{Z}}_{k+1}\left(\mathbf{v}\right))\right\}$$

(36)

where $\varphi (\mathbf{v},p,\mathbf{Z})$ is the real-valued reward function associated with the control action $\mathbf{v}$. (36) results in the predicted PDF p and the future measurement set $\mathbf{Z}$. Online sensor control via (36) tries to obtain the maximum reward based on a single future step. This is done by anticipating possible future measurements. To find the optimal control action for the anti-submarine ship to take next, we should predict the system states if the control action ${\mathbf{u}}_k$ has been chosen by using the RFS-based simulation model. In addition, we should generate the predicted measurement set before actually receiving the measurement set ${\mathbf{Z}}_{k+1}$. Hence the calculation of the expected value of the Rényi divergence for each possible control action is closely related to the data-assimilation process.

The future measurement set ${\mathbf{Z}}_{k+1}\left(\mathbf{v}\right)$ supports the computation of the reward function $\varphi$. Since ${\mathbf{Z}}_{k+1}\left(\mathbf{v}\right)$ is obtained after the control action has been executed, this will create uncertainty. To overcome the impact of uncertainty, (36) employs the expectation operator $\mathbb{E}$. The reward function $\varphi ({\mathbf{u}}_k,p,\mathbf{Z})$ in (36) is adopted as the Rényi divergence between:

• the predicted PDF $f_{k+1|k}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k},{\mathbf{u}}_{0:k})$ given by (34) which is based on action ${\mathbf{u}}_k$, and
• the predicted future posterior $f_{k+1|k+1}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k+1},{\mathbf{u}}_{0:k})$ given by (35), obtained by using the new measurement set ${\mathbf{Z}}_{k+1}$ collected after the sensor has been controlled to take action ${\mathbf{u}}_k$.

We simplify the reward function $\varphi$ by suppressing its second and third argument. Depending on (33), the reward function $\varphi$ can be represented by:

$$\varphi \left({\mathbf{u}}_k\right)=\frac{1}{\alpha -1}log\int {\left[f_{k+1|k+1}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k+1},{\mathbf{u}}_{0:k})\right]}^{\alpha }{\left[f_{k+1|k}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k},{\mathbf{u}}_{0:k})\right]}^{1-\alpha }\delta {\mathbf{X}}_{k+1}$$

(37)

5.2. Data-Assimilation-Based Computation

The expected reward function $\mathbb{E}\left[\varphi \left({\mathbf{u}}_k\right)\right]$ does not have the closed-form analytic solution. Thus, we employ numerical approximate method and data-assimilation-based solution for it. This makes the SMC-based implementation of data assimilation become valuable. By adopting the SMC-based implementation, the optimal control action selection can be implemented quickly. Furthermore, this makes the data-assimilation process and online sensor control to be interdependent and interoperable.

Since the predicted PDF and the updated PDF in the RFS-based data-assimilation process are Bernoulli PDFs, let $f_{k+1|k}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k},{\mathbf{u}}_{0:k})$ and $f_{k+1|k+1}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k+1},{\mathbf{u}}_{0:k})$, which feature in (37), be specified by the pairs ($q_{k+1|k}$, $s_{k+1|k}\left(\mathbf{x}\right)$) and ($q_{k+1|k+1}$, $s_{k+1|k+1}\left(\mathbf{x}\right)$) respectively. So, the predicted PDF and the updated PDF can be written as follows:

$$f_{k+1|k}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k},{\mathbf{u}}_{0:k})=\left\{\begin{array}{cc}1-q_{k+1|k}\hfill & \mathrm{if}{\mathbf{X}}_{k+1}=\varnothing ,\hfill \\ q_{k+1|k}·s_{k+1|k}\left(\mathbf{x}\right)\hfill & \mathrm{if}{\mathbf{X}}_{k+1}=\left\{\mathbf{x}\right\},\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(38)

$$f_{k+1|k+1}\left({\mathbf{X}}_{k+1}\right|{\mathbf{Z}}_{1:k+1},{\mathbf{u}}_{0:k})=\left\{\begin{array}{cc}1-q_{k+1|k+1}\hfill & \mathrm{if}{\mathbf{X}}_{k+1}=\varnothing ,\hfill \\ q_{k+1|k+1}·s_{k+1|k+1}\left(\mathbf{x}\right)\hfill & \mathrm{if}{\mathbf{X}}_{k+1}=\left\{\mathbf{x}\right\},\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(39)

According to the rules of set integral, the reward function defined in (37) can be simplified to:

$$\begin{array}{cc}\hfill \varphi \left({\mathbf{u}}_k\right)=& \frac{1}{\alpha -1}log\{{[1-q_{k+1|k}]}^{1-\alpha }{[1-q_{k+1|k+1}]}^{\alpha }\hfill \\ & +{\left[q_{k+1|k}\right]}^{1-\alpha }{\left[q_{k+1|k+1}\right]}^{\alpha }\int {\left[s_{k+1|k}\left(\mathbf{x}\right)\right]}^{1-\alpha }·{\left[s_{k+1|k+1}\left(\mathbf{x}\right)\right]}^{\alpha }d\mathbf{x}\}\hfill \end{array}$$

(40)

According to (36), the optimal control action can be selected as the expected value:

$${\mathbf{u}}_k=arg\underset{\mathbf{v}\in {\mathbb{U}}_k}{max}\mathbb{E}\left\{\varphi \left(\mathbf{v}\right)\right\}$$

(41)

Now we use SMC to obtain the numerical implementation of (41). First, we obtain the values of the reward functions for the predicted future measurement set sequence ${\mathbf{Z}}_{k+1}\left(\mathbf{v}\right)$. Then we compute the expected value of $\varphi \left(\mathbf{v}\right)$ by calculating the sample mean of the obtained values. Here ${\mathbf{Z}}_{k+1}\left(\mathbf{v}\right)$ is obtained after taking control action $\mathbf{v}\in {\mathbb{U}}_k$. Each realization of ${\mathbf{Z}}_{k+1}\left(\mathbf{v}\right)$ is generated from the predicted PDF represented by $\left(q_{k+1|k},s_{k+1|k}\left(\mathbf{x}\right)\right)$ by using the RFS-based measurement model.

Basing on the SMC-based implementation of the proposed data-assimilation algorithm, we select M predicted submarine states from ${\left\{w_{k+1|k}^{\left(i\right)},{\mathbf{x}}_{k+1|k}^{\left(i\right)}\right\}}_{i=1}^{N+B_k}$ with probability $w_{k+1|k}^{\left(i\right)}$. So, we can get M predict ideal measurements for the computation of the reward function. The key of computing (40) is the computation of the following integral

$$Q=\int {\left[s_{k+1|k}\left(\mathbf{x}\right)\right]}^{1-\alpha }·{\left[s_{k+1|k+1}\left(\mathbf{x}\right)\right]}^{\alpha }d\mathbf{x}$$

(42)

Depending on the SMC-based implementation of the proposed data assimilation, it can be computed as follows. Let $s_{k+1|k}\left(\mathbf{x}\right)$ be approximated by ${\{w_{k+1|k}^{\left(i\right)},{\mathbf{x}}_{k+1|k}^{\left(i\right)}\}}_{i=1}^{N+B_k}$. By taking control action $\mathbf{v}$, we can obtain a sample of the future measurement set ${\left\{{\mathbf{Z}}_{k+1}^{\left(i\right)}\left(\mathbf{v}\right)\right\}}_{i=1}^M$ with M submarine originated noiseless measurement set. According to the proposed RFS-based data-assimilation algorithm, $s_{k+1|k+1}\left(\mathbf{x}\right)$ can be approximated by the particle system ${\{w_{k+1|k+1}^{\left(i\right)},{\mathbf{x}}_{k+1|k+1}^{\left(i\right)}\}}_{i=1}^{N+B_k}$, where $w_{k+1|k+1}^{\left(i\right)}$ is computed according to (31)–(32). Equation (42) relies on the measurement set ${\left\{{\mathbf{Z}}_{k+1}^{\left(i\right)}\left(\mathbf{v}\right)\right\}}_{i=1}^M$ and $w_{k+1|k+1}^{\left(i\right)}$, and is approximated by:

$$Q\approx \sum _{i=1}^{N+B_k}{\left(w_{k+1|k}^{\left(i\right)}\right)}^{1-\alpha }·{\left(w_{k+1|k+1}^{\left(i\right)}\right)}^{\alpha }$$

(43)

In conclusion, the computation of the reward function assigned to every control action $\mathbf{v}\in {\mathbb{U}}_k$ for the RFS-based online sensor control method is as follows:

$$\begin{array}{cc}\hfill \varphi \left(\mathbf{v}\right)=& \frac{1}{\alpha -1}log\{{\left[1-q_{k+1|k}\right]}^{1-\alpha }{\left[1-q_{k+1|k+1}\right]}^{\alpha }\hfill \\ & +{\left[q_{k+1|k}\right]}^{1-\alpha }{\left[q_{k+1|k+1}\right]}^{\alpha }\sum _{i=1}^{N+B_k}{\left(w_{k+1|k}^{\left(i\right)}\right)}^{1-\alpha }·{\left(w_{k+1|k+1}^{\left(i\right)}\right)}^{\alpha }\}\hfill \end{array}$$

(44)

followed by the application of (41).

6. Simulation Experiments

To verify the effectiveness of the proposed digital twin and RFS-based approach for sensor control in ASW, we have carried out three groups of experiments: the first one is related to the data-assimilation algorithm, the second one is related to the online sensor control with single submarine, the last one is related for the online sensor control with multiple submarines. The data-assimilation experiment is used to verify the correctness and effectiveness of the proposed RFS-based data-assimilation algorithm for assimilating the online measurements to the simulation system. The following two groups of online sensor control experiments are-based the data-assimilation experiment, and used to verify the proposed sensor control method for ASW. They use the prediction results of the data-assimilation experiment to control the real sensor of the anti-submarine ship.

6.1. Data-Assimilation Experiment

As it is shown in Figure 5, we adopt the identical-twin experiment to evaluate the proposed data-assimilation algorithm [34,46]. In this experiment, the RFS-based simulation model of enemy submarine with designated initial setting is first running, and the measurements corresponding to the simulation results are recorded. These simulation results are regarded as the real states of the physical ASW. And the measurements recorded here are regarded as the real-time measurements generated by the real sensor. We assimilate the measurements by using the proposed data-assimilation algorithm, and then compare the assimilated simulation states with the obtained real states.

As it is shown in Figure 5, we use three terms in the experiment: the real submarine state, the assimilated one, and the simulated one, to present the experimental result. A real submarine state is the simulated one from which the measurements are recorded. To reflect the fact that the submarine simulation execution usually depends on the biased initial parameters as compared with the real submarine in ASW, here the simulated submarine state is the simulation result based on some biased initial parameters, for example, imprecise process noise intensity. Here “biased” means in the sense that the parameters are different from those used in the real submarine state. Finally, an assimilated submarine state is the data-assimilation enhanced simulation result based on the same biased initial parameters as in the simulated one. The goal of this experiment is to prove that the assimilated submarine state is more accurate than the simulated one by assimilating measurements.

6.1.1. Experimental Setup

In the real ASW, the submarine moves at a speed of approximately 5 knots. The scan repetition time of the sensor on the anti-submarine ship is 30 s. The probability of detection is assumed to be Gaussian distributed with mean 0 and covariance ${\sigma }_D=5000$. The number of clutters per scan is assumed to be Poisson distributed with the mean value $\lambda =1$. The parameters of the data-assimilation process are as follows: particle number $N=5000$, birth probability $p_b=0.01$. The initial parameters for the enemy submarine and the anti-submarine ship are presented in

Table 1. Parameters settings in the experiment.

Name	Parameter	True Value	Biased Value
Enemy submarine	Initial speed	5 knot	3 knot
	Process noise intensity $\varpi$	0.0	0.2
	Initial position	(10,000, 1000) m	(9500, 1000) m
	Initial heading	−135 deg	−90 deg
	Survival Probability $p_s$	0.99	-
Anti-submarine ship	Initial speed	4 knot	-
	Initial heading	−50 deg	-
	Initial position	(0, 0) m	-
	Detection probability $P_D$	$\mathcal{N}(\mathbf{x},0,5000)$	-
	Max detection probability $p_{D,Max}$	0.98	-
	Measurement standard deviation ${\sigma }_{\theta }$	1 deg	-

.

The performance measure of experimental results is the positional root mean square (RMS) error defined as follows:

$${\epsilon }_k=\sqrt{\frac{1}{P}{\sum }_{p=1}^P\left[{({\widehat{x}}_{k|k}^{\left(p\right)}-x_k)}^2+{({\widehat{y}}_{k|k}^{\left(p\right)}-y_k)}^2\right]}$$

(45)

where P is the total number of Monte Carlo runs, $({\widehat{x}}_{k|k}^{\left(p\right)},{\widehat{y}}_{k|k}^{\left(p\right)})$ is the assimilated (or simulated) submarine state at time k in the pth run, and $(x_k,y_k)$ is the ground truth.

6.1.2. Experimental Results

Figure 6a displays the real submarine trajectory, simulated one, and assimilated one by averaging over 500 Monte Carlo runs. Figure 6b displays the RMS error curves by averaging 500 Monte Carlo runs. This experiment tests the effectiveness of the proposed data-assimilation algorithm when the process noise intensity, initial speed, heading, and position are biased. From Figure 6 we can see that the simulated one has large deviations from the real submarine state because of the erroneous initial parameters. However, the assimilated submarine state is much closer to the real one. By using the proposed data-assimilation algorithm, the assimilated submarine state overcomes the problem of erroneous initial parameters, and matches the real submarine state with much smaller errors.

Figure 7a illustrates a typical result of a single run of the submarine’s existence probability obtained by the data-assimilation algorithm with $\varpi =0.2$. The red dotted line $q_{k|k}=1$ is the ground truth of the submarine’s existence probability which means that the enemy submarine exists in the real ASW all the way. The submarine’s existence probability shown in Figure 7a grows to 1 after some time steps and remains high throughout the simulation execution. If the detection of the submarine is missing, it drops but is still bigger than 0.8. Figure 7b shows the submarine’s existence probability averaged over 500 Monte Carlo simulations. We can find that as time involves, the assimilated submarine’s existence probability gradually approaches to 1. The occasional missed detections and false detections could not affect markedly the performance of the data-assimilation algorithm for this application.

The results of data-assimilation experiment prove that the proposed data-assimilation algorithm can successfully assimilate the online measurements to the running simulation model of ASW. In the following section, we will analyze the sensitivity of the proposed data-assimilation algorithm.

6.1.3. Sensitivity Analysis

The influence of particle number N on the overall performance of the proposed data-assimilation algorithm is studied by using different particle numbers. Figure 8a shows the RMS position errors averaged over 500 Monte Carlo simulations for different particle number N. From Figure 8a we can find that if the particle number N increases, the RMS position error will decrease. However, if the particle number N is larger than a certain degree, the influence of the particle number N on the RMS position error will be every small. This is consistent with standard particle filter theory. Figure 8b shows the influence of the particle number N on the estimated probability of submarine existence $q_{k|k}$ for $k=1,2,\dots ,80$. From Figure 8b we can find that the influence of particle number N on the estimated probability of submarine existence is very limited for this application.

To estimate the influence of the mean value $\lambda$ of clutters on the performance of the data-assimilation algorithm, Figure 9a,b shows the RMS position errors and assimilated probability of submarine existence curves obtained from different $\lambda$. From Figure 9a we can find that if $\lambda$ is small than a certain value, its influence on RMS position error is very limited. However, if $\lambda$ increases, the RMS position error will gradually increase, too. Figure 9b shows that if the mean value $\lambda$ of clutters increases, the error of estimated probability of submarine existence will also increase. In addition, it will take more time for assimilated probability of submarine existence to approach the ground truth. The results agree with Equation (22), (23) and (31). When $\lambda$ increases, the updated weight w for the particles representing the true submarine states will decrease, this leads to the increasing of RMS position error.

We also compare the performance of the data-assimilation algorithm for the different settings of max detection probability $p_{D,Max}$. The results are as shown in Figure 10a,b. From these figures, we can find that the bigger $p_{D,Max}$ is, the smaller the errors of RMS position and probability of submarine existence are. If $p_{D,Max}$ increases, according to (31), the updated weights for the submarine generated measurements will increase. This also makes the accuracy of data-assimilation increase.

The influence of $p_s$ on the performance of the data-assimilation algorithm is shown in Figure 11a,b. From Figure 11a,b we can find that the increase of $p_s$ will improve the accuracy of estimated submarine states and probability of submarine existence. From Equation (19), we know that if $p_s$ increases, the predicted existence probability q will also increase. From Equation (28), we know that if $p_s$ increases, the predicted weights for the survival particles will increase, and this leads to the decreasing of the RMS position error.

6.2. Online Sensor Control Experiment with Single Submarine

To verify the effectiveness of the proposed online sensor control method, we use a scenario where the anti-submarine ship trajectory consists of only two constant velocity motion legs. The online sensor control is conducted at the end of the first leg, when the choice is between different turns at different headings. In this experiment, there is only one enemy submarine in the operational area.

6.2.1. Experimental Setup

During the first leg, the speed of the anti-submarine ship is 4 knots and the course during this leg is $-{50}^{\circ }$. As it is shown in Figure 12, at the end of the first leg ($k=50$), the anti-submarine ship needs to choose a new course for the second leg. We verify the performance of the proposed online sensor control method by using detection parameters $p_{D,Max}=0.98$, ${\sigma }_D=5000$ and $\lambda =5$. To find the best option for the anti-submarine ship heading among the 24 second leg options, we need to obtain the RMS position error at time step $k=51$ for each admissible course $\theta =-{170}^{\circ },{100}^{\circ },\dots ,{175}^{\circ }$. We obtained these RMS position errors by fixing the value of $\theta$ and conducting 500 Monte Carlo runs for each $\theta$ to compute the averaged RMS position errors.

The set of admissible control actions ${\mathbb{U}}_k$ is determined as follows. If the current position of the anti-submarine ship is ${\mathbf{u}}_k={\left[{\chi }_k{\psi }_k\right]}^T$, its future admissible locations are:

$${\mathbb{U}}_k=\left\{({\chi }_k+V_{ship}·cos(l{▵}_{\theta }+{\theta }_0),{\psi }_k+V_{ship}·sin(l{▵}_{\theta }+{\theta }_0));l=1,\dots ,N_{\theta }\right\}$$

(46)

where ${▵}_{\theta }=2\pi /N_{\theta }$ is a selected course step size, ${\theta }_0=-{50}^{\circ }$ is the initial course of the anti-submarine ship, $N_{\theta }=24$, and $V_{ship}=$ 4 knots. The anti-submarine ship can move in its current course (case $l=24$) or move in other courses. 24 control actions are considered at the end of the first motion leg.

We first get the RMS position errors for all the courses at the time step at which the control actions are executed. Then we test the online sensor control method by comparing the number of times out of 500 Monte Carlo runs that it has chosen for each particular course.

6.2.2. Experimental Results

The RMS position errors for different courses are plotted in Figure 13. It indicates that the second leg course decisions of ${70}^{\circ }$ and ${85}^{\circ }$ are preferred in this experiment. After finding out the best control actions for the anti-submarine ship’s second leg heading, we let the anti-submarine ship make its own decision by using the proposed online sensor control method. Figure 14 shows the number of times out of 500 Monte Carlo runs the online sensor control method has chosen for each particular course. From Figure 13 and Figure 14, we can know that the proposed online sensor control method can give a suitable control decision by using the digital twin-based framework of sensor control in ASW.

The results successfully verify the correctness and effectiveness of the proposed digital twin and RFS-based framework. In the following section, we will study the influence of the number M of prediction measurements on the results.

6.2.3. Sensitivity Analysis

We analyze the influence of the number M of prediction measurements on the proposed online sensor control method.

Table 2. Number of times (out of 500 Monte Carlo runs) a particular second leg course $\theta$ has been chosen for different values of parameter M.

$\mathit{\theta }\setminus \mathit{M}$	1	5	10	20	50	80	100	200	500	800	900	1000	1500	2000
$-{170}^{\circ }$	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$-{155}^{\circ }$	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$-{140}^{\circ }$	1	0	0	0	0	0	0	0	0	0	0	0	0	0
$-{125}^{\circ }$	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$-{110}^{\circ }$	8	0	0	0	0	0	0	0	0	0	0	0	0	0
$-{95}^{\circ }$	25	11	11	3	1	0	0	0	0	0	0	0	0	0
$-{80}^{\circ }$	45	23	37	28	22	15	17	11	9	5	7	6	8	2
$-{65}^{\circ }$	45	42	27	25	32	30	32	33	30	29	30	30	31	32
$-{50}^{\circ }$	22	30	29	22	20	16	20	19	19	26	18	15	20	17
$-{35}^{\circ }$	19	17	14	15	8	12	10	7	12	12	12	10	8	12
$-{20}^{\circ }$	2	7	10	7	7	8	6	8	6	4	5	8	6	5
$-5^{\circ }$	7	9	5	9	6	6	6	8	6	6	5	4	5	4
${10}^{\circ }$	4	3	6	4	3	5	3	2	3	6	1	7	4	7
${25}^{\circ }$	9	12	9	9	9	8	6	4	5	4	4	3	6	5
${40}^{\circ }$	18	10	12	19	14	13	13	9	11	9	16	13	12	11
${55}^{\circ }$	26	30	28	24	39	32	32	35	30	34	39	28	26	22
${70}^{\circ }$	65	90	82	96	108	120	110	108	104	102	92	106	94	91
${85}^{\circ }$	111	125	142	162	175	192	207	225	242	256	263	266	276	290
${100}^{\circ }$	64	77	81	76	56	43	38	31	23	7	8	4	4	2
${115}^{\circ }$	22	13	7	1	0	0	0	0	0	0	0	0	0	0
${130}^{\circ }$	5	1	0	0	0	0	0	0	0	0	0	0	0	0
${145}^{\circ }$	1	0	0	0	0	0	0	0	0	0	0	0	0	0
${160}^{\circ }$	0	0	0	0	0	0	0	0	0	0	0	0	0	0
${175}^{\circ }$	1	0	0	0	0	0	0	0	0	0	0	0	0	0

shows the results for different values of M. From Table 2, we know that when M increases, the number of times for the good control actions (including courses ${70}^{\circ }$ and ${85}^{\circ }$) increases, and the number of times for the bad control actions decreases. However, the influence extent is very limited, since M increases quickly and the performance of the method is slightly improved.

6.3. Online Sensor Control Experiment with Multiple Submarines

In this experiment, we verify the effectiveness of the proposed online sensor control method by a scenario where the anti-submarine ship tracks multiple submarines. Here we control the range-only sensor of the anti-submarine ship by using the proposed digital twin and RFS-based method.

6.3.1. Experimental Setup

The state of single submarine at time step k is represented by ${\mathbf{x}}_{\mathbf{k}}={\left[{\mathbf{p}}_{\mathbf{k}}^{\mathbf{T}}{\mathbf{v}}_{\mathbf{k}}^{\mathbf{T}}\right]}^T$, here ${\mathbf{p}}_{\mathbf{k}}={\left[x_k{y}_k\right]}^T$ is the position, ${\mathbf{v}}_{\mathbf{k}}={\left[v_{k,x}{v}_{k,y}\right]}^T$ is the velocity. The position of the anti-submarine ship is represented by ${\mathbf{u}}_k^o={\left[x_k^o{y}_k^o\right]}^T$. The probability of detecting the enemy submarine at position ${\mathbf{p}}_{\mathbf{k}}={\left[x_k{y}_k\right]}^T$ is computed as follows:

$$p_D\left({\mathbf{p}}_k\right)=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}\parallel {\mathbf{p}}_k-{\mathbf{u}}_k^o\parallel \le R_0\hfill \\ max\{0,1-(\parallel {\mathbf{p}}_k-{\mathbf{u}}_k^o\parallel -R_0\left)\hslash \right\},\hfill & \mathrm{if}\parallel {\mathbf{p}}_k-{\mathbf{u}}_k^o\parallel >R_0\hfill \end{array}\right.$$

(47)

here $\parallel {\mathbf{p}}_k-{\mathbf{u}}_k^o\parallel =\sqrt{{(x_k-x_k^o)}^2+{(y_k-y_k^o)}^2}$ is the distance between the sensor and the submarine at ${\mathbf{p}}_k$. In this experiment, the operational area is a square of sides $s=1200$ m, $R_0=300$ and $\hslash =0.0002$ m${}^{-1}$. The measurement is generated by $z=h\left({\mathbf{x}}_{\mathbf{k}}\right)=\parallel {\mathbf{p}}_k-{\mathbf{u}}_k^o\parallel +\omega$, where $\omega$ is zero-mean white Gaussian measurement noise, with deviation ${\sigma }_{\omega }={\sigma }_0+\beta {\parallel {\mathbf{p}}_k-{\mathbf{u}}_k^o\parallel }^2$. In this experiment, ${\sigma }_0=1$ m and $\beta =5\times {10}^{-5}$ m${}^{-1}$. The clutters are modeled as the Poisson RFS. The intensity of clutters is modeled by the uniform density$\kappa \left(z\right)=\lambda ·c\left(z\right)$ with mean $\lambda =5$.

In this experiment, there are 5 moving enemy submarines in the operational area. As it is shown in (46), we control the sensor by finding the optional course. To prove the validity of the proposed method, we compare the proposed control method with the other method that the control vector is randomly selected from the set ${\mathbb{U}}_k$. Each method runs 10 times, and the averaged OSPA errors (order parameter $p=2$ and cutoff $c=100$ m) are compared.

6.3.2. Experimental Results

In this experiment, the anti-submarine ship controls the sensor and runs the sensor control method every $T_c=5$ time steps. We compare the performance by using the OSPA error at every time step. The mean OSPA errors of two methods are given in Figure 15. We can see that the proposed digital twin and RFS-based online sensor control method can effectively reduce the OSPA error. We can also find that as time evolves, the OSPA error of the proposed method also gradually reduces in this application.

Figure 16 gives the estimated submarines’ numbers of 10 Monte Carlo simulations generated by using the proposed data-assimilation algorithm and the different online sensor control methods. The black line represents the truth of submarines’ number, and the data points of different shapes represent the estimated submarines’ number at each time step of 10 Monte Carlo simulations. Figure 16a is the result of the proposed sensor control method, and Figure 16b is the result of the random control method. We can find that the proposed sensor control method performs much better than the random control method on estimating submarines’ number.

Figure 17 gives the estimated submarines’ states of 10 Monte Carlo simulations generated by using the proposed data-assimilation algorithm and the different online sensor control methods. The black lines represent the truth of submarines’ states, and the data points of different shapes represent the truth of the submarines’ states. Figure 17a is the result of the proposed sensor control method, and Figure 17b is the result of the random control method. We can find that the proposed sensor control method also performs much better than the random control method for estimating submarines’ states.

The paths of the anti-submarine ship for two methods are shown in Figure 18. Different colors represent different Monte Carlo simulations, and the red points represent the truth of the submarines’ trajectories. We can see that the proposed sensor control method can successfully guide the sensor to move close to the enemy submarines to get more accurate and more reliable measurements.

6.3.3. Sensitivity Analysis

In this section, the performance of the prosed online sensor control method is further analyzed by sensitivity analysis. We analyze some parameters’ influence on the performance of the proposed sensor control method. These parameters are the time interval $T_c$ of two sensor control actions, factor $\alpha$, the number M of prediction measurements, and particle number N for each submarine.

From Figure 19, we can see that as $T_c$ increases, the mean OSPA error also increases. This means that decreasing $T_c$ can improve the performance of the prosed online sensor control method. Figure 20 gives the OSPA error for various values of parameter $\alpha$. We can find that $\alpha$ has quite little influence on the performance of the proposed sensor control method in this application. The influence of the number M of prediction measurements on the online sensor control method is shown in Figure 21. We can see that the influence of M on the performance of the proposed sensor control method is little. The reason is that the performance of the proposed method does not only depend on the reward function, but also depends on the data-assimilation algorithm. As it is shown in Figure 22, the particle number N can affect the performance of the proposed method. The increase of N can ensure more reasonable choice of the control action.

7. Conclusions

In this paper, we studied the digital twin-based framework of sensor control in ASW. We firstly combine the simulated ASW with the real ASW by employing the digital twin theory. Then we proposed an RFS-based data-assimilation algorithm to dynamically incorporate online measurements generated from the real ASW. At last, we also fostered the ability of the simulation system to control the sensor in ASW by deriving and implementing the data-assimilation-based reward function. The proposed data-assimilation algorithm has the potential to overcome the limitations of the conventional vector-based algorithms. It can jointly estimate the number of enemy submarines and the state of each enemy submarine. We tested the proposed data-assimilation algorithm by using the identical-twin experiment. The results prove that the proposed algorithm can assimilate the input measurements and improve the accuracy of simulation results. We tested the proposed online sensor control method with two group of experiments, including the scenario with single submarine and the scenario with multiple submarines. The results showed that the proposed online sensor control method can effectively control the sensor. This paper can be regarded as an application of digital twin in ASW and the methods can also be applied to other applications.