The Situation Assessment of UAVs Based on an Improved Whale Optimization Bayesian Network Parameter-Learning Algorithm

Abstract

To realize unmanned aerial vehicle (UAV) situation assessment, a Bayesian network (BN) for situation assessment is established. Aimed at the problem that the parameters of the BN are difficult to obtain, an improved whale optimization algorithm based on prior parameter intervals (IWOA-PPI) for parameter learning is proposed. Firstly, according to the dependencies between the situation and its related factors, the structure of the BN is established. Secondly, in order to fully mine the prior knowledge of parameters, the parameter constraints are transformed into parameter prior intervals using Monte Carlo sampling and interval transformation formulas. Thirdly, a variable encircling factor and a nonlinear convergence factor are proposed. The former and the latter enhance the local and global search capabilities of the whale optimization algorithm (WOA), respectively. Finally, a simulated annealing strategy incorporating Levy flight is introduced to enable the WOA to jump out of the local optimum. In the experiment for the standard BNs, five parameter-learning algorithms are applied, and the results prove that the IWOA-PPI is not only effective but also the most accurate. In the experiment for the situation BN, the situations of the assumed mission scenario are evaluated, and the results show that the situation assessment method proposed in this article is correct and feasible.

1. Introduction

The modern working environments for UAVs are extremely complex and filled with a wide variety and a large number of entities, and these bring difficulties and challenges for UAV situation assessment. Situation assessment is the process of perceiving the attributes, states, and behaviors of entities in an environment, understanding an environment’s information, inferring entities’ intentions, and finally predicting entities’ future short-term actions [1]. Intention inference is the core of situation assessment, and the mission is the expression of the intention. Therefore, the situation is embodied as the entities’ mission in this article. The general methods of situation assessment include the analytic hierarchy process (AHP) [2,3], the technique for order preference by similarity to ideal solution (TOPSIS) [4,5], neural networks [6,7], fuzzy logic [8,9], Bayesian networks [10,11], etc.

All the above methods, except for Bayesian networks, achieve situation assessment to a certain degree, but each of them has some shortcomings. The AHP relies on the human subjective experience too much, and inaccurate experiences lead to a large deviation in situation assessment results. The TOPSIS is unable to solve the problem of single-target situation assessment. In the case of multiple targets, the distance from the optimal solution to the positive ideal solution may approximate the distance to the negative ideal solution, which makes it possible to obtain the opposite assessed results. Neural networks have black box characteristics that make it hard to accurately determine the effect on the output when certain inputs change. It is difficult to determine the logic operations and inferring methods of fuzzy logic. In addition, both neural networks and fuzzy logic have only a single output for known inputs and cannot provide other feasible solutions.

A Bayesian network (BN) is a probabilistic graphical model that combines the abilities to express and infer uncertain knowledge and is capable of handling multivariate information. BNs consist of a structure and parameters [12]. The former abstracts the dependencies between the situations and its associated elements into a visual network, and the latter transforms the degree of association between the elements into a conditional probability table (CPT). The construction process of BNs is consistent with the human cognitive habits of situations, which is more conducive to commanders’ inference and application. Compared to the AHP, thanks to structure and parameter learning, BNs alleviate the influence of human subjective experiences. In comparison with the TOPSIS, BNs solve the problem of single-target and multi-target situation assessment. Speaking of neural networks and fuzzy logic, BNs not only have a clear structure with white box characteristics, but also follow mature Bayesian criteria in probability calculation and inference. In addition, BNs output the probabilities of different states of multiple nodes and provide decision-makers with multiple choices.

Despite the many advantages of BNs, the network parameters are sometimes difficult or even impossible to obtain. Combined with expert knowledge and sample data, the parameters are acquired via parameter learning. The main parameter-learning methods include constrained optimization methods and Bayesian estimation methods.

The constrained optimization methods represent expert knowledge as constraints and then treat parameter learning as an optimization problem. Niculescu [13], Campos [14], and Hou [15] transform parameter learning into an optimization problem with constraints. They use the logarithmic or entropy function as the objective function and optimize the objective function convexly in the feasible domain, limited by constraints. Altendorf [16] and Liao [17] transformed parameter learning into an optimization problem without constraints. They constructed a penalty function with constraints and then summed the logarithmic likelihood function and the penalty function to obtain the augmented objective function to be optimized.

The Bayesian estimation methods transform expert knowledge into the ranges of the values of the parameters, calculate the hyper-parameters of the prior distribution of the parameters within these ranges, and finally combine the hyper-parameters and the sample data to compute the Bayesian maximum a posteriori (MAP) to obtain the BN parameters. Ren [18] and Di [19] assume that the parameters obey a uniform distribution in the feasible domain, and Chai [20] proposes to represent the approximate equality constraints based on the normal distribution. They use Beta distribution to approximate the uniform and normal distributions, respectively. Since the Beta distribution is the marginal distribution of the Dirichlet distribution, the hyper-parameters of the Dirichlet distribution are derived from the Beta distribution parameters. In the end, the BN parameters are obtained using the MAP formula. Gao [21] proposes a constrained Bayesian estimation (CBE) algorithm that enhances learning accuracy by introducing expert criteria. Di [22] proposes a constrained adjusted MAP (CaMAP) algorithm by choosing a reasonable equivalent sample size. The qualitative maximum a posteriori estimation (QMAP) algorithm proposed by Chang [23] performs Monte Carlo (MC) sampling on the feasible domain of the parameters determined based on constraints. The pseudo prior counts obtained are functionally equivalent to the hyper-parameters of the Dirichlet distribution, and then the parameters are computed using the MAP formula. Guo [24] made an improvement to the QMAP algorithm and proposes a further constrained QMAP algorithm.

The above two methods have the following defects: firstly, the learning results of the convex optimization algorithm in the studies [13,14,15] are often located at the boundary of the feasible domain of the parameters [25], which leads to the degradation of inequality constraints into equality constraints. This indicates a failure to fully utilize expert knowledge. In addition, directly optimizing within the limited range of ordinary constraints often fails to fully utilize prior knowledge in the constraints. Secondly, the penalty function method in studies [16,17] needs to design specific penalty functions for different constraints, and the penalty factor is manually determined, which is sometimes inaccurate. Thirdly, for some sample data, the learning results of QMAP [23] may violate some parameter constraints. This means a failure to make full use of the sample data. In addition, study [19] only applies to monotonic constraints, and study [20] only applies to approximate equality constraints. The prior distribution of [21] is set to be the BDeu priors rather than the transferred priors, which are more meaningful. When no parameter constraints are available, the CaMAP [22] is inferior to the MAP. The FC-QMAP algorithm [24] only outperforms QMAP when dealing with small datasets.

In order to overcome the shortcomings of the above algorithms, in this article, an improved whale optimization algorithm based on parameter prior intervals (IWOA-PPI) for parameter learning is proposed, and the parameters learned by the IWOA-PPI are substituted into the situation assessment BN to evaluate the situation. In this algorithm, expert knowledge is transformed into parameter qualitative constraints; then, parameter prior intervals (PPIs) are calculated depending on the qualitative constraints, and finally, the improved WOA is used to optimize within the PPIs to obtain the optimal parameters. Due to the integration of the advantages of the constrained optimization methods and Bayesian estimation methods, the new algorithm proposed in this article can not only mine the information in the expert knowledge as much as possible, but can also fully utilize sample data. The main contributions of this article are that (1) the concept of PPIs is proposed. The PPIs not only contain the prior knowledge of the parameters but also narrow the feasible domain of the parameters, thus improving the search accuracy of subsequent IWOA. (2) A new variable encircling factor for the WOA is proposed. The variable encircling factor can properly shrink the local search area with the operation of the algorithm, and the local search ability of the WOA is enhanced.

The rest of the article is organized as follows: Section 2 introduces the preliminary knowledge, including BN parameter learning, parameter constraints, and the WOA. Section 3 establishes the structure of the BN for situation assessment. Section 4 proposes the improved whale optimization algorithm based on parameter prior intervals. Section 5 designs the simulation experiments. The parameters of four BNs are learned by five algorithms in the experiment for standard BNs, and the results prove that the IWOA-PPI proposed in this article is of the highest accuracy. The situation of the assumed mission scenario is evaluated in the experiment for the situation assessment BN, and the results verify that the BN established in this article can correctly infer the target intention. Finally, conclusions are drawn in Section 6.

2. Preliminaries

2.1. BN Parameter Learning

A BN is a direct acyclic graph with probabilities, which is generally denoted as $B=\left(G,\Theta \right)$. $G=\left(V,E\right)$ is the graph structure, the set of nodes ($V$) denotes the set of random variables, and the set of directed arcs ($E$) denotes the dependencies between the random variables. $\Theta$ is the set of network parameters denoting the conditional probability distributions between states of nodes. $\Theta$ is also denoted as $P\left(X_i|pa\left(X_i\right)\right)$, where $pa\left(X_i\right)$ denotes the set of parents of node $X_i$. When the parents are given, node $X_i$ is conditionally independent of its non-descendant nodes. According to the Markov condition, the joint probability distribution of a BN can be represented in Equation (1):

(1)$P\left(X_1,X_2,\cdots ,X_n\right)={\displaystyle \mathbf{\prod }_{i=1}^{n}P\left(X_i|pa\left(X_i\right)\right)}$

When the structure $G$ is known, the parameter learning of a BN refers to the estimation of the true parameters of the BN from a given sample dataset, $D$, according to a certain criterion. The BN parameter ${\theta }_{ijk}$ is denoted in Equation (2):

(2)$\begin{array}{l}{\theta }_{ijk}=P\left(X_i=k|pa\left(X_i\right)=j\right),\\ 1\le i\le n,1\le j\le q_i,1\le k\le r_i\end{array}$

