Experimental validation of propulsion fault detection system using moving horizon estimation in quad-maran unmanned vessel

ABSTRACT A quad-maran unmanned vessel is a new type of unmanned surface vessel for the automatic collection of various environmental data in aquaculture fields. We propose a fault detection and control scheme for a possible thrust decrease failure of quad-maran vessels. The fault detection and control are based on moving horizon estimation and model predictive control. The proposed estimation and control scheme is experimentally validated.


Introduction
Fault detection and fault-tolerant control are important aspects of practical control systems. Faults in systems cause a loss of control performance and/or lead to catastrophic damage to the components of control systems. For example, a sudden loss or damage of propellers in multirotor copters may result in crashing into the ground. Faults can occur during operation, and fault detection needs to be performed in real time based on the measurements of a plant. Various fault detection and fault-tolerant control methods have been proposed and applied to many different systems in different industry sectors [1]. One approach to fault detection is to use moving horizon estimation (MHE) [2][3][4], and several MHE-based fault detection/control methods have been reported [5][6][7][8][9]. The MHE approach uses an estimation horizon that moves along with the time evolution, and the estimation is performed by solving an optimization problem.
The authors have been engaged in the control of a new type of unmanned surface vessel, a quad-maran unmanned vessel [10][11][12][13]. This vessel is mainly developed for environmental monitoring [14] in aquaculture fields, and an observation project has been conducted at oyster farms in Nanao Bay, Ishikawa Pref., Japan ( Figure 1); initial observation results are analysed and reported in [15]. There have been problems such as shellfish dying in aquaculture farms recently because of changes in the water environment [16]. The vessel is expected to make it possible to obtain frequent and high-resolution observation data, such as temperature, oxygen level, and salinity, in aquaculture fields. Such data can play important roles in providing highlevel services for aquaculture; for example, predicting the occurrence of hypoxia water mass. Quad-maran unmanned vessels have a feature that is suitable for operation in aquaculture fields; namely, four independently rotating hulls with propellers. This capability allows for good sailing and dynamic positioning ability in aquaculture fields where many obstacles such as buoys and rafts for farming are present. While fault detection and fault-tolerant control are essential features required for the control system of such vessels, they have yet to be implemented. During a sea trial, the propeller of the hull (Figure 2(a)) became entangled in a floating object (Figure 2(b)). The failed propeller could no longer generate thrust, and this failure could occur without significantly affecting rotation speed and motor current. We were thus unable to detect this type of fault simply by observing propeller system data; the motion of the vessel also needs to be taken into account for fault detection. This fault results in control performance deterioration or complete stalling depending on its severity.
In this paper, we propose a fault detection system for quad-maran vessels to identify faults in the propellers. We use a simple model of the vessel during dynamic positioning and formulate the fault detection problem based on MHE. We then develop a controller using model predictive control (MPC) where the vessel model and tuning parameters are adjusted according to the current estimated fault condition. The fault detection and control are established by solving optimization problems, a mixed-integer quadratic program (MIQP), and a quadratic program (QP), in real time. The MHEbased fault detection with the MPC-based controller is experimentally tested and validated using quad-maran vessels.  The contributions of the paper are twofold: (1) Most existing studies on real-time fault detection and fault-tolerant control [1,[5][6][7][8][9] use simulations for evaluation, but there is always a gap between a plant model and an actual plant (such as model mismatch and disturbances). Experimental validation is necessary to asses whether or not the methods actually work in a real environment. We propose a fault detection system and carry out an experiment, serving as an experimental case study in this field, yielding information on implementation and insight from the experimental results. (2) The proposed fault detection scheme adds an essential feature to quad-maran vessels. The feature serves to increase the reliability and safety of the vessel operation in aquaculture fields. This ultimately contributes to the IoT/ICT-assisted advanced aquaculture and fishery industry [17].
The preliminary version of this paper was presented at a conference [18]. Details of the model identification, MPC-based controller, and additional experimental results have been added to this paper.
The paper is structured as follows. In Section 2, we briefly explain the quad-maran vessel and introduce a model of the vessel. In Section 3, we represent possible fault conditions based on the vessel model and propose an MHE-based fault detection system with an MPC-based controller. In Section 4, we show the experimental results of the proposed method.

Overview
A quad-maran unmanned vessel is composed of four hulls and a central structure that joins them ( Figure 3). Each hull can change its orientation individually and has a propeller at its rear. The rotation direction and speed of the propellers can be changed individually. These features contribute to the dynamic positioning capability and good maneuverability during sailing. For high-speed sailing, all the hulls are parallel, as shown in Figure 3(a). When performing dynamic positioning, the hulls take an X form, as shown in Figure 3(b). In the configuration of Figure 3(b), the vessel can move in all directions, making it suitable for positioning at an observation point.
The position of the vessel is obtained from a GNSS (Global Navigation Satellite System) device, and the azimuth can be obtained from a magnetometer. Anemometers can also be installed if necessary. Details of the vessel equipment are found in [19].

Linear state-space model
A model is required for MHE-based fault detection. As a simple model, we use the following linear state-space model for the X configuration in Figure 3(b).
x := p x p yṗxṗy T , u := u 1 u 2 , where x ∈ R 4 is a state vector, u ∈ R 2 is a control input vector (propeller speed reference), and y ∈ R 2 is an observation vector (vessel's position). The vessel's coordinate system is shown in Figure 4. p x indicates the vessel position parallel to the longitudinal axis of hulls No. 1 and No. 3. The p y axis is orthogonal to the p x axis. The two propellers of the two aligned hulls are paired, and the same control input is given to them.

System identification experiment
The parameter values of a, d, e, and h in A c and B c are required. We perform a system identification experiment and identify the parameters a, d, e, h, w 1 , and w 2 in the model, (3) This model comes from (1) with the wind-effect terms. w x [m/s] is the wind speed component aligned with the p x axis. w y is the same but aligned with the p y axis. w x and w y are obtained by decomposing the wind speed vector to each axis, p x and p y . In the model (3), we consider the motion of each axis to be independent, and thus we identified each axis p x and p y independently. In the identification, we use a discrete-valued input signal, taking −5.0 V and 5.0 V alternately every 8.0 s. The identification data on the p x axis are shown in Figure 5. Figure 5(b) shows that the wind speed w x , and it is mostly positive, giving a positive wind thrust. With these data, we estimated the parameters a, e, w 1 by using the MATLAB ssest function (see Appendix A for an excerpt of MATLAB script used in the identification).   For validation, other data not used in the identification were compared with the behaviour of the model with the estimated parameters a, e, w 1 . The comparison results are shown in Figure 6, and the model output resembles the measured data for the first 20.0 s. Even after 20.0 s, the model can reproduce similar behaviour. The same experiment was performed on the p y axis, and the following parameters were obtained: We performed the similar validation test for p y axis and obtained the result shown in Fig. 7, which indicates that the model output resembles the measured data.
Although the same identification method is used for both p x and p y axis models, only the validation result of p x axis shows the off-set deviation (Fig. 6), occurring at around 25 s. We suspect the gust happened somewhere between 20-25 s as the cause of the error. The anemometer used in the system identification experiment has some smoothing property, and it is not suited to measure the sudden change in the wind like gust. When the gust is not captured well in the anemometer, an off-set like error occurs between the measured data and model output; in Fig. 6, an off-set occurs at around 25 s and the amount of the off-set is kept at almost the same level after that.
The identification input in Fig. 5(a) is a singlefrequency pulse, but not that it contains only a single frequency component. The MALAB function pexcit with default tolerance parameters is used to compute the degrees of persistence of excitation of the identification inputs, and they are 16 (u 1 and u 2 ) and 20 (w x and w y ). Thus it is not fully inappropriate for identifying the second order system with these identification inputs.

Model discretization
We discretize the model (1), (2), giving where w t ∈ R 4 and v t ∈ R 2 are newly introduced, and they represent the process and observation noises, respectively. Note that the wind disturbance terms in (3) are not explicitly included in (4), but instead, the model (4) has the general process noise term. Our fault detection method, described in the following sections, uses the model (4).

Representation of fault situations
When a thrust decrease occurs during the dynamic positioning configuration shown in Figure 3(b), this fault can be regarded as a change in the input matrix B in (4). We consider two fault conditions for each propeller pair (u 1 : #1 and #3, u 2 : #2 and #4). We thus have the following three conditions indexed by j = 0, 1, 2 for each propeller pair: (0) Normal: there is no thrust decrease.
(1) Partial fault: One propeller in a pair has a significant thrust decrease. (2) Near-total fault: Both propellers in a pair have a significant thrust decrease.
A model incorporating these fault conditions into (4) is given by where B j ∈ R 4×2 (B 0 = B) is the input matrix corresponding to the fault condition j (= 0, 1, 2), and δ ij t is the 0-1 integer variable that takes 1 when the assumed fault j at time t occurs on the control input i and 0 when it does not. The equality constraint (8) means that only one fault condition can arise for each control input. For example, when δ 10 t = δ 20 t = 1 and δ ij t = 0 (i = 1, 2, j = 1, 2), the model (5) is equivalent to (4), which indicates no faults. When δ 12 t = δ 21 t = 1 and the other 0-1 variables are 0, the fault conditions j = 2 and j = 1 occur on the control inputs u 1 t and u 2 t , respectively.

Fault detection with MHE
We use MHE (Moving Horizon Estimation) to identify the current fault condition in real time based on the data set of the past outputs and inputs over a certain time horizon [5] (see Figure 8). At time step k, the state variables x k−N+1 , . . . , x k and δ ij k−N , . . . , δ ij k−1 are estimated from u and y between time steps k−1 and k−N. We estimate the states and 0-1 fault variables by solving the following optimization problem.
This optimization problem is a mixed-integer quadratic programming (MIQP) problem. The cost function J(x,ŵ,v,δ) is defined as follows: whereδ j t is defined aŝ The matrix Q −1 ∈ R 4×4 is a weight for process noise, R −1 ∈ R 2×2 is a weight for observation noise, π −1 ∈ R 4×4 is a weight for the estimation error of the state at time step k−N, and δw j ∈ R 2×2 is a weight that adjusts the detection sensitivity of the fault condition j. The form of cost function J is typical in MHE problems. The state estimates do not appear explicitly in (10), but they are evaluated because the state-space equation is included in the constraints of the optimization problem. At each time step, the optimization problem (9) is solved, and the fault condition is identified. The optimization problem requires the state variablex k−N at the starting point of the estimation horizon. For the first estimation,x k−N is initialized with an appropriate value, e.g. the zero vector. After the second estimation,x k−N is initialized with a state estimate obtained at the previous estimation, i.e., being k as the current time,x k−N =x k−N|k−1 wherex k−N|k−1 denotes the estimated state for k − N step made at the time k − 1.
A fault condition is estimated at each time step over the estimation horizon. Because there are unmeasurable disturbances and observation noise, there are errors in estimating the fault conditions. For this reason, the estimated results of a single point in the estimation horizon alone cannot determine whether a fault is actually occurring. Therefore, the actual fault is considered to have occurred if a specific fault condition is identified more than a given ratio T h (0 < T h ≤ 1) over the estimation horizon: F ij t is set to 1 if a fault condition j occurs at the control input u i at time t, and 0 if it does not. A fault indicator variable, F i t , is set to j if a fault condition j occurs at the control input u i at time t as follows: For example, if F i t is zero, it indicates that the control input u i is normal at time t. If F i t is 1, it indicates that the control input u i is in a fault condition 1 at time t.
In practice, the threshold T h is considered as a tuning parameter for adjusting the sensitivity of fault detection, and should be chosen from an interval (0.5, 1). T h should be larger than 0.5 in order to ensure that two fault conditions are not detected at the same time, and T h should be less than 1; here 1 is not allowable because 1 indicates that there are no estimation errors in δ ij t over the entire estimation horizon, which is unrealistic in practical environment. There could be a systematic way for choosing appropriate T h in advance if a certain property on the estimation performance under disturbances and model mismatch is available. However, currently no such method is established, and T h needs to be chosen by trial and error.
The discrete-time model (5)- (8) incorporating the fault conditions is in fact described by the form of switched linear systems where σ t ∈ {1, 2, . . . , 9} denotes the discrete-valued variable representing the fault condition. This is because we have three conditions (normal, partial fault, and near-total fault) for each of two input channels, which totals nine fault patterns. Since we assume that a fault is characterized by a change in the B matrix, nine fault patterns can be represented by  B(1), B(2), . . . , B(9). Thus the MHE problem (9) can be alternatively posed as an estimation problem for the switched linear system where we seek to find the continuous-valued state x t and discrete-valued switching variable σ t over the estimation horizon based on the measurements (past outputs and inputs). Theoretical analysis on the moving horizon estimation of switched linear systems has been conducted in [20]; the main result (Theorem 2 in [20]) is that the estimation error of the continuous-valued state at t − N + α (α > 0) is bounded from above under some assumptions. Note that the formula was presented for the switched linear system without input term, which is not directly applicable to our case. Also the estimation performance on the current state x t and discrete-valued variable have not been made clear in [20]. In our paper, the estimation performance of the proposed scheme is not fully supported by the theoretical results; however, we have verified the estimation performance of the current state and 0-1 variable δ ij t by numerical simulations (Appendix B).

Dynamic positioning system with fault detection
We use model predictive control (MPC) for a controller that generates the control input u t such that the vessel is regulated to a reference state x f . In MPC, control inputs over a future time horizon are computed by solving an optimization problem. We use the first component of the computed inputs over the horizon as an actual input.
The model used in MPC is adjusted to the current identified fault condition F i t . We split the matrices B 0 (= B), B 1 , B 2 in (5) as follows: as the input coefficient matrix in the dynamics equation in MPC. This is updated at each sampling time with the newly estimated fault condition F i t . We solve the following optimization problem at each sampling time: where N c − 1 is the length of the prediction horizon and the cost function J(x, u) is defined as follows: The estimated current state by MHE (denoted bȳ x in the above optimization problem) is used as an initial condition. This optimization problem is a quadratic programming (QP) problem. The inequality constraints in (17) are the upper and lower limits of the voltage reference for the propeller rotation speed. The matrix Q MPC ∈ R 4×4 in (18) is a weight for the error from the reference state x f , and R MPC ∈ R 2×2 is a weight for the control inputs. The weight matrices are adjusted to the identified fault condition. Figure 9 shows the overall structure of the MHEbased fault detection system incorporated into the MPC-based dynamic positioning system.

Experimental system and configuration
The experiment was conducted in a pool at Osaka Prefecture University ( Figure 10). The size of the pool is 50 m × 20 m. The average wind speed on the date of the experiment was 2.6 m/s. The dynamic positioning system with fault detection in Figure 9 was implemented on a laptop PC (CPU : Intel(R) Core(TM) i7-8550U 1.80 GHz, RAM : 16.0 GB) and Raspberry Pi 3B on the central structure of the vessel, as shown in Figure 11. The PC performs the numerical optimizations for MHE (MIQP) and MPC (QP). Raspberry Pi collects observation data and applies the control input to the propeller system. Raspberry Pi sends the current position [p x , p y ] to the PC. The PC receives this value and calculates the control inputs [u 1 , u 2 ] for the next time step. Raspberry Pi receives the control inputs and controls the propellers. The PC and Raspberry Pi are connected by Wi-Fi and communicate by TCP socket connection. The programs on the PC and Raspberry Pi are written in Python. The numerical optimizations for MHE and MPC are performed by calling IBM ILOG CPLEX Optimization Studio library functions from Python. The estimation horizon N was set to 14 steps. Because the sample time is 0.5 s, 14 steps are equivalent to 7.0 s. The weights Q −1 , R −1 , π −1 , δw 1 , and δw 2 used in MHE were set as follows: The threshold T h in (12) was set to 0.75. This is chosen through the process of trial and error and observing the fact that 0.75 is a middle value over (0.5, 1). The length N c of the prediction horizon for MPC was set to 20 steps, which is equivalent to 10.0 s. In the no-fault condition, the following weights were used: When the fault is detected, the diagonal component of R MPC corresponding to the failed propeller is changed to 0.01 for partial fault and 0.001 for near-total fault. This choice is intended to make the propeller compensate for the decreased thrust. The reference state x f for MPC was set to be a zero vector.

Experimental results
In this experiment, two fault situations were assumed. One situation assumed in this experiment was that the propeller of the No. 1 hull, corresponding to the control input u 1 , has a fault (fault condition, j = 1). This situation was simulated by artificially stopping the propeller of the No. 1 hull regardless of the value of the control input u 1 . The experimental results are shown in Figure 12. The fault starts at t = 10.0 s and continues. Figure 12(a) shows that after the fault occurs, the vessel moves in the direction of the reference by the thrust of the No. 3 hull propeller. Afterward, p x is regulated to around 0 m. Figure 12(b) shows that the control input u 1 takes 5.0 V to get the vessel back to the origin from t = 0.0 s to t = 53.0 s. Afterward, u 1 is controlled to keep p x at 0 m. Figure 12(c) shows that F 1 t intermittently takes 1. This shows that the fault condition, which corresponds to the propeller fault on hull No. 1 or No. 3, was detected intermittently. Comparing F 1 t (partial fault) to F 2 t (no fault) in Figure 12(c), there is a different tendency between them. The result suggests that whether the fault is actually occurring can be identified based on the relatively long interval data rather than the short interval (or a single point in time) data.
Another situation assumed in this experiment was that the propellers of hulls No. 2 and No. 4, which correspond to the control input u 2 , have a fault (fault condition, j = 2). This situation was simulated by artificially stopping the propellers regardless of the value of the control input u 2 . The experimental results are shown in Figure 13. The fault occurs at t = 10.0 s and continues. Figure 13(a) shows that for a few seconds after the fault occurs, the vessel moves in the positive p y direction because of inertia. Afterward, p y stops increasing at around 17.0 s. As time progressed further, p y began to decrease and went in the opposite direction of the reference. Figure 13(b) shows that the control input u 2 always takes 5.0 V trying to get the vessel back to the origin. Note that our fault scenario here is complete thrust loss on the control input 2, making the control effort useless after the fault occurs at 10 s. Figure 13(c) shows that F 2 t intermittently takes 2. This shows that the fault condition, which corresponds to the propeller fault on hulls No. 2 and No. 4, was detected intermittently. From t = 43.0 s to t = 58.0 s, F 2 t continuously takes 0, and the propellers had been erroneously estimated to be normal. This estimation error was caused by the vessel moving toward the origin because of wind disturbance. Comparing the fault estimation result F 1 t for the normal propellers (#1 and #3 hulls) to F 2 t , the tendency of each is clearly different.

Conclusion
In this study, a fault detection system using MHE in a quad-maran unmanned vessel was proposed. First, the model of the vessel in the X configuration for dynamic positioning was obtained by system identification experiments. Using the obtained model, the MHE-based fault detection system with an MPC-based? controller was experimentally tested, and it was confirmed that the proposed scheme works for identifying the fault. The future challenge is to distinguish behaviours caused by actual faults and disturbances such as wind and waves to increase detection reliability. Also investigating the estimation performance of the MHE from the theoretical point of view is a challenge from the control theory side.