A nonlinear optimal control approach to stabilization of a macroeconomic development model

1.   Introduction

The ocean is the frontier of major scientific researches such as the origin of life, the earth's evolution, and climate change. It is also the most important trade channel under the background of economic globalization ^[1]. Underwater vehicles play an indispensable role as the key equipment for ocean exploration and development ^[2,3,4]. Thrusters are important propulsive devices of an underwater vehicle, working in high-pressure and high-load environments for a long time. Once one or more thrusters of an underwater vehicle fail, the smooth execution of underwater tasks will not be guaranteed. It may even be impossible to recover the underwater vehicle, which will result in a heavy loss. Therefore, it is of urgent need and great significance to carry out research on thruster fault diagnosis to enhance the safety and reliability of underwater thrusters ^[5,6,7,8].

Data-driven fault diagnostic methods have been widely used because of the difficulty in establishing dynamic models of underwater vehicles and their thrusters ^[9,10]. They can detect faults by directly extracting fault features from observed data, without the mathematical model of the object to be diagnosed. Zhu et al. introduced a deep belief network into multisensory information fusion model to identify uncertain and continuously changing fault patterns of deep-sea human occupied vehicle thruster ^[11]. Wang et al. improved the performance of support vector domain description model for fault diagnosis of AUV thrusters by describing the distribution forms and rules of the mapping data in a high dimensional feature space and optimizing the parameters of kernel function ^[12]. Chu et al. proposed a data argumentation method based on random forest regression for fault samples to solve the problem of low classification accuracy caused by extremely imbalanced datasets collected from actual underwater environments ^[13]. More researches which have been developed or methods which may be applied on underwater thrusters can be obtained in ^[14,15,16]. However, the time dependency and cross-correlation in observed data of actual thrusters under different operating conditions are not fully considered because most of the existing data-driven methods only focus on current or voltage at a fixed rotational speed. Additionally, output signals of a thruster contain nonlinear dynamic disturbances caused by operating environments and load changes, etc. It is difficult to rapidly distinguish whether the changes in output signals are caused by thruster faults or disturbances outside, so that the occurrence time of faults cannot be precisely estimated.

The failure process of a thruster generally includes the following stages: normal operation, performance degradation, abnormal output and complete failure. Each stage of the failure process may also contain several unobservable ("hidden") states because of the high dynamicity of the thruster system and its data. HMM is able to establish the relationship between the thruster data and hidden states by extracting intrinsic information from the data. It can model the failure process of the thruster, and is suitable for fault diagnosis of underwater thrusters ^[17]. Zhang et al. extracted fault feature set of mixed domains from bearing vibration signal by principal component analysis, and trained HMMs on part of the reduced feature set ^[18]. Soleimani et al. presented a methodology for fault detection, prediction and isolation of automotive exhaust gas aftertreatment system based on the integration of HMM and Bayesian networks ^[19]. Arpaia et al. proposed a fault detection method exploiting HMMs for fluid machinery without adequate a-priori information about faulty conditions ^[20]. However, the local optima found by manually setting initial model parameters are not always good enough in the training of traditional HMMs.

In this paper, we explore the HMM for fault diagnosis of underwater thrusters. Overall, our contributions are as follows:

1) We design a three-step data preprocessing method, including smoothing output signals in thruster data by Savitzky-Golay filters, aligning each control signal and its output signal via correlation-based latency estimation, and vector-quantizing the data through k-means algorithm.

2) We propose a fault classification method based on HMMs trained by Particle Swarm Optimization (PSO). The HMMs trained by PSO have better parameter estimations and classification performance than vanilla HMMs.

3) We build a fault emulation experimental platform for fault diagnosis of underwater thrusters and verify the effectiveness of the data preprocessing method and fault classification method using the samples collected from the experiments.

2.   Methods

In Section 2, we introduce the theoretical basis of HMM-based thruster fault diagnosis, design the process of preprocessing thruster data, propose the theory of training HMM by PSO, and then describe the framework of the PSO-HMM-based fault diagnostic approach.

2.1.   Fault diagnostic method based on HMM

HMM describes a Markov process with hidden states. It is suitable for the modeling of dynamic system and can be used for pattern classification. An HMM can be described as:

where $\pi$ is the initial state probability matrix, $A$ is the state transition probability matrix, $B$ denotes a set of emission probability distributions.

In Eq (1), $B$ may be a probability matrix or a set of probability density functions, which depends on whether the observations are discrete (Discrete HMM, DHMM) or continuous (Continuous HMM, CHMM). It is generally assumed that the observations of a CHMM follow multiple Gaussian mixture distributions, that is, each hidden state represents a Gaussian mixture model (GMM).