where $i$ denotes the serial number of $X_i$, $j$ denotes the state of $pa\left(X_i\right)$, and $k$ denotes the state of $X_i$.

There are two basic methods for BN parameter learning: maximum likelihood estimation (MLE) and Bayesian estimation. MLE is usually used for a sample set with a sufficient sample size and no missing data. The logarithmic likelihood function for BNs is defined in Equation (3):

(3)$L\left(\theta |D\right)=\log P\left(D|\theta \right)={\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^{q_i}{\displaystyle {\sum }_{k=1}^{r_i}{N}_{ijk}\log {\theta }_{ijk}}}}$

where $N_{ijk}$ denotes the number of samples that satisfy $X_i=k$ and $pa\left(X_i\right)=j$. The ${\theta }_{ijk}$ that maximizes $L\left(\theta |D\right)$ is called MLE, and is shown in Equation (4):

(4)${\theta }_{ijk}=\frac{N_{ijk}}{{\displaystyle \sum _{k=1}^{r_i}{N}_{ijk}}}=\frac{N_{ijk}}{N_{ij}}$

where $N_{ij}$ denotes the number of samples that satisfy $pa\left(X_i\right)=j$.

MAP is usually used for a sample set with an insufficient sample size and no missing data. Assuming that the prior distribution, $P\left(\theta \right)$, of the BN parameter satisfies the Dirichlet distribution $D\left({\alpha }_{ij1},{\alpha }_{ij1},\cdots {\alpha }_{ij{r}_i}\right)$, and since the Dirichlet distribution is a conjugate family of polynomial distributions, the posterior distribution $P\left(\theta |D\right)$ also satisfies the Dirichlet distribution, which is expressed in Equation (5):

(5)$P\left(\theta |D\right)\propto {\displaystyle \mathbf{\prod }_{i=1}^n{\displaystyle \mathbf{\prod }_{j=1}^{q_i}{\displaystyle \mathbf{\prod }_{k=1}^{r_i}{\theta }^{N_{ijk}+{\alpha }_{ijk}-1}}}}$

At this point, the logarithmic likelihood function for the BN is defined in Equation (6):

(6)$\log P\left(\theta |D\right)=\log P\left(\theta \right)P\left(D|\theta \right)-c$

where $c$ is a constant. The ${\theta }_{ijk}$ that maximizes $P\left(\theta |D\right)$ is called MAP, and is shown in Equation (7):

(7)${\theta }_{ijk}=\frac{N_{ijk}+{\alpha }_{ijk}}{{\displaystyle \sum _{k=1}^{r_i}{N}_{ijk}+{\alpha }_{ijk}}}$

where ${\alpha }_{ijk}$ is the hyperparameter of the Dirichlet distribution.

2.2. Parameter Constraints

Parameter constraints refer to the constraint relationships between the conditional probabilities of a BN given by domain experts. There are generally five types of parameter constraints.

• Axiomatic constraint.

This constraint defines a special sum-of-parameters relationship and is given in Equation (8). The sum of the probabilities of the various states of the child nodes with constant parent combinations is 1.

(8)${\displaystyle \sum _{k=1}^{r_i}{\theta }_{ijk}}=1,0\le {\theta }_{ijk}\le 1,\forall i,j,k$

• 2.. Range constraint.

This constraint defines the upper and lower bounds of the parameter and is given in Equation (9).

(9)$0\le {\alpha }_{ijk}\le {\theta }_{ijk}\le {\beta }_{ijk}\le 1$

where ${\alpha }_{ijk}$ and ${\beta }_{ijk}$ are the bounds of the parameter.

• 3.. Approximate equality constraint.

This constraint defines the approximate equality relationship between the parameters and is given in Equation (10).

(10)${\theta }_{ijk}\approx {\theta }_{i^{\prime }{j}^{\prime }{k}^{\prime }},\forall ijk\ne i^{\prime }{j}^{\prime }{k}^{\prime }$

Equation (10) is not easy to apply, so the form of Equation (11) is usually adopted:

(11)$\left|{\theta }_{ijk}-{\theta }_{i^{\prime }{j}^{\prime }{k}^{\prime }}\right|<\epsilon ,\forall i\ne i^{\prime },j\ne j^{\prime },k\ne k^{\prime }$

where $\epsilon$ is a very small positive rational number.

• 4.. Inequality constraint.

This constraint defines the inequality relationship between the parameters and includes three types that are given in Equations (12)–(14):

• Intra-distribution constraint:

(12)${\theta }_{ijk}\le {\theta }_{ij{k}^{\prime }},\forall k\ne k^{\prime }$

• Cross-distribution constraint:

(13) ${\theta }_{ijk}\le {\theta }_{i{j}^{\prime }k},\forall j\ne j^{\prime }$

• Inter-distribution constraint:

(14) ${\theta }_{ijk}\le {\theta }_{i^{\prime }{j}^{\prime }{k}^{\prime }},\forall i\ne i^{\prime },j\ne j^{\prime },k\ne k^{\prime }$

• 5.. Synergy constraint.

This constraint defines the inequality relationship between the sum or product of multiple parameters and includes two types that are given in Equations (15) and (16):

• Additive synergy constraint:

(15)${\theta }_{i{j}_{1}k}+{\theta }_{i{j}_{2}k}\le {\theta }_{i{j}_{3}k}+{\theta }_{i{j}_{4}k},\forall j_1\ne j_2\ne j_3\ne j_4$

• Product synergy constraint:

(16) ${\theta }_{i{j}_{1}k}\cdot {\theta }_{i{j}_{2}k}\le {\theta }_{i{j}_{3}k}\cdot {\theta }_{i{j}_{4}k},\forall j_1\ne j_2\ne j_3\ne j_4$

2.3. Whale Optimization Algorithm

The whale optimization algorithm (WOA) is a meta-heuristic intelligent optimization algorithm proposed by Mirjalili [26], inspired by the natural phenomenon of the predatory behavior of humpback whales. This algorithm has the advantages of a high optimization accuracy and few control factors, and has been applied in many fields such as photovoltaic power generation [27], resource scheduling [28], and path planning [29]. In addition, some scholars have made improvements to the whale algorithm [30,31,32].

The WOA is divided into three parts, namely encircling prey, bubble-net attacking (exploitation phase), and the search for prey (exploration phase).

• 1.. Encircling prey:

The WOA assumes that the position of the individual that reaches the optimal solution in the current population is the prey position, and then other whales will move towards that position for the purpose of encirclement. The mathematical model for encircling prey is expressed in Equations (17)–(20):

(17)$\overrightarrow{D}=\left|\overrightarrow{C}\cdot \overrightarrow{X^{*}}\left(t\right)-\overrightarrow{X}\left(t\right)\right|$

(18)$\overrightarrow{X}\left(t+1\right)=\overrightarrow{X^{*}}\left(t\right)-\overrightarrow{A}\cdot \overrightarrow{D}$

(19)$\overrightarrow{A}=2\overrightarrow{a}{\overrightarrow{r}}_1-\overrightarrow{a}$

(20)$\overrightarrow{C}=2{\overrightarrow{r}}_2$

where $t$ is the iteration number, $\overrightarrow{X^{*}}\left(t\right)$ is the optimal position, and $\overrightarrow{X}\left(t\right)$ is the individual position. $\overrightarrow{A}$ and $\overrightarrow{C}$ are the coefficient vectors determined by $\overrightarrow{a}$, ${\overrightarrow{r}}_1$, and ${\overrightarrow{r}}_2$. $\overrightarrow{a}$ decreases linearly from 2 to 0 with each iteration. ${\overrightarrow{r}}_1$ and ${\overrightarrow{r}}_2$ are random vectors in [0, 1].

• 2.. Bubble-net attacking (exploitation phase)

Bubble-net attacking consists of two parts: shrinking the encircling mechanism and the spiral updating position. When $\left|\overrightarrow{A}\right|\le 1$, the former is realized in Equations (17) and (18). The mathematical model of the latter is expressed in Equations (21) and (22):

(21)$\overrightarrow{D^{\prime }}=\left|\overrightarrow{X^{*}}\left(t\right)-\overrightarrow{X}\left(t\right)\right|$

(22)$\overrightarrow{X}\left(t+1\right)=\overrightarrow{X^{*}}\left(t\right)+\overrightarrow{D^{\prime }}\cdot e^{bl}\cdot \cos \left(2\pi l\right)$

where $b$ defines the shape of the logarithmic helix, and $l$ is a random number in [−1, 1].

When the whales implement bubble-net attacking, the shrinking of the encirclement and the spiral update occur simultaneously. To model simultaneity, it is assumed that there is a 50% probability of selecting one of the behaviors each time, as shown in Equation (23):