Both types of HMMs can be used to model continuous signals. The signals need to be transformed into discrete observations through vector quantization for modeling with DHMMs, while they can be directly used as observations for modeling with CHMMs. However, the amount of computation of GMM is much larger than that of vector quantization, which will increase the training time and real-time prediction latency. Therefore, DHMMs are used in the proposed fault classification method for lower computational cost and time requirement.

Figure 1 shows the general architecture of an HMM. The random variable $s\left(t\right)$ denotes the hidden state of thruster at time $t$. The state of the thruster at initial time depends on the matrix $\pi$, and the transition probability from each state to another is determined by the matrix $A$. The random variable $o\left(t\right)$ denotes the observed thruster data at time $t$. The distribution of the variable $o\left(t\right)$, one of the distributions in $B$, is only governed by the value of the hidden variable $s\left(t\right)$.

Figure 1.  General architecture of an HMM.

DownLoad: Full-Size Img PowerPoint

2.2.   Data preprocessing

The preprocessing method of thruster data can be divided into three steps. First, the output signals in thruster data are smoothed by Savitzky-Golay filters to reduce the impact of underwater environmental noises on signal precision. Second, each control signal and its output signal are aligned via latency estimation based on correlation to alleviate dynamic latency between them. Third, the processed data are vector-quantized through k-means algorithm to meet the modeling requirement of DHMM.

2.2.1.   Signal smoothing

Output signals in actual data of underwater thrusters usually contain a lot of environmental noises. These noises such as sensor errors may seriously affect the accuracy of data-driven diagnostic methods. Therefore, it is necessary to smooth the output signals first. A Savitzky-Golay filter is a convolution filter that smooths data by fitting successive sub-sets of adjacent data points with a low-degree polynomial via the method of linear least squares ^[21]. The outliers that deviate from the normal trend are pulled back to a certain extent during the fitting process. Suppose that the length of the signal to be smoothed is $N$, the size $W$ of the sub-set is $2M+1$, and the degree $O$ of the polynomial is $P$. The smoothed signal is obtained by:

where $x[n-m]$ is the ${(n-m)}^{th}$ sampling point, $y\left[n\right]$ is the n^th smoothed output value, and $h\left[m\right]$ is the ${m}^{th}$ convolution coefficient which is given by:

where ${a}_{k}$ is the ${k}^{th}$ polynomial coefficient which is given by:

Where $\mathit{\boldsymbol{d}} = (\mathrm{0, 0}, ..., \mathrm{0, 1}, 0, ..., \mathrm{0, 0}{)}^{T}$ is a column vector with the size of (2M + 1, 1), and ${\mathit{\boldsymbol{A}}}^\mathit{\boldsymbol{T}}$ is a matrix with the size of (P + 1, 2M + 1):

It can be seen that the values of convolution coefficients depend only on $M$ and $P$ and is independent of the input samples, that is, the same coefficients will be obtained when the size W and the degree $O$ are fixed. Moreover, both the extent of noise reduction and signal distortion increase as $M$ increases, while they decrease as $P$ increases ^[22]. Therefore, the size $W$ and the degree $O$ should be carefully chosen to achieve a balance between noise reduction and signal distortion.

2.2.2.   Latency estimation

There is a dynamic latency between each control signal and its output signal in thruster data, which will affect the results of on-line data-driven fault diagnosis of underwater thrusters. The output signal in thruster data is composed of rotational speed and current. The latency between the control signal and rotational speed signal is usually large, while the latency between the control signal and current signal is relatively small. We adopt a kind of correlation-based method for latency estimation ^[23]. The method is described as follows:

1) Assume that the latency between the control signal and rotational speed (or current) signal of a certain sample is $t$ (a non-negative integer) sampling periods.

2) Shift the rotational speed (or current) signal by the assumed number of periods $t$.

3) Calculate the correlation coefficient $r(U, W)$ (or $r(U, I)$) between the control signal $U$ and aligned rotational speed (or current) signal $W$ (or $I$) according to Eq (7) (or Eq (8)).

where $\stackrel{̄}{u}$, $\stackrel{̄}{w}$ and $\stackrel{̄}{i}$ is arithmetic mean of control signal, rotation speed signal and current signal, respectively.

1) Repeat the above three steps with different value of $t$ until the maximum correlation coefficient is found.

2) Consider the number of periods $t$ corresponding to the maximum correlation coefficient as the estimate of the latency between the control signal and rotational speed (or current) signal.

Each sample is processed by the above steps to alleviate the latency between the control signal and output signal.

2.2.3.   Vector quantization

The data of underwater thrusters belong to continuous signals. They must be vector-quantized to meet the modeling requirement of DHMM. Vector quantization is a kind of data compression technique with high efficiency and large compression ratio. K-means clustering is one of the classical vector quantization algorithms.

K-means algorithm is a simple clustering method based on distance with fast convergence speed in most cases. It partitions sampling points into $k$ clusters by measuring distance between sampling point and cluster center. The algorithm is as follows: First, $k$ sampling points are randomly selected from all sampling points as cluster centers. All sampling points are then partitioned into the cluster in which the distance between the sampling point and the cluster center calculated by Eq (9) is the smallest. After partition, new center of each cluster is calculated and selected. The above two steps are repeated iteratively until a desired level of convergence. Finally, the vector quantized data are generated by representing the sampling points in each cluster via the symbol of the cluster.

where ${y}_{i}$ denotes the ${i}^{th}$ sampling point, ${u}_{j}$ denotes the ${j}^{th}$ cluster center.

2.3.   PSO-HMM

There are forward/backward algorithm, Viterbi algorithm and Baum-Welch algorithm for the evaluation task, decoding task and learning task in HMMs, respectively. The learning task in HMMs is a problem of parameter estimation. Assuming that there are $Q$ kinds of hidden states and $Z$ kinds of observations in the HMM, the goal is to find the maximum likelihood estimate of the parameters $\lambda = \left(\pi, A, B\right)$ of the HMM given the observed sequences $O$, which is described as:

Where ${\alpha }_{ij}$ is one of the elements in the state transition probability matrix $A$, and ${\beta }_{ij}$ is one of the elements in the emission probability matrix $B$.

Baum-Welch algorithm is usually used to solve the learning task to a certain extent. It is a special case of the expectation-maximization algorithm. Baum-Welch algorithm generally finds a local optimum because of the manual setting of initial model parameters, although it does guarantee the monotonic increasing of the likelihood. In fact, the parameter learning task in HMMs can be regarded as a general optimization problem with constraints. The objective function to be optimized is $\mathit{log}P\left(O\left|\lambda \right.\right)$. PSO, an optimization algorithm based on swarm cooperation, is able to effectively reduce the gap between the found solution and the global optimal solution. It is easy to program the algorithm and implement parallel processing ^[24]. Therefore, we introduce PSO into HMM for better parameter estimation.

First, the constraints of PSO are defined by the equality constraints in Eq (9). Second, n sets of parameters of HMM subject to the constraints are randomly initialized as the particles with their own positions and velocities. The fitness value of each particle represents the log probability that the sequences $O$ are observed given the set of parameters. Third, the current global best-known solution of the whole swarm and current local best-known solution of each particle are determined by calculating fitness values. Then, the position and velocity of each particle is updated via Eq (11). The above two steps are repeated iteratively until a certain level of convergence. Lastly, the particle with the maximum fitness value is considered as the optimal solution of the parameter learning task.

where ${x}_{i}^{k}$ is the position of the ${i}^{th}$ particle during the ${k}^{th}$ iteration, ${v}_{i}^{k}$ is the velocity of the ${i}^{th}$ particle during the ${k}^{th}$ iteration, ${p}_{best\hspace{0.33em}i}^{k}$ is the current local best-known position of the ${i}^{th}$ particle after $k$ iterations, ${g}_{best}^{k}$ is the current global best-known position of the whole swarm after $k$ iterations, ${c}_{1}$ denotes the maximum step size for particles to move towards ${p}_{best\hspace{0.33em}i}^{k}$, ${c}_{2}$ denotes the maximum step size for particles to move towards ${g}_{best}^{k}$, ${r}_{1}$ and ${r}_{2}$ are random numbers in the interval [0.0, 1.0], $w$ denotes the ability of particles to retain their original velocities.

2.4.   Framework of the approach

Figure 2 shows the framework of the thruster fault diagnostic approach built in this study. The whole process of the fault diagnosis is as follows:

Figure 2.  Framework of the PSO-HMM-based fault diagnostic approach.

DownLoad: Full-Size Img PowerPoint

1) A fault emulation experimental platform is built. The experimental samples of the thruster under normal state and three kind of fault states: driver damage, propeller loss and propeller winding, are obtained from fault emulation experiments.

2) The samples are preprocessed by three main steps. First, they are smoothed by Savitzky-Golay filters and segmented into relatively short samples through sliding window. Then, the samples are aligned via latency estimation. Lastly, they are scaled by min-max normalization and vector quantized by k-means algorithm.