(23)$\overrightarrow{X}\left(t+1\right)=\left\{\begin{array}{cc}\overrightarrow{X^{*}}\left(t\right)-\overrightarrow{A}\cdot \overrightarrow{D}& if\begin{array}{c}\end{array}p<0.5\\ \overrightarrow{X^{*}}\left(t\right)+\overrightarrow{D^{\prime }}\cdot e^{bl}\cdot \cos \left(2\pi l\right)& if\begin{array}{c}\end{array}p\ge 0.5\end{array}\right.$

where $p$ is a random number in [0, 1].

• 3.. Search for prey (exploration phase)

When $\left|\overrightarrow{A}\right|>1$, whales search randomly according to each other’s position. The mathematical model is expressed in Equations (24) and (25):

(24)$\overrightarrow{D}=\left|\overrightarrow{C}\cdot \overrightarrow{X_{rand}}\left(t\right)-\overrightarrow{X}\left(t\right)\right|$

(25)$\overrightarrow{X}\left(t+1\right)=\overrightarrow{X_{rand}}\left(t\right)-\overrightarrow{A}\cdot \overrightarrow{D}$

where $\overrightarrow{X_{rand}}\left(t\right)$ is the position vector of a random individual in the current population.

3. Structure Establishment of Situation Assessment BN

In this section, the influencing factors of the target intention are enumerated, and the structure of the situation assessment BN is established based on their dependencies.

When a target performs a mission, its intention can be reflected by its attributes and state to a certain extent. Firstly, by restricting the range of missions that the target can perform, the target type is an important attribute to infer the intention. Secondly, relative motion is the movement trend of the target relative to our UAVs, which can indicate that the target is approaching or leaving our UAVs. Finally, the relative velocity and relative height between the target and our UAVs are able to serve as an auxiliary basis for inferring the target intention. Therefore, the BN for intention inference is shown in Figure 1.

Similar to target intention inference, target type recognition requires obtaining target information, which generally includes the velocity, height, radar cross-section (RCS), and radar frequency band variability (RFV) of the target. Therefore, the BN for type recognition is shown in Figure 2.

The relative motion trend of the target is determined based on the distance and bearing relative of the target to our UAVs. Therefore, the BN for relative motion determination is shown in Figure 3.

By connecting the three sub-networks in Figure 1, Figure 2 and Figure 3 into a multi-layer BN, the structure of the situation assessment BN is established, as shown in Figure 4. Among them, the light-colored nodes such as velocity, height, etc., are the base nodes, and their states can be obtained via observation. The three dark-colored nodes, namely target intention, target type, and relative motion, are the hidden nodes that use base nodes as their child nodes. In conjunction with the CPTs and the states of the child nodes, the probabilities of each state of the hidden nodes are calculated and acquired. The relevant information on each node in the network is shown in Table 1.

4. IWOA-PPI for Parameter Learning

The establishment of the structure for a situation assessment BN is introduced in Section 3, and the IWOA-PPI for parameter learning is discussed in detail in this section. The parameter learning problem is treated as an unconstrained optimization problem in the PPIs, and the objective function can be acquired using Equation (3). The idea of the IWOA-PPI is to transform parameter constraints into PPIs, and then use the improved WOA to search for the optimal parameters. In this way, the algorithm has the advantages of Bayesian estimation, which fully utilizes prior information, and constrained optimization, which fully utilizes sample data. Therefore, the key to the algorithm is the PPIs and the improvements to the WOA. The former has the ability to mine the parameter prior information embedded in the expert knowledge, thereby improving the accuracy of the parameter learning. For the latter, since the original WOA has the shortcomings of a slow convergence speed, an inability to jump out of the local optima, and difficulty in reaching the global optima, it is necessary to make appropriate improvements to accelerate the search speed and enhance the global optimization ability. In this section, the establishment of the PPIs is described in Section 4.1, and the improvements to the WOA are explained in Section 4.2, Section 4.3 and Section 4.4.

4.1. Parameter Prior Interval

Expert knowledge is defined as the qualitative parameter relationships that are summarized by domain experts based on physical phenomena or objective laws. The knowledge contains prior information of the parameters and is typically represented as parameter constraints. In order to make full use of prior knowledge, the concept of PPIs is proposed.

PPIs are essentially the intervals bounded by the upper and lower bounds of the parameters, formally equivalent to the range constraint of Equation (9), but fundamentally different from general parameter constraints. Firstly, the PPI is a concise upper and lower bound form, but the parameter constraints are embodied in a variety of more complex forms. Secondly, the PPI is transformed by extracting the prior information of the parameters in the parameter constraints, so the prior knowledge can be fully utilized. Finally, compared to parameter constraints, PPIs narrow the parameters’ feasible domain to a more precise region, which accelerates the subsequent IWOA speed and improves its accuracy.

PPIs are obtained in three steps:

• By performing MC sampling on the parameter space delimited by the parameter constraints, the parameters without samples are obtained using Equation (26):

  (26)${\theta }_{ijk}^{MC}=P\left(X_i=k,pa\left(X_i\right)=j|\Omega \right)=\frac{{\displaystyle \sum _{l=1}^S{P}_l\left(X_i=k,pa\left(X_i\right)=j|\Omega \right)}}{S}$

  where $\Omega$ denotes the parameters constraints, $S$ denotes the number of samples of MC sampling, $P_l\left(X_i=k,pa\left(X_i\right)=j|\Omega \right)$ denotes the value of ${\theta }_{ijk}^l$ for a single sample, and ${\theta }_{ijk}^{MC}$ denotes the mean value of all ${\theta }_{ijk}^l$ for $S$ samples.

  ${\theta }_{ijk}^{MC}$ is the parameter obtained via MC sampling without considering the samples. Its numerical value is close to the true parameter to some extent, and the degree of closeness is positively related to the number of constraints.

• Using the interval transform formulas, the uncorrected PPIs are obtained. ${\theta }_{ijk}^{MC}$ ($k=1,2,\dots ,r_i$) are denoted from small to large, as shown in Equation (27):

  (27)${\theta }_{ij{k}_1}^{MC}\le {\theta }_{ij{k}_2}^{MC}\le \dots \le {\theta }_{ij{k}_m}^{MC}\le \dots \le {\theta }_{ij{k}_{r_i}}^{MC}$

  where $k_m\in \left\{1,2,\dots ,r_i\right\}$.

  $d_{k_m}^{low}$ is defined as the lower bound interval of ${\theta }_{ij{k}_m}^{MC}$, $d_{k_m}^{up}$ is defined as the upper bound interval of ${\theta }_{ij{k}_m}^{MC}$, and $[{\theta }_{ij{k}_m}^{MC}-d_{k_m}^{low},{\theta }_{ij{k}_m}^{MC}+d_{k_m}^{up}]$ is the uncorrected PPI of ${\theta }_{ij{k}_m}^{MC}$. The determination of $d_{k_m}^{low}$ and $d_{k_m}^{up}$ needs to be categorized into two cases:

• When ${\theta }_{ij{k}_1}^{MC}<{\theta }_{ij{k}_2}^{MC}<\dots <{\theta }_{ij{k}_m}^{MC}<\dots <{\theta }_{ij{k}_{r_i}}^{MC}$, the interval transform formulas of ${\theta }_{ij{k}_m}^{MC}$ are implemented as Equations (28) and (29):

  (28)$d_{k_m}^{low}=\frac{{\theta }_{ij{k}_m}^{MC}-{\theta }_{ij{k}_{m-1}}^{MC}}{{\theta }_{ij{k}_m}^{MC}/{\theta }_{ij{k}_{m-1}}^{MC}}$

  (29)$d_{k_m}^{up}=\frac{{\theta }_{ij{k}_{m+1}}^{MC}-{\theta }_{ij{k}_m}^{MC}}{{\theta }_{ij{k}_{m+1}}^{MC}/{\theta }_{ij{k}_m}^{MC}}$

  For ${\theta }_{ij{k}_1}^{MC}$, $d_{k_1}^{up}$ is calculated first, and then $d_{k_1}^{low}=d_{k_1}^{up}$. For ${\theta }_{ij{k}_{r_i}}^{MC}$, $d_{k_{r_i}}^{low}$ is calculated first, and then $d_{k_{r_i}}^{up}=d_{k_{r_i}}^{low}$.

• When ${\theta }_{ij{k}_p}^{MC}={\theta }_{ij{k}_{p+1}}^{MC}=\dots ={\theta }_{ij{k}_{p+q}}^{MC}$, the interval transform formulas of ${\theta }_{ij{k}_m}^{MC}$ are implemented as Equation (30):

  (30)$d_{k_m}^{low}=d_{k_m}^{up}=\frac{\omega \cdot \left({\theta }_{ij{k}_p}^{MC}+{\theta }_{ij{k}_{p+1}}^{MC}+\dots +{\theta }_{ij{k}_{p+q}}^{MC}\right)}{q}$

  where $\omega$ is the weight coefficient that takes a small value.

• 3.. Combining the parameter constraints, the PPIs are obtained. The uncorrected PPIs may violate some parameter constraints, so the intersection of the parameter constraints and the uncorrected PPIs is sought to obtain the PPIs. For example, in order to satisfy the axiomatic constraint, the lower bound interval, $d_{k_1}^{low}$, of ${\theta }_{ij{k}_1}^{MC}$ is $\max (0,{\theta }_{ij{k}_1}^{MC}-d_{k_1}^{low})$ and the upper bound interval, $d_{k_1}^{up}$, of ${\theta }_{ij{k}_1}^{MC}$ is $\min ({\theta }_{ij{k}_{r_i}}^{MC}+d_{k_i}^{up},1)$.

4.2. Variable Encircling Factor

Aimed at the encircling prey behavior, a variable encircling factor is proposed to enhance the local search ability of the WOA.

The value range of $\overrightarrow{C}\cdot \overrightarrow{X^{*}}(t)$ in Equation (17) is $[\overrightarrow{0},2\overrightarrow{X^{*}}(t)]$. The physical meaning of $[\overrightarrow{0},2\overrightarrow{X^{*}}(t)]$ refers to a high-dimensional neighborhood. This neighborhood is centered on the current optimal position, $\overrightarrow{X^{*}}(t)$, and the range of each dimension is $[0,2X_i^{*}(t)]$. This neighborhood is embodied as a rectangle in two dimensions, a rectangular body in three dimensions, and a hyper-cuboid in high dimensions. The encircling prey behavior of Equations (17) and (18) is a random selection of a point, $\overrightarrow{C}\cdot \overrightarrow{X^{*}}(t)$, within the neighborhood, $[\overrightarrow{0},2\overrightarrow{X^{*}}(t)]$, and then all whale individuals shrink and move from their current position, $\overrightarrow{X}(t)$, to $\overrightarrow{C}\cdot \overrightarrow{X^{*}}(t)$.

$\overrightarrow{C}$ is defined as the encircling factor, and its value range determines the size of $[0,2X_i^{*}\left(t\right)]$. In the original WOA, all dimensional variables of ${\overrightarrow{r}}_2$ in Equation (20) obey the uniform distribution, $U\left(0,1\right)$; i.e., all dimensional variables of $\overrightarrow{C}$ obey the uniform distribution, $U\left(0,2\right)$. As the algorithm runs, $\overrightarrow{X^{*}}(t)$ gradually approaches the optimal position. Since the neighborhood $[\overrightarrow{0},2\overrightarrow{X^{*}}(t)]$ is too large, this leads to a large number of invalid local searches. Therefore, a variable encircling factor, $\overrightarrow{C^V}$, obeying the normal distribution, is proposed as Equation (31):

(31)$p\left(C_i^V\right)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\left(C_i^V-1\right)}^2}{2{\sigma }^2}}$

where $C_i^V$ is a one-dimensional variable of $\overrightarrow{C^V}$.

The variability of the variable encircling factor, $\overrightarrow{C^V}$, is reflected in $\sigma$: firstly, the value of $\sigma$ is very large at the beginning of the algorithm. At this moment, the distribution of $C_i^V$ degenerates from the normal distribution, $N\left(1,{\sigma }^2\right)$, to the uniform distribution, $U\left(0,2\right)$, and it means a degeneration towards the original WOA. Secondly, after certain iterations, make $3\sigma$ equal to 1, and then $C_i^V$ obeys the normal distribution, $N\left(1,1/3\right)$. Then, the $3\sigma$ area of $C_i^V$ is $[0,2]$, which significantly increases the probability of $\overrightarrow{C}\cdot \overrightarrow{X^{*}}(t)$ approaching $\overrightarrow{X^{*}}(t)$. Finally, $3\sigma$ decreases gradually to 0.1 at the later stage of the algorithm, and the $3\sigma$ area of $C_i^V$ becomes $[0.9,1.1]$ in the end. The process of changing $\sigma$ shrinks the local search area to accelerate the local search speed.

4.3. Nonlinear Convergence Factor

The WOA performs a global search in the exploration phase and a local search in the exploitation phase, and switches between the two phases based on control coefficient $\overrightarrow{A}$, determined as $a$. As $a$ decreases linearly, both global and local searches are performed in the early stages ($a>1$), and only a local search is performed in the later stages ($a\le 1$). However, when dealing with certain optimization problems, the linear convergence factor $a$ makes the algorithm enter the later phase so early that a sufficient global search has not been performed in the early stage. Therefore, the global optimal solution cannot be found. To overcome this problem, a nonlinear convergence factor, $a_{nl}$, is proposed as Equation (32):

(32)$a_{nl}=2\cdot \cos \left(\frac{\pi }{2}\cdot \frac{t}{t_{\max }}\right)$

where $t$ is the iteration number and $t_{\max }$ is the maximum iteration number.

The curves of $a$ and $a_{nl}$ with iterations are shown in Figure 5. As shown in Figure 5, $a_{nl}$ has two characteristics: firstly, from the perspective of the number of iterations, the number of iterations where $a>1$ is 250, while the number of iterations where $a_{nl}>1$ is extended to about 330. This indicates that $a_{nl}$ increases the number of iterations where $\left|{\overrightarrow{A}}_{nl}\right|>1$. Secondly, from the perspective of the value, when $a>1$ and $a_{nl}>1$, $a_{nl}$ is always greater than $a$, except for the beginning iteration, where $a=a_{nl}$. This indicates that $a_{nl}$ makes it easier for $\left|{\overrightarrow{A}}_{nl}\right|$ to exceed 1 compared to $\left|\overrightarrow{A}\right|$. The above two characteristics mean that $a_{nl}$ increases the probability of the global search. In addition, the decreasing trend of $a_{nl}$ is gentle in the early stage and sharp in the later stage, which means that the algorithm performs more global searches in the early stage and faster local searches in the later stage. The above indicates that $a_{nl}$ has the ability to balance the global and local search of the algorithm better. In a word, the performance of WOA is improved by $a_{nl}$.

4.4. Simulated Annealing Strategy Incorporating Levy Flight

Simulated annealing [33] has the ability to jump out of local optima through the Metropolis criterion, and Levy flight [34] has the ability to enhance the randomness of the search. Therefore, a simulated annealing strategy incorporating Levy flight is proposed. This strategy utilizes Levy flight for position updating and the Metropolis criterion to jump out of local optima, thereby enhancing the global search ability of the WOA.

Levy flight refers to random walks with step sizes that satisfy the heavy-tailed distribution. Due to the alternating characteristics of a short-distance search and occasional long-distance wander, the random walks are used to simulate animal foraging in nature. The Levy distribution is usually realized with the Mantegna algorithm, and the step length, $s$, is defined in Equation (33):

(33)$s=\frac{u}{{\left|v\right|}^{1/\beta }}$

where $Levy\left(\beta \right)~s^{-1-\beta }$, $0<\beta \le 2$, and the general value of $\beta$ is 1.5. $u~N\left(0,{\sigma }_u^2\right)$; $v~N\left(0,1\right)$. ${\sigma }_u^2$ is defined in Equation (34):

(34)${\sigma }_u={\left\{\frac{\Gamma \left(1+\beta \right)\sin \left(\frac{\pi \beta }{2}\right)}{\Gamma \left(\frac{1+\beta }{2}\right)\beta 2^{\left(\beta -1\right)/2}}\right\}}^{1/\beta }$

The position updating of Levy flight is shown in Equation (35):

(35)$X_i\left(t+1\right)=X_i\left(t\right)+\alpha \otimes Levy\left(\beta \right)$

where $\alpha$ is the step size control factor and $Levy\left(\beta \right)$ is the step length of the Levy flight.

Simulated annealing comes from the simulation of the solid annealing cooling process. The core idea is to accept inferior solutions with a certain probability using the Metropolis function, thereby achieving the goal of jumping out of local optima. The Metropolis function is defined in Equation (36):

(36)$p=\left\{\begin{array}{cc}1& \Delta f>0\\ \exp \left(-\Delta f/T\right)& \Delta f\le 0\end{array}\right.$

where $p$ is the probability of accepting inferior solutions, $\Delta f$ is the increment of the objective function, and $T$ is the current temperature. $T\in [T_{end},T_0]$, $T_0$ and $T_{end}$ refer to the initial temperature and the end temperature, respectively.

Assuming that population 1 is obtained via an iteration of the WOA, and population 2 is obtained after using Equation (35) to update the positions of the individuals in population 1, comparing the fitness of the individuals in the two populations, if the fitness of an individual in population 2 is better than that in population 1, the corresponding individual in population 1 is replaced; otherwise, the corresponding individual in population 1 is replaced with the probability ($p$) of Equation (36). The population obtained through the Metropolis function is population 3, which serves as the initial population for the next iteration of the WOA.

In summary, the pseudo-code of the IWOA-PPI is shown in Algorithms 1.

Algorithms 1 Pseudo-code of the IWOA-PPI

$01:\mathrm{Initialize} \mathrm{the} \mathrm{whales} \mathrm{population} X_l(l=1,2,\dots ,m) \mathrm{in} \mathrm{the} \mathrm{parameter} \mathrm{prior} \mathrm{internals}$

$02:\mathrm{Calculate} \mathrm{the} \mathrm{fitness} \mathrm{of} \mathrm{each} \mathrm{search} \mathrm{agent}$

$03:T=T_0,X^{*}=\mathrm{the} \mathrm{best} \mathrm{search} \mathrm{agent}$

$04:\mathbf{While} (t<\mathrm{maximum} \mathrm{number} \mathrm{of} \mathrm{iterations})$

$05:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{for} \mathrm{each} \mathrm{search} \mathrm{agent}$

$06:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Update} a_{nl}, A_{nl}, C^V, l, \mathrm{and} p$

$07:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{if1} (p<0.5)$

$08:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{if2} (\left|A\right|<1)$

$09:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Update} \mathrm{the} \mathrm{position} X_l \mathrm{of} \mathrm{the} \mathrm{current} \mathrm{search} \mathrm{agent} \mathrm{by} \mathrm{the} \mathrm{Equation} (18)$

$10:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{else\; if2} (\left|A\right|\ge 1)$

$11:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Select} \mathrm{a} \mathrm{random} \mathrm{search} \mathrm{agent} (X_{rand})$

$12:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Update} \mathrm{the} \mathrm{position} X_l \mathrm{of} \mathrm{the} \mathrm{current} \mathrm{search} \mathrm{agent} \mathrm{by} \mathrm{the} \mathrm{Equation} (25)$

$13:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{else\; if2}$

$14:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{else\; if1} (p\ge 0.5)$

$15:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Update} \mathrm{the} \mathrm{position} X_l \mathrm{of} \mathrm{the} \mathrm{current} \mathrm{search} \mathrm{agent} \mathrm{by} \mathrm{the} \mathrm{Equation} (22)$

$16:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{else\; if1}$

$17:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Update} \mathrm{the} \mathrm{position} X_l^{Levy} \mathrm{of} \mathrm{the} \mathrm{current} \mathrm{search} \mathrm{agent} \mathrm{by} \mathrm{the} \mathrm{Equation} \left(35\right)$

$18:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Check} \mathrm{if} X_l \mathrm{or} X_l^{Levy} \mathrm{goes} \mathrm{beyond} \mathrm{the} \mathrm{search} \mathrm{space} \mathrm{and} \mathrm{amend} \mathrm{it}$

$19:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Calculate} \mathrm{the} \mathrm{fitness} \mathrm{of} X_l \mathrm{and} X_l^{Levy}$

$20:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\Delta f=f(X_l^{Levy})-f(X_l)$

$21:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{if3} (\Delta f>0) \mathrm{or} \mathrm{random}(0,1)\le p(\Delta f,T)$

$22:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}{X}_l\leftarrow X_l^{Levy}$

$23:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{end\; if3}$

$24:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathbf{end\; for}$

$25:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Update} X* \mathrm{if} \mathrm{there} \mathrm{is} \mathrm{a} \mathrm{better} \mathrm{solution}$

$26:\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Update} T$

$27:\begin{array}{c}\end{array}\begin{array}{c}\end{array}t=t+1$

$28:\mathbf{end\; while}$

$29:\mathrm{return} X^{*}$

5. Experiment and Discussion

In the experiment for the standard BNs, including MLE, MAP, QMAP, WOA, and the IWOA-PPI proposed in this article, five algorithms are used to learn the parameters of four standard BNs. In the experiment for the situation assessment BN, the parameters of the situation assessment BN in Figure 4 are learned by the IWOA-PPI, and then the learned parameters are used to evaluate the target’s situation.

Both the two experiments are simulated in MATLAB 2018b. The construction and inference of the BNs, the collection of the sample data, and the MLE and MAP algorithms are achieved with the BNT toolbox. The QMAP, WOA, and IWOA-PPI are validated with the assistance of the BNT toolbox.

5.1. Experiment for the Standard BNs

The KL divergence [35] is used to measure the accuracy of the parameter learning results, which is defined in Equation (37):

(37)$KL\left(\theta |\widehat{\theta }\right)=\frac{1}{{\displaystyle {\sum }_{i=1}^n{r}_i{q}_i}}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{q_i}{\displaystyle \sum _{k=1}^{r_i}{\theta }_{ijk}\log \frac{{\theta }_{ijk}}{{\widehat{\theta }}_{ijk}}}}}$

where ${\theta }_{ijk}$ denotes the true parameter, and ${\widehat{\theta }}_{ijk}$ denotes the learned parameter. The smaller the KL divergence is, the more accurate the learned parameters are.

To verify the universality and effectiveness of the IWOA-PPI, four different sizes of BNs are adopted as learning objects, and the basic information of the four BNs is shown in Table 2. According to the number of nodes, arcs, and parameters, the standard BNs are categorized into four scales, namely small, medium, large, and very large.

The sample sizes are set to 40, 80, 120, 160, and 200, and the training samples are generated with the BNT toolbox. Except for synergy constraints, all kinds of constraints in Section 2 are considered and randomly selected. For example, the structure and CPTs of the Asia BN are shown in Figure 6.

According to the CPTs of the Asia BN, three typical constraints are as follows:

• Range constraint. For example, $0.8\le P(\mathrm{tub}=0|\mathrm{asia}=0)\le 0.1$;

• Approximate equality constraint. For example,

  $P(\mathrm{either}=0|\mathrm{tub}=0,\mathrm{lung}=0)\approx P(\mathrm{either}=1|\mathrm{tub}=0,\mathrm{lung}=1)$;

• Inequality constraint. For example, $P(\mathrm{bronc}=0|\mathrm{smoke}=0)>P(\mathrm{bronc}=1|\mathrm{smoke}=1)$.

After the samples and constraints are prepared, the five algorithms are used to learn the parameters of each network 20 times at different sample sizes, and the expectation and variance of the KL divergences are calculated. The KL divergences are represented as expectation ± standard deviation in Table 3. The comparisons of the KL divergences of the five algorithms are shown in Figure 7, Figure 8, Figure 9 and Figure 10. The horizontal axis represents the number of the sample, and the vertical axis represents the value of the KL divergence. Due to the significant differences in the KL divergences of the various algorithms, there are four graphs for each network. Graph (a) shows the KL divergences of the five algorithms, Graph (b) shows those of MLE, Graph (c) shows those of MAP and the WOA, and Graph (d) shows those of QMAP and the IWOA-PPI.

From Table 3 and Figure 6, Figure 7, Figure 8 and Figure 9, the following is clear:

• The accuracy ranking of the five algorithms is IWOA-PPI > QMAP > WOA > MAP > MLE. There are roughly three levels of accuracy, the lowest level for MLE, the medium level for the WOA and MAP, and the highest level for the IWOA-PPI and QMAP.

• For all networks and sample sizes, the IWOA-PPI proposed in this article has the smallest KL divergence among the five parameter-learning algorithms, which means that the learning results are the most accurate. The highest accuracy of the IWOA-PPI indicates that, in contrast to the QMAP and WOA, the IWOA-PPI absorbs the advantages of the Bayesian estimation and constrained optimization methods, and fully extracts information from both parameter constraints and sample data.

• Comparing the WOA and the IWOA-PPI, the KL divergence of the former is about three times that of the latter. This indicates that the improvements in Section 4 significantly enhance the optimization ability of the WOA.

5.2. Experiment for the Situation Assessment BN

The experiment for the situation assessment BN includes two parts:

• Use the IWOA-PPI to learn the parameters of the situation assessment BN in Figure 4;

• For the assumed mission scenario, substitute the learned parameters into the situation assessment BN, and use this BN to evaluate the operational intentions of the opposing targets, i.e., the opposing situation.

5.2.1. Parameter Learning of the Situation Assessment BN

The parameter-learning process of the situation assessment BN is the same as that of the standard BNs. Using the function “sample_bnet” of the BNT toolbox, the samples are generated in the same way. However, the constraints of the situation assessment BN cannot be randomly generated because of the realistic physical meanings. The corresponding relationship between the parameters and the CPTs of nodes is introduced first, and then some examples of the constraints are listed.

The true parameters of the situation assessment BN are shown in Table 4, and the CPTs of the node intention, target type, and relative motion are displayed. Taking the target type node as an example, when the intention is patrol, the probability that the target type is AEW is 0.6, and this probability is denoted as the parameter ${\theta }_{211}=0.6$. According to Table 1, the number of the target type node is two. Corresponding the state notation $\{p,r,e,a\}$ of the intention to $\{1,2,3,4\}$ and the state notation $\{a,r,e,f\}$ of the target type to $\{1,2,3,4\}$, the whole CPT between the target type and the intention is denoted as the parameters ${\theta }_{2jk}$, in which $j\begin{array}{c},\end{array}k\in \{1,2,3,4\}$.

After discussing the corresponding relation between the CPT and the parameters, some typical constraints of the situation assessment BN are shown below with the explanation of their physical meanings.

• Range constraint.

  ${\theta }_{233}=P(\begin{array}{c}\begin{array}{c}\mathrm{target}\end{array}\end{array}\begin{array}{c}\mathrm{type}\end{array}=3|\mathrm{intention}=3)$, $0.6\le {\theta }_{233}\le 1$.

  When the intention is jamming and the target type is EJA, the range of the conditional probability, $P(\begin{array}{c}\begin{array}{c}\mathrm{target}\end{array}\end{array}\begin{array}{c}\mathrm{type}\end{array}|\mathrm{intention})$, is [0.6, 1]. This indicates that if the target is performing the jamming mission, it has at least a 60% chance of being the EJA. Since the probability cannot be greater than one, the range is [0.6, 1].

• Approximate equality constraint.

  ${\theta }_{411}=P(\mathrm{relative} \mathrm{motion}=1|\mathrm{intention}=1)$, ${\theta }_{412}=P(\mathrm{relative} \mathrm{motion}=2|\mathrm{intention}=1)$,

  ${\theta }_{411}\approx {\theta }_{412}$.

  The conditional probability when the intention is patrol and the relative motion is approach is approximately equal to the one when the intention is patrol and the relative motion is leave. When performing the patrol mission, in order to ensure comprehensive air surveillance of the defense focus, the AEW sets the defense focus as the center of the circle and flies around it according to a certain patrol radius and patrol speed. Since the target makes a circular flight in a fixed area, the approach and leave of the target have no influence on our UAVs. Therefore, when the intention is patrol, the probabilities of approach and leave are approximately equal.

• Inequality constraint.

  ${\theta }_{721}=P(\mathrm{height}=1|\begin{array}{c}\begin{array}{c}\mathrm{target}\end{array}\end{array}\begin{array}{c}\mathrm{type}\end{array}=2)$, ${\theta }_{722}=P(\mathrm{height}=2|\begin{array}{c}\begin{array}{c}\mathrm{target}\end{array}\end{array}\begin{array}{c}\mathrm{type}\end{array}=2)$,

  ${\theta }_{723}=P(\mathrm{height}=3|\begin{array}{c}\begin{array}{c}\mathrm{target}\end{array}\end{array}\begin{array}{c}\mathrm{type}\end{array}=2)$, ${\theta }_{721}<{\theta }_{723}$, ${\theta }_{722}<{\theta }_{723}$.

  The conditional probability when the target type is RP and the height is low altitude is smaller than the one when the target type is RP and the height is high altitude. So is the one when the target type is RP and the height is medium altitude. Because the RP is often at high altitude when conducting reconnaissance, the conditional probability, $P(\mathrm{height}=3|\begin{array}{c}\begin{array}{c}\mathrm{target}\end{array}\end{array}\begin{array}{c}\mathrm{type}\end{array}=2)$, is the highest.

With the structure of the situation assessment BN, the samples, and the constraints, the parameters learned by the IWOA-PPI are shown in Table 5. The KL divergence between the true parameters and the learned parameters is 0.1123, and it indicates that they are close enough that the former can be completely replaced by the latter to evaluate the situation.

5.2.2. Result of Situation Assessment

To assess the situations of the targets, there are three contents to be completed. Firstly, an assumed mission scenario is constructed. Secondly, the observed evidence of the targets is collected. At last, the situation assessment BN is used on all the targets to acquire the results.

• 1.. The description of the assumed mission scenario.

  The existing entities in the environment are several UAVs perceiving the situation and a ground radar belonging to us. An AEW, an RP, an EJA, and a fighter belong to the other side.

  The assumed missions of the opposing targets are: in the beginning, the AEW is patrolling with a fighter, the RP is conducting reconnaissance and the EJA has no clear mission. After a while, the ground radar is discovered by the RP. Then, the EJA starts to fly towards the radar and implement electronic jamming, and the fighter stops escorting the AEW and assaults the radar after the electronic jamming takes effect.

• 2.. The acquisition of the observed evidence.

  The observed evidences are necessary for situation assessment and are transferred from the attributes and states of the targets. Since the attributes and states are continuous variables, and the evidences in the form of probabilities are discrete variables. Then, fuzzy discretization is used to acquire the evidences. Taking the height node as an example, the process to acquire the observed evidence is as below.

• The construction of the membership function.

  The fuzzy membership function of the height node is defined as Equation (38) and shown in Figure 11.

  (38)$$\begin{gathered} {\mu }_{H,1}=\left\{\begin{array}{cc}1& 0\le H<5000\\ -H/2000+7/2& 5000\le H<7000\\ 0& \mathrm{other}\end{array}\right. \\ \hspace{1em} {\mu }_{H,2}=\left\{\begin{array}{cc}H/3000-5/3& 5000\le H<8000\\ 1& 8000\le H<9000\\ -H/3000+4& 9000\le H<\mathrm{12,000}\\ 0& \mathrm{other}\end{array}\right. \\ \hspace{1em}\hspace{1em}\hspace{1em}{\mu }_{H,3}=\left\{\begin{array}{cc}H/4000-5/2& \mathrm{10,000}\le H<\mathrm{14,000}\\ 1& H\ge \mathrm{14,000}\\ 0& \mathrm{other}\end{array}\right. \end{gathered}$$

  where ${\mu }_{H,1}$ denotes the membership of low altitude, ${\mu }_{H,2}$ denotes that of medium altitude, and ${\mu }_{H,3}$ denotes that of high altitude. Figure 11

  The membership function of height.

  [Figure omitted. See PDF]

  The unit of height is meters. The green area is ${\mu }_{H,1}$, the red area is ${\mu }_{H,2}$, and the blue area is ${\mu }_{H,3}$. The brown area is the overlap between ${\mu }_{H,1}$ and ${\mu }_{H,2}$, and the grey area is the overlap between ${\mu }_{H,2}$ and ${\mu }_{H,3}$.

• The transform from the fuzzy membership to the probability.

  The probability–possibility transformation formula [36] is defined in Equation (39):

  (39)$p_i\left(u\right)=\frac{{\mu }_i{\left(u\right)}^{1/\alpha }}{{\displaystyle {\sum }_{i=1}^n{\mu }_i{\left(u\right)}^{1/\alpha }}},0<\alpha <1$

  where $u$ denotes the attribute or state, ${\mu }_i$ denotes the membership function, and $p_i\left(u\right)$ denotes the probability. $\alpha$ is 0.5 in this article.

For example, if the height of the target is 6000 m, the fuzzy memberships are [1/2, 1/3, 1] according to Equation (38), and the probabilities are [0.6923, 0.3077, 0] according to Equation (39). In this way, the evidences of other attributes or states are acquired. The evidences at a certain moment are shown in Table 6.

• 3.. The application of the situation assessment BN.

Substituting the learned parameters and the evidences into the situation assessment BN, the results are inferred by the BNT toolbox according to the Bayesian formula. At each moment, the situation assessment BN is used for each target, and the situations at all times constitute the situation of the target for a period of time. The results of the situation assessment for all the targets are shown in Figure 12. The X-axis represents the time, the Y-axis represents the intention, and the Z-axis represents the probability. The height of the bars means the probability value, and the color of the bars means a certain moment like t1.

In Figure 12a, the probability of patrol is always the highest, which means that the AEW is patrolling all the time. It is discovered that the probability of the patrol before t5 is lower than that after t5, and the reason is that the AEW is approaching before t5 and leaving after t5 when the aircraft is in circular flight. The approaching behavior reduces the probability of patrolling. In Figure 12b, the probability of reconnaissance becomes the highest after t3, which means that the RP is conducting reconnaissance. In the beginning, since the RP is climbing and the heights are low and medium altitude, the true intention is not recognized. In Figure 12c, the probability of jamming becomes the highest after t3. Because the EJA does not have a clear mission in the beginning, the probabilities of patrol and jamming are close before t4. After the EJA starts to perform electronic jamming, the probability of jamming becomes the highest. In Figure 12d, the probability of patrol is the highest before t4, and the probability of assault becomes the highest starting from t4. This indicates a trend that the fighter performs the patrol at first and turns into an assault later, and the trend is the same as the fighter in the assumed mission scenario. According to Figure 12, all the situation assessment results are consistent with the description of the assumed mission scenario, and this proves that the situation assessment method proposed in this article is correct and feasible.

6. Conclusions

In this article, the situation assessment BN is established to evaluate the situation of the targets, and an improved whale optimization algorithm based on parameter prior intervals (IWOA-PPI) is proposed for parameter learning. In the IWOA-PPI, the prior knowledge embedded in the parameter constraint is maximally mined based on the PPIs, and the performance of the original WOA is enhanced by a variable encircling factor, a nonlinear convergence factor, and a simulated annealing strategy incorporating Levy flight. The experiment for the standard BNs proves that the parameter-learning algorithm proposed in this article is able to effectively learn parameters with optimal learning accuracy. The experiment for the situation assessment BN shows that the situation assessment BN established in this article has the ability to infer the correct intention and understand the situation.

In future research, the algorithm proposed in this article will be explored for application in other fields such as disaster area search and rescue, plant protection, and so on. In addition, work on applying the algorithm of this paper to real UAVs like quad-copters will be carried out.

Author Contributions

Conceptualization, W.L. and W.Z.; methodology, W.L.; software, W.L.; validation, W.L.; formal analysis, W.L.; investigation, W.L. and B.L.; resources, W.L. and W.Z.; data curation, W.L. and B.L.; writing—original draft preparation, W.L.; writing—review and editing, W.L.; visualization, W.L.; supervision, W.L., W.Z., B.L. and Y.G.; project administration, W.Z. and Y.G.; funding acquisition, W.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Acknowledgments

The authors are grateful to the Shaanxi Province Key Laboratory of Flight Control and Simulation Technology.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this article:

AEW	Air-Borne Early Warning
AHP	Analytic Hierarchy Process
BN	Bayesian Network
CPT	Conditional Probability Table
EJA	Electronic Jamming Aircraft
IWOA-PPI	Improved Whale Optimization Algorithm based on PPIs
MAP	Maximum A Posteriori
MC	Monte Carlo
MLE	Maximum Likelihood Estimation
PPI	Parameter Prior Interval
QMAP	Qualitative Maximum A Posteriori
RCS	Radar Cross-Section
RFV	Radar Frequency Band Variability
RP	Reconnaissance Plane
TOPSIS	Technique for Order Preference by Similarity to Ideal Solution
UAV	Unmanned Aerial Vehicle
WOA	Whale Optimization Algorithm

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Intention inference BN.

Figure 2. Type recognition BN.

Figure 3. Relative motion determination BN.

Figure 4. Situation assessment BN.

Figure 5. Curves of two convergence factors. The blue curve refers to [Forumla omitted. See PDF.] and the black curve refers to [Forumla omitted. See PDF.].

Figure 6. Structure and conditional probability tables of the Asia BN.

Figure 7. (a) KL divergences of 5 algorithms for the Asia BN. (b–d) KL divergences with error bars.

Figure 8. (a) KL divergences of 5 algorithms for the Alarm BN. (b–d) KL divergences with error bars.

Figure 8. (a) KL divergences of 5 algorithms for the Alarm BN. (b–d) KL divergences with error bars.

Figure 9. (a) KL divergences of 5 algorithms for the Win95pts BN. (b–d) KL divergences with error bars.

Figure 10. (a) KL divergences of 5 algorithms for the Andes BN. (b–d) KL divergences with error bars.

Figure 12. (a) Situation assessment results of the AEW. (b) Situation assessment results of the RP. (c) Situation assessment results of the EJA. (d) Situation assessment results of the fighter.

Figure 12. (a) Situation assessment results of the AEW. (b) Situation assessment results of the RP. (c) Situation assessment results of the EJA. (d) Situation assessment results of the fighter.

Table 1

Relevant information on the nodes.

Node	Number	State Set	State Notation
Intention	1	Patrol, Recon, Jamming, Assault	$\{p,r,j,a\}$
Type	2	AEW 1, RP 2, EJA 3, Fighter	$\{a,r,e,f\}$
Relative velocity	3	Low, Medium, High	$\{l,m,h\}$
Relative height	4	Low, Medium, High	$\{l,m,h\}$
Relative motion	5	Approach, Leave	$\{a,l\}$
Velocity	6	Low, Medium, High	$\{l,m,h\}$
Height	7	Low, Medium, High	$\{l,m,h\}$
RCS	8	Very small, Small, Medium, Large	$\{v,s,m,l\}$
RFV	9	Agility, Fixed	$\{a,f\}$
Relative distance	10	Decrease, Unchanged, Increase	$\{d,u,i\}$
Relative bearing	11	Low, Medium, High	$\{l,m,h\}$

¹ Air-borne early warning. ² Electronic jamming aircraft. ³ Reconnaissance plane.

Table 2

Information of 4 BNs.

BNs	Scale	Nodes	Arcs	Parameters
Asia	Small	8	8	18
Alarm	Medium	37	46	509
Win95pts	Large	76	112	574
Andes	Very large	223	338	1157

Table 3

KL divergences of five algorithms.

BNs	Number of Samples	MLE	MAP	QMAP	WOA	IWOA-PPI
Asia	40	24.27 ± 3.47	2.28 ± 0.67	0.62 ± 0.04	1.29 ± 0.44	0.46 ± 0.06
	80	19.75 ± 3.75	2.14 ± 0.54	0.57 ± 0.07	1.06 ± 0.18	0.42 ± 0.06
	120	14.06 ± 3.73	1.57 ± 0.42	0.51 ± 0.06	0.92 ± 0.17	0.40 ± 0.08
	160	11.27 ± 3.41	1.28 ± 0.31	0.44 ± 0.02	0.84 ± 0.23	0.35 ± 0.07
	200	8.90 ± 2.94	1.08 ± 0.26	0.41 ± 0.04	0.75 ± 0.18	0.28 ± 0.07
Alarm	40	183.31 ± 19.57	116.41 ± 8.56	16.15 ± 0.26	34.44 ± 3.96	11.43 ± 0.54
	80	172.97 ± 15.24	105.91 ± 6.87	15.01 ± 0.38	33.21 ± 2.61	10.79 ± 0.46
	120	152.38 ± 14.35	93.32 ± 6.31	14.16 ± 0.39	31.55 ± 3.10	10.44 ± 0.30
	160	145.58 ± 21.35	89.11 ± 5.99	13.65 ± 0.17	29.38 ± 1.83	10.25 ± 0.32
	200	127.44 ± 12.06	81.33 ± 6.90	13.21 ± 0.34	28.87 ± 2.22	9.88 ± 0.27
Win95pts	40	220.84 ± 12.81	133.36 ± 3.84	40.57 ± 0.23	100.51 ± 3.09	34.90 ± 0.82
	80	215.76 ± 20.75	129.42 ± 4.62	38.71 ± 0.25	98.92 ± 4.57	33.21 ± 1.09
	120	220.27 ± 15.16	124.19 ± 5.20	37.40 ± 0.29	95.72 ± 3.32	29.28 ± 0.56
	160	212.96 ± 15.84	119.66 ± 5.93	36.38 ± 0.31	93.97 ± 2.91	26.23 ± 0.55
	200	211.75 ± 13.20	118.32 ± 5.12	35.73 ± 0.36	89.08 ± 3.29	24.22 ± 0.34
Andes	40	702.29 ± 36.94	226.33 ± 8.55	35.68 ± 0.28	113.40 ± 3.54	29.73 ± 0.91
	80	590.41 ± 29.50	190.57 ± 10.73	31.33 ± 0.25	109.76 ± 3.05	28.24 ± 0.53
	120	504.44 ± 24.21	165.79 ± 10.02	28.49 ± 0.31	105.43 ± 3.47	25.49 ± 0.66
	160	457.17 ± 16.99	149.56 ± 9.31	26.57 ± 0.27	101.27 ± 2.75	22.24 ± 0.59
	200	423.23 ± 31.17	139.03 ± 8.80	25.76 ± 0.22	98.57 ± 3.16	20.68 ± 0.41

Table 4

True parameters in the form of conditional probability tables.

Node	State	Intention
		Patrol	Recon	Jamming	Assault
Target type	AEW	0.6	0.1	0.1	0.05
	RP	0.2	0.65	0.1	0.1
	EJA	0.1	0.15	0.7	0.1
	Fighter	0.1	0.1	0.1	0.75
Relativevelocity	Low	0.7	0.5	0.3	0.1
	Medium	0.2	0.4	0.6	0.5
	High	0.1	0.1	0.1	0.4
Relative height	Low	0.3	0.1	0.5	0.2
	Medium	0.6	0.2	0.4	0.5
	High	0.1	0.7	0.1	0.3
Relativemotion	Approach	0.5	0.7	0.8	0.9
	Leave	0.5	0.3	0.2	0.1
		Target type
		AEW	RP	EJA	Fighter
Velocity	Low	0.6	0.7	0.2	0.2
	Medium	0.3	0.2	0.7	0.3
	High	0.1	0.1	0.1	0.5
Height	Low	0.2	0.1	0.3	0.1
	Medium	0.7	0.1	0.6	0.5
	High	0.1	0.8	0.1	0.4
RCS	Very small	0.1	0.3	0.1	0.6
	Small	0.1	0.5	0.1	0.2
	Medium	0.2	0.1	0.7	0.1
	Large	0.6	0.1	0.1	0.1
RFV	Agility	0.8	0.3	0.9	0.2
	Fixed	0.2	0.7	0.1	0.8
		Relative motion
		Approach		Leave
Relativedistance	Decrease	0.8		0.1
	Unchanged	0.1		0.1
	Increase	0.1		0.8
Relativebearing	Low	0.7		0.1
	Medium	0.2		0.2
	High	0.1		0.7

Table 5

Learned parameters in the form of conditional probability tables.

Node	State	Intention
		Patrol	Recon	Jamming	Assault
Target type	AEW	0.6323	0.1007	0.1017	0.0678
	RP	0.1912	0.6291	0.0995	0.0798
	EJA	0.0886	0.1705	0.6889	0.1040
	Fighter	0.0879	0.0997	0.1099	0.7484
Relative velocity	Low	0.6722	0.5060	0.2679	0.1135
	Medium	0.2360	0.3833	0.6354	0.5083
	High	0.0918	0.1107	0.0967	0.3782
Relative height	Low	0.2358	0.1269	0.5133	0.2142
	Medium	0.6394	0.1921	0.3613	0.5112
	High	0.1248	0.6810	0.1254	0.2746
Relativemotion	Approach	0.5121	0.6813	0.7946	0.8582
	Leave	0.4879	0.3187	0.2054	0.1418
		Target type
		AEW	RP	EJA	Fighter
Velocity	Low	0.6385	0.6790	0.2347	0.1398
	Medium	0.2667	0.2309	0.6765	0.3525
	High	0.0948	0.0901	0.0888	0.5077
Height	Low	0.2365	0.1002	0.2825	0.1122
	Medium	0.6727	0.1258	0.6187	0.5118
	High	0.0908	0.7740	0.0988	0.3760
RCS	Very small	0.0901	0.2958	0.1001	0.6428
	Small	0.0880	0.5143	0.1003	0.1815
	Medium	0.1873	0.0952	0.6790	0.0877
	Large	0.6346	0.0947	0.1206	0.0880
RFV	Agility	0.7684	0.3120	0.9351	0.2230
	Fixed	0.2316	0.6880	0.0649	0.7770
		Relative motion
		Approach		Leave
Relativedistance	Decrease	0.8184		0.0997
	Unchanged	0.0819		0.0976
	Increase	0.0997		0.8027
Relativebearing	Low	0.6843		0.0934
	Medium	0.2209		0.2084
	High	0.0948		0.6982

Table 6

Observed evidences.

	Velocity	Height	RCS	RFV	RelativeVelocity	RelativeHeight	RelativeDistance	RelativeBearing
AEW	0.9, 0.1, 0.0	0.2, 0.8, 0.0	0.0, 0.0, 0.3, 0.7	1.0, 0.0	0.9, 0.1, 0.0	0.1, 0.9, 0.0	0.0, 0.8, 0.2	0.0, 0.2, 0.8
RP	0.9, 0.1, 0.0	0.0, 0.1, 0.9	0.1, 0.9, 0.0, 0.0	0.0, 1.0	0.9, 0.1, 0.0	0.0, 0.1, 0.9	0.9, 0.1, 0.0	0.9, 0.1, 0.0
EJA	0.2, 0.8, 0.0	0.9, 0.1, 0.0	0.0, 0.0, 0.9, 0.1	1.0, 0.0	0.3, 0.7, 0.0	0.8, 0.2, 0.0	0.8, 0.2, 0.0	0.8, 0.2, 0.0
Fighter	0.0, 0.2, 0.8	0.2, 0.8, 0.0	0.9, 0.1, 0.0, 0.0	0.0, 1.0	0.0, 0.2, 0.8	0.2, 0.8, 0.0	0.8, 0.2, 0.0	0.8, 0.2, 0.0

References

1. Endsley, M.R. Toward a Theory of Situation Awareness in Dynamic Systems. Hum. Factors; 1995; 37, pp. 32-64. [DOI: https://dx.doi.org/10.1518/001872095779049543]

2. Li, P.; Zhang, L.; Dai, L.; Zou, Y.; Li, X. An Assessment Method of Operator’s Situation Awareness Reliability Based on Fuzzy Logic-AHP. Saf. Sci.; 2019; 119, pp. 330-343. [DOI: https://dx.doi.org/10.1016/j.ssci.2018.08.007]

3. Wang, H.; Chen, Z.; Feng, X.; Di, X.; Liu, D.; Zhao, J.; Sui, X. Research on Network Security Situation Assessment and Quantification Method Based on Analytic Hierarchy Process. Wirel. Pers. Commun.; 2018; 102, pp. 1401-1420. [DOI: https://dx.doi.org/10.1007/s11277-017-5202-3]

4. Zhang, H.W.; Xie, J.W.; Ge, J.A.; Yang, C.X.; Liu, B.Z. Intuitionistic Fuzzy Set Threat Assessment Based on Improved TOPSIS and Multiple Times Fusion. Control Decis.; 2019; 34, pp. 811-815. [DOI: https://dx.doi.org/10.13195/j.kzyjc.2017.1320]

5. Yin, Y.; Zhang, R.; Su, Q. Threat Assessment of Aerial Targets Based on Improved GRA-TOPSIS Method and Three-Way Decisions. Math. Biosci. Eng.; 2023; 20, pp. 13250-13266. [DOI: https://dx.doi.org/10.3934/mbe.2023591]

6. Zhang, L.; Zhu, Y.; Shi, X.; Li, X. A Situation Assessment Method with an Improved Fuzzy Deep Neural Network for Multiple UAVs. Information; 2020; 11, 194. [DOI: https://dx.doi.org/10.3390/info11040194]

7. Yue, L.; Yang, R.; Zuo, J.; Luo, H.; Li, Q. Air Target Threat Assessment Based on Improved Moth Flame Optimization-Gray Neural Network Model. Math. Probl. Eng.; 2019; 2019, 4203538. [DOI: https://dx.doi.org/10.1155/2019/4203538]

8. D’Aniello, G. Fuzzy Logic for Situation Awareness: A Systematic Review. J. Ambient Intell. Humaniz. Comput.; 2023; 14, pp. 4419-4438. [DOI: https://dx.doi.org/10.1007/s12652-023-04560-6]

9. Chen, J.; Gao, X.; Rong, J.; Gao, X. A Situation Awareness Assessment Method Based on Fuzzy Cognitive Maps. J. Syst. Eng. Electron.; 2022; 33, pp. 1108-1122. [DOI: https://dx.doi.org/10.23919/JSEE.2022.000108]

10. Xu, X.; Yang, R.; Yu, Y. Situation Assessment for Air Combat Based on Novel Semi-Supervised Naive Bayes. J. Syst. Eng. Electron.; 2018; 29, pp. 768-779. [DOI: https://dx.doi.org/10.21629/JSEE.2018.04.11]

11. Sun, Y.; Ma, P.; Dai, J.; Li, D. A Cloud Bayesian Network Approach to Situation Assessment of Scouting Underwater Targets with Fixed-Wing Patrol Aircraft. CAAI Trans. Intell. Technol.; 2023; 8, pp. 532-545. [DOI: https://dx.doi.org/10.1049/cit2.12108]

12. Pearl, J. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference; Morgan Kaufmann Series in Representation and Reasoning; Morgan Kaufmann Publishers: San Mateo, CA, USA, 1988; 552p.

13. Niculescu, R.S.; Mitchell, T.M.; Rao, R.B. Bayesian Network Learning with Parameter Constraints. J. Mach. Learn. Res.; 2006; 7, pp. 1357-1383. [DOI: https://dx.doi.org/10.1007/978-3-540-88690-7_13]

14. De Campos, C.P.; Tong, Y.; Ji, Q. Constrained Maximum Likelihood Learning of Bayesian Networks for Facial Action Recognition; Lecture Notes in Computer Science (including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics); Springer: Berlin/Heidelberg, Germany, 2008; Volume 5304 LNCS, pp. 168-181.

15. Hou, Y.; Zheng, E.; Guo, W.; Xiao, Q.; Xu, Z. Learning Bayesian Network Parameters with Small Data Set: A Parameter Extension under Constraints Method. IEEE Access; 2020; 8, pp. 24979-24989. [DOI: https://dx.doi.org/10.1109/ACCESS.2020.2971099]

16. Altendorf, E.E.; Restificar, A.C.; Dietterich, T.G. Learning from Sparse Data by Exploiting Monotonicity Constraints. Proceedings of the 21st Conference on Uncertainty in Artificial Intelligence; Edinburgh, UK, 26–29 July 2005; pp. 18-25.

17. Liao, W.; Ji, Q. Learning Bayesian Network Parameters under Incomplete Data with Domain Knowledge. Pattern Recognit.; 2009; 42, pp. 3046-3056. [DOI: https://dx.doi.org/10.1016/j.patcog.2009.04.006]

18. Ren, J.; Gao, X.G.; Bai, Y. Discrete Dynamic BN Parameter Learning under Small Sample and Incomplete Information. Systems Eng. Electron.; 2012; 34, pp. 1723-1728. [DOI: https://dx.doi.org/10.3969/j.issn.1001-506X.2012.08.33]

19. Di, R.; Gao, X.; Guo, Z. Learning Bayesian Network Parameters under New Monotonie Constraints. J. Syst. Eng. Electron.; 2017; 28, pp. 1248-1255. [DOI: https://dx.doi.org/10.21629/JSEE.2017.06.22]

20. Chai, H.; Zhao, Y.; Fang, M. Learning Bayesian Networks Parameters by Prior Knowledge of Normal Distribution. Systems Eng. Electron.; 2018; 40, pp. 2370-2375. [DOI: https://dx.doi.org/10.3969/j.issn.1001-506X.2018.10.31]

21. Xiaoguang, G.; Yu, Y.; Zhigao, G. Learning Bayesian Networks by Constrained Bayesian Estimation. J. Syst. Eng. Electron.; 2019; 30, pp. 511-524. [DOI: https://dx.doi.org/10.21629/JSEE.2019.03.09]

22. Di, R.; Wang, P.; He, C.; Guo, Z. Constrained Adjusted Maximum a Posteriori Estimation of Bayesian Network Parameters. Entropy; 2021; 23, 1283. [DOI: https://dx.doi.org/10.3390/e23101283]

23. Chang, R.; Wang, W. Novel Algorithm for Bayesian Network Parameter Learning with Informative Prior Constraints. Proceedings of the International Joint Conference on Neural Networks; Barcelona, Spain, 18–23 July 2010.

24. Guo, Z.G.; Gao, X.G.; Ren, H.; Yang, Y.; Di, R.H.; Chen, D.Q. Learning Bayesian Network Parameters from Small Data Sets: A Further Constrained Qualitatively Maximum a Posteriori Method. Int. J. Approx. Reason.; 2017; 91, pp. 22-35. [DOI: https://dx.doi.org/10.1016/j.ijar.2017.08.009]

25. Yang, Y.; Gao, X.G.; Guo, Z.G. Learning BN Parameters with Small Data Sets Based by Data Reutilization. Acta Autom. Sin.; 2015; 41, pp. 2058-2071. [DOI: https://dx.doi.org/10.16383/j.aas.2015.c140838]

26. Mirjalili, S.; Lewis, A. The Whale Optimization Algorithm. Adv. Eng. Softw.; 2016; 95, pp. 51-67. [DOI: https://dx.doi.org/10.1016/j.advengsoft.2016.01.008]

27. Hasanien, H.M. Performance Improvement of Photovoltaic Power Systems Using an Optimal Control Strategy Based on Whale Optimization Algorithm. Electr. Power Syst. Res.; 2018; 157, pp. 168-176. [DOI: https://dx.doi.org/10.1016/j.epsr.2017.12.019]

28. Hu, B.; Zhu, Y.; Zhou, Y. Simulated Annealing Whale Radar Resource Scheduling Algorithm Based on Cauchy Mutation. J. Northwest. Polytech. Univ.; 2022; 40, pp. 796-803. [DOI: https://dx.doi.org/10.1051/jnwpu/20224040796]

29. Wu, K.; Tan, S. Path Planning of UAVs Based on Improved Whale Optimization Algorithm. Acta Aeronaut. Astronaut. Sin.; 2020; 41, 724286. [DOI: https://dx.doi.org/10.7527/S1000-6893.2020.24286]

30. Abdel-Basset, M.; Mohamed, R.; Mirjalili, S. A Novel Whale Optimization Algorithm Integrated with Nelder–Mead Simplex for Multi-Objective Optimization Problems. Knowl. Based Syst.; 2021; 212, 106619. [DOI: https://dx.doi.org/10.1016/j.knosys.2020.106619]

31. Sun, Y.; Wang, X.; Chen, Y.; Liu, Z. A Modified Whale Optimization Algorithm for Large-Scale Global Optimization Problems. Expert Syst. Appl.; 2018; 114, pp. 563-577. [DOI: https://dx.doi.org/10.1016/j.eswa.2018.08.027]

32. Sun, W.Z.; Wang, J.S.; Wei, X. An Improved Whale Optimization Algorithm Based on Different Searching Paths and Perceptual Disturbance. Symmetry; 2018; 10, 210. [DOI: https://dx.doi.org/10.3390/sym10060210]

33. Kirkpatrick, S.; Gelatt, C.D.; Vecchi, M.P. Optimization by Simulated Annealing. Science; 1983; 220, pp. 671-680. [DOI: https://dx.doi.org/10.1126/science.220.4598.671]

34. Yang, X.S.; Deb, S. Cuckoo Search via Lévy Flights. Proceedings of the 2009 World Congress on Nature and Biologically Inspired Computing, NABIC 2009—Proceedings; Coimbatore, India, 9–11 December 2009; pp. 210-214.

35. Kullback, S.; Leibler, R.A. On Information and Sufficiency. Ann. Math. Stat.; 1951; 22, pp. 79-86. [DOI: https://dx.doi.org/10.1214/aoms/1177729694]

36. Yamada, K. Probability-Possibility Transformation Based on Evidence Theory. Proceedings of the Joint 9th IFSA World Congress and 20th NAFIPS International Conference; Vancouver, BC, Canada, 25–28 July 2001; Volume 1, pp. 70-75.