3) Four HMMs are trained by PSO on the preprocessed samples of normal state, driver damage state, propeller loss state and propeller winding state, respectively. These four HMMs comprise the whole thruster fault classifier.

4) The current observed thruster data are input into the thruster fault classifier after the same preprocessing. The output probabilities of the data under the four HMMs are calculated by forward algorithm. The state corresponding to the HMM with the maximum output probability is taken as the current state of the thruster.

The thruster fault diagnosis is implemented through the above process.

3.   Experiments

In Section 3, we introduce the details of thruster fault emulation experiments, show the effect of the designed data preprocessing method, and analyze the results of the PSO-HMM-based fault diagnosis. We also discuss the advantages of the proposed approach by comparing with some other methods.

3.1.   Fault emulation experiments of underwater thruster

In order to test the classification performance of the PSO-HMM-based thruster fault diagnostic approach, we choose a thruster commonly used in the industry as the experimental object. Table 1 shows the specifications of the thruster.

Table 1.  Specifications of the thruster.

Parameters	Values
Product model	T561–150
Rated voltage	150 VDC
Rated speed	2600 rpm
Rated power	1100 W
Control signal	-5–5 V

---

 | Show Table

DownLoad: CSV

Figure 3 shows the structure of the experimental platform. The platform is used to collect the data of the thruster under normal state and three kind of fault states. During one experiment, the thruster is arranged on a fixed support which is in a pool. The PC starts sending control signal to the thruster through the STM32 after the system is powered on. The thruster starts working according to control signal. Meanwhile, output signal consisting of rotational speed signal and current signal is collected by the STM32 and current sensor, and stored by the PC. We collected the experimental samples of the four kinds of thruster states in order by emulating normal operation and different faults over the thruster.

Figure 3.  Structure of the experimental system.

DownLoad: Full-Size Img PowerPoint

3.2.   Results of data preprocessing

As described in Section 2.2, we firstly performed smoothing over the output signals in the samples via Savitzky-Golay filters. The degree $O$ was set to 3. The size $W$ was set to 9 and 19 for rotational speed signal and current signal, respectively. Take the sampling points 1–2000 of the thruster under normal state as an example, the results of smoothing the rotational speed signal and current signal are shown in Figure 4. We observed that the noise in the smoothed signal is significantly less than that in the original signal and the signal distortion is not very serious. It shows the effectiveness of this method.

Figure 4.  Results of smoothing the output signal.

DownLoad: Full-Size Img PowerPoint

After smoothing the signals, we segmented the samples into relatively short samples through sliding window. Each short sample consists of 320 three-dimensional sampling points. The three dimensions denote control signal, rotational speed signal, and current signal. Thus, the size of each sample is 320 × 3. There are 100 samples for each thruster state, i.e., a total of 400 samples. We randomly selected 40 samples of each thruster state as test set and used all remaining samples for training data, i.e., a total of 240 training samples and 160 test samples.

Then, we alleviated the time delay between each control signal and its output signal via the correlation-based latency estimation. Take the sampling points 1–320 of normal state as an example, the correlation coefficients between the control signal and rotational speed (or current) signal with assumed latency from 0 to 15 sampling periods are shown in Figure 5. It can be seen that the number of periods $t$ corresponding to the maximum correlation coefficient between the control signal and rotational speed (or current) signal is 2 (or 1). Thus, we shifted the rotational speed signal and current signal by 2 and 1 periods, respectively, to align the data.

Figure 5.  Correlation coefficients with different assumed latencies.

DownLoad: Full-Size Img PowerPoint

Lastly, we performed min-max normalization and vector quantization over the samples. The range of each feature was scaled into the interval [0.0, 1.0]. The number of clusters $k$ was set to 8. Figure 6 shows the 3D view of the sampling points and clustering result. Each color in Figure 6(A) represents a kind of thruster state, while each color in Figure 6(B) represents a cluster. All the sampling points are partitioned into eight clusters. We then generated the vector quantized samples by representing the sampling points in each cluster via the symbol of the cluster. For example, a normalized sampling point (Control: 0.385, Speed: 0.391, Current: 0.126) which belongs to Cluster 4 is transformed into a scalar 4, while a normalized sampling point (Control: 0.869, Speed: 0.631, Current: 0.668) which belongs to Cluster 6 is transformed into a scalar 6.

Figure 6.  3D view of the sampling points and result of k-means clustering.

DownLoad: Full-Size Img PowerPoint

3.3.   Results of fault diagnosis

We trained four HMMs, including HMM of normal state, HMM of driver damage state, HMM of propeller loss state, and HMM of propeller winding state, on the preprocessed training samples of the four thruster states by the proposed PSO-HMM algorithm. The number of the states in HMMs was set to 3. The parameters of PSO are shown in Table 2. We set the values of ${c}_{1}$ and ${c}_{2}$ by experience and tuned the size of the swarm and $w$ by grid search. The explored grid is the cross-product of the size values ranging in [100,200,300, …, 1000] and w values in [0.3, 0.4, 0.5, …, 1.5]. We chose the values displayed in Table 2, considering both performance and computational cost. The iterative processes of the training are shown in Figure 7. All the training processes reach the convergence in about 50 steps.

Table 2.  Parameters of PSO.

Parameters	Value
Maximum iterations	100
Size of the swarm	500
${c}_{1}$	2
${c}_{2}$	2
$w$	0.6
${x}_{min},$ ${x}_{max}$	0, 1
${v}_{min}$, ${v}_{max}$	-0.05, 0.05

---

 | Show Table

DownLoad: CSV

Figure 7.  Training processes of the four PSO-HMMs.

DownLoad: Full-Size Img PowerPoint

We input the test samples into the four HMMs for evaluating the performance of the PSO-HMM-based fault classifier. Figure 8 and Table 3 show the results of the fault diagnosis. Note that we transformed negative infinity into $-1\times 1{0}^{4}$ to plot all the log probability values into the figure. It can be seen that the samples of driver damage fault and propeller loss fault are completely distinguished from the others, while some of the samples of propeller winding fault are diagnosed as normal samples. We conjecture that the confusion is mainly because the propeller winding fault emulated in Section 3.1 is not serious and the fault features are not obvious when the control voltage is not high enough. We also observed that a few of the normal samples are misclassified. The reason may be the loss of information due to the high compression level of the vector quantization.

Figure 8.  Results of the samples of each thruster state.

DownLoad: Full-Size Img PowerPoint

Table 3.  Confusion matrix and accuracy of the fault classification.

	Normal	Driver damage	Propeller loss	Propeller winding	Accuracy
Normal	35	0	4	1	87.5%
Driver damage	0	40	0	0	100%
Propeller loss	0	0	40	0	100%
Propeller winding	12	0	0	28	70%

---

 | Show Table

DownLoad: CSV

In order to prove the advantages of the proposed approach, we trained a tsfresh-SVM (part of the RFR-SVM that is proposed in ^[13]) and a fault classifier consisting of four vanilla HMMs for comparison. The results are shown in Table 4. It can be seen that the accuracy of the PSO-HMM is higher than that of the vanilla HMM on each thruster state except propeller loss state, showing the effectiveness of training parameters of HMM via PSO. Although the accuracy of the tsfresh-SVM on normal state is little higher than that of the PSO-HMM, the performance of the PSO-HMM on propeller winding state is much better than that of the tsfresh-SVM. Moreover, the total accuracy of the PSO-HMM on the test samples is 89.375%, further illustrating the effectiveness of the proposed diagnostic approach for driver damage fault, propeller loss fault, and propeller winding fault.

Table 4.  Comparison of the accuracy of different fault classifiers.

Class of test samples	Quantity	Accuracy
		tsfresh-SVM	vanilla HMM	PSO-HMM
Normal	40	95%	45%	87.5%
Driver damage	40	100%	95%	100%
Propeller loss	40	100%	100%	100%
Propeller winding	40	37.5%	47.5%	70%
All	160	83.125%	71.875%	89.375%

---

 | Show Table

DownLoad: CSV

4.   Conclusions

In this paper, we provided a fault diagnostic approach based on PSO-HMM for underwater thrusters. Four kinds of thruster states, including normal state, driver damage state, propeller loss state and propeller winding state, can be effectively classified using the approach. The experiments show that the PSO-based parameter learning method for HMM improves the accuracy of thruster fault diagnosis by 17.5% compared with the classical method, proving the effectiveness of the method. Furthermore, the proposed data preprocessing method can be a reference for fault diagnosis of underwater thrusters in future because of its usefulness to thruster data for noise reduction and latency alleviation. In the next step, we will pay more attention on increasing the accuracy of classification on a highly imbalanced thruster dataset and the authenticity of fault emulation experiments so that we can finally apply fault classification on real thruster datasets.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (No. 51839004, No. 52171313, and No. 51909039), the High tech ship innovation project (CY04N20) and the Key Laboratory Fund Project 6142215200305.

Conflict of interest

The authors declare there is no conflicts of interest.