Estimation of Observer Parameters for Dynamic Positioning Ships

Abstract

Considering the problem of dynamic positioning systems for the slowly varying disturbances, a parametrically adaptive observer is presented. The peak frequency of observer is adjusted on-line by autoregressive (AR) spectral estimation; other parameters of observer are optimized using particle swarm optimization (PSO). The peak frequency can be calculated by spectral analysis of the pitch, roll, and heave measurements. In the spectral estimation, Levinson-Durbin algorithm is used to solve the Yule-Walker equations. Finally, the computer simulation is given to demonstrate the effectiveness of the proposed method.

1. Introduction

Dynamic positioning systems keep floating structures in fixed position or predetermined track for marine operation purposes exclusively by means of active thrusters [1].

Filtering and state estimation are important features of a DP system. In most cases, measurements of the vessel velocities are not available. Hence, estimates of the velocities must be computed from noisy position and heading measurements through a state observer. Unfortunately, the position and heading measurements are corrupted with colored noise due to wind, waves, and ocean currents as well as sensor noise. However, only the slowly varying disturbances should be counteracted by the propulsion system, whereas the oscillatory motion due to the waves (first-order wave disturbances) should not enter the feedback loop. This is done by using the so-called wave filtering techniques, which separate the position and heading measurements into a low-frequency and wave frequency position and heading estimate.

The conventional observer filters out the WF motions from the measured position and estimates the LF position and velocity. In the early studies, Balchen et al. [2] and Sørensen et al. [3] used the Kalman filter to filter the WF motion. Later, Fossen and Strand [4] introduced the nonlinear passive observer. It should be noted that the two types of observers are based on a priori knowledge of the sea state to filter the WF motions; this means that peak frequency is assumed to be known. However, over a longer time frame, the sea state may change, and therefore peak frequency in general is not known. So, Strand and Fossen [5] improved the nonlinear passive observer with recursively adaptive WF filtering; however this design has two main drawbacks. Firstly, only the wave model parameters (not the observer gains) are adapted. This means that a priori knowledge of the sea state is required to choose observer gains, such that a notch filtering effect is achieved to remove the first-order wave components. Secondly, the observer tuning can be quite difficult. Guttorm Torsetnes gives a gain scheduled observer. The gain scheduled wave filtering is achieved by measuring the slowly varying wave model parameters on-line and therefore requires no a priori knowledge of the sea state.

In this study, particle swarm optimization (PSO) is implemented in conjuncture with an AR spectrum estimation to adajust the parameters of observer. Some numerical simulation results are given to demonstrate the effectiveness of the proposed method.

2. Ship Model

A useful model describing the dynamics of a surface ship sailing in a horizontal plane having 3 degrees of freedom is given in [1], which is common to separate the model into a kinematic model and dynamic model.

2.1. Kinematic Model

In dynamic positioning, the motions and state variables of the control system are defined and measured with respect to some reference frames or coordinate systems as shown in Figure 1. The Earth-fixed reference frame is denoted by 𝑂_𝐸−𝑋_𝐸𝑌_𝐸𝑍_𝐸, and the body frame is 𝑂−𝑋𝑌𝑍 (see [1]).

Figure 1 

Earth-fixed, reference-parallel, and body-fixed frame.

The state vector is defined by 𝜼=[𝑥,𝑦,𝜓]^𝑇, (𝑥,𝑦)∈ℜ² is the position of the ship given in an inertial frame, and 𝜓∈[0,2𝜋) is the heading angle of the ship relative to geographic North.

If only surge, sway, and yaw (3-DOF) are considered, the kinematics and the state vectors are reduced to 𝜼̇=𝐑𝜓𝐯, (1) where 𝐑(𝜓) is defined as 𝐑𝜓=cos𝜓−sin𝜓0sin𝜓cos𝜓0001. (2)

2.2. Dynamic Model

In the mathematical modeling of ship dynamics, it is common to separate the model into a low-frequency model and wave-frequency model. The WF motion of the ship is due to first-order wave loads. The nonlinear LF equation of motion is driven by second-order mean and slowly varying wave, current and wind loads, and thrust forces.

2.2.1. Low-Frequency Model

The equations of motion of a large class of surface ships can be described by the following model: 𝐌𝐯̇+𝐃𝐯=𝝉+𝐑^𝑇𝜓𝐞. (3)

𝝉=[𝜏_𝑥,𝜏_𝑦,𝜏_𝜓]^𝑇 is the control vector consisting of forces and moments produced by the thruster system; 𝐞∈ℜ^3×1 is the environmental disturbance; 𝐌∈ℜ^3×3 is the system inertia matrix including added mass; 𝐃∈ℜ^3×3 is the damping matrix.

2.2.2. Wave-Frequency Model

In the controller design synthetic white-noise-driven processes consisting of uncoupled harmonic oscillators with damping will be used to model the WF motions. The synthetic WF model can be written in state-space form according to 𝝃̇=𝐀_𝑤𝝃+𝐄_𝑤𝐰_𝑤,𝜼_𝑤=𝐂_𝑤𝝃. (4)𝜼_𝑤=[𝑥_𝑤,𝑦_𝑤,𝜓_𝑤]^𝑇 is the position and orientation measurement vector and 𝜉∈ℜ⁶, 𝐰_𝑤∈𝑅^3×1 is a zero-mean Gaussian white noise vector. Consider 𝐀_𝑤=𝟎_3×3𝐈_3×3𝐀_𝑤²¹𝐀_𝑤²²,𝐄_𝑤=𝟎_3×3𝚺,𝐂_𝑤=𝟎_3×3𝐈_3×3,𝐀_𝑤²¹=diag−𝜔²_𝑜1−𝜔²_𝑜2−𝜔²_𝑜3,𝐀_𝑤²²=diag−2𝜉₁𝑤_𝑜1−2𝜉₂𝑤_𝑜2−2𝜉₃𝑤_𝑜3,𝚺=diag𝑘₁𝑘₂𝑘₃. (5)

The relative damping ratio 𝜉_𝑖 will typically be in the range [0.05–0.1]; 𝜔_𝑜𝑖 is the parameter vector of the unknown dominating wave frequencies in surge, sway, and yaw.

The bias model may also be modelled as 𝐞̇=𝐄_𝑏𝐰_𝑏. (6)

3. Observer Design

Fossen and Strand have proposed a passive nonlinear observer; see [4]. The observer includes wave filtering properties, bias state estimation, reconstruction of the LF motion components, and noise-free estimates of the nonmeasured vessel velocities. The proposed observer has been proven to be passive and globally exponentially stable (GES). Based on (1), (3), (4), and (6) the observer takes the following form: 𝝃̇=𝐀_𝑤𝝃+𝐊₁̇𝜓𝑦,𝜼=𝐑_𝑦𝐯+𝐊₂̇̂𝐲,𝐞=−𝐓_𝑏⁻¹𝐞̂+𝐾₃𝐌̇𝐲,𝐯=−𝐃𝐯+𝐑^𝑇𝜓_𝑦𝐞̂+𝝉+𝐑^𝑇𝜓_𝑦𝐊₄𝐲,𝐲=𝜼+𝐂_𝑤𝝃,𝐯̇^𝑇, (7) where 𝐊₁∈ℜ^6×3, 𝐊₂, 𝐊₃ and 𝐊₄∈ℜ^3×3 are observer gain matrices; 𝐲=𝐲−𝐲 is the estimation error. Consider 𝐊₁=diag𝑘₁,𝑘₂,𝑘₃diag𝑘₄,𝑘₅,𝑘₆,𝐊₂=diag𝑘₇,𝑘₈,𝑘₉,𝐊₃=diag𝑘₃₁,𝑘₃₂,𝑘₃₃,𝐊₄=diag𝑘₄₁,𝑘₄₂,𝑘₄₃.𝑘_𝑖=−2𝜁_𝑛𝑖−𝜁_𝑖𝜔_𝑐𝑖𝜔_𝑖𝑘,𝑖=1,2,3,_𝑖=2𝜔_𝑖𝜁_𝑛𝑖−𝜁_𝑖𝜔_𝑐𝑖𝜔_𝑖,𝑖=4,5,6,𝜔_𝑐𝑖=1.22𝜔_𝑖,𝑘_𝑖=𝜔_𝑐𝑖𝑖=7,8,9, (8)(9) where 𝜁_𝑖 is the desired damping ratio (typically 1) and 𝜔_𝑖 is often set equal to the wave peak frequency. 𝜔_𝑐𝑖=1.22𝜔_𝑖 is the cutoff frequency of the notch filter. 𝜁_𝑛𝑖>𝜁_𝑖 is a tuning parameter to be set between 0.1 and 1.0. 𝐊₃ and 𝐊₄ should be sufficient high to ensure proper bias estimation. 𝐊₃ and 𝐊₄ are simply chosen so that 1/𝐓_𝑏𝑖<𝐊_3𝑖/𝐊_4𝑖<𝜔_𝑖<𝜔_𝑐𝑖 is fulfilled. We can see that 𝜔_𝑖, 𝐊₃, and 𝐊₄ are the unknown parameters in the wave model.

4. Measuring 𝜔_𝑖 by AR Estimation

The dominating wave frequency 𝜔_𝑖 is assumed to be the peak frequency of the unbiased Power Density Spectrum (PDS) based on pitch, roll, and heave response of the vessel. At this point, it is worth noting that the pitch, roll, and heave motions are usually not controlled in DP systems for ships, and they will therefore not be affected by the controller in the same degree as the three horizontal motions; see [6]. The power spectrum estimation method used here is the parametric method of Autoregressive (AR) Spectral Estimation, which has several advantages compared to nonparametric methods; see [7]. Parametric methods are recommended since they give better estimates for shorter finite-length data records.

An autoregressive (AR) model is a possible method of modeling a stationary stochastic process. Given the signal 𝑥(𝑛), 𝑥(𝑛−1), 𝑥(𝑛−2),…,𝑥(𝑛−𝑝), the autoregressive model can be written as 𝑥(𝑛)=𝑤(𝑛)−^𝑝_𝑘=1𝑎_𝑘𝑥(𝑛−𝑘), (10) where 𝑤(𝑛) is an independent identically distributed stochastic sequence with zero mean and variance. The coefficients 𝑎_𝑘 are called the autoregressive parameters of the model, and is 𝑝 the model order. Once such a model of the process is available, its spectrum can be easily computed: 𝑃_𝑥𝑥𝑒^𝑗𝑤=𝜎²_𝑤1+∑^𝑝_𝑘=1𝑎_𝑘𝑒^−𝑗𝑤2. (11)

The values of the parameters 𝑎_𝑘 and 𝜎²_𝑤 are unknown. The problem is to estimate the unknown parameters. 𝑎_𝑘 and 𝜎²_𝑤 can be calculated by solving the following: 𝜙_𝑥𝑥(𝑚)[𝑥(𝑛)𝑥(𝑥(𝑛)=𝐸𝑛+𝑚)]=−𝐸^𝑝_𝑘=1𝑎_𝑘𝑥(𝑛−𝑘+𝑚)+𝑤(𝑛+𝑚)=−^𝑝_𝑘=1𝑎_𝑘𝜙_𝑥𝑥(𝑚−𝑘)+𝐸[𝑥(𝑛)𝑤(𝑛+𝑚)]. (12)

Since 𝐸[𝑥(𝑛)𝑤(𝑛+𝑚)]=_𝜎^0𝑚>02_𝑤𝑚=0, we have 𝜙_𝑥𝑥(𝑚)=−∑^𝑝_𝑘=1𝑎_𝑘𝜙_𝑥𝑥()−𝑚−𝑘,𝑚>0,_𝑝𝑘∑₌₁𝑎_𝑘𝜙_𝑥𝑥(−𝑘)+𝜎²_𝑤,𝑚=0. (13)

From (12) and (13), we obtain 𝜙_𝑥𝑥(0)𝜙_𝑥𝑥(−1)⋯𝜙_𝑥𝑥𝜙−𝑝_𝑥𝑥(1)𝜙_𝑥𝑥(0)⋯𝜙_𝑥𝑥−𝜙𝑝−1⋮⋮⋯⋮_𝑥𝑥𝑝𝜙_𝑥𝑥𝑝−1⋯𝜙_𝑥𝑥(0)1𝑎₁⋮𝑎_𝑝=𝜎²_𝑤0⋮0. (14)

Equation (14) is called the Yule-Walker equations and describes the 𝑝 unknown AR coefficients in terms of the 𝑝+1 estimated autocorrelation coefficients. Solving the Yule-Walker equations for 𝑎_𝑘 is termed the autocorrelation method of AR parameter estimation and can be accomplished using a recursive technique such as the Levinson-Durbin algorithm. The Levinson-Durbin algorithm is given by 𝑘_𝑚𝜙=−^𝑥𝑥(𝑚)+∑^𝑚−1_𝑘=1𝑎_𝑚−1(𝑘)𝜙_𝑥𝑥(𝑚−𝑘)𝜎²_𝑚−1,𝑎_𝑚(𝑘)=𝑎_𝑚−1(𝑘)+𝑘_𝑚𝑎_𝑚−1(),𝜎𝑚−𝑘_2𝑚=𝜎²_𝑚−11−𝑘²_𝑚. (15)

The coefficients 𝑎_𝑚(𝑘) are generated recursively from their (𝑚−1)-order values. 𝜎²_𝑚 is the minimum backward prediction errors of 𝑚-order model, and 𝑘_𝑚 is the reflection coefficient.

5. Optimizing K₃ and K₄

The tuning of the parameters is a necessary part of design observer, and the appropriate observer can guarantee the system obtaining fast dynamic response and robust. In this section, the matrics 𝐊₃ and 𝐊₄ of observer are optimised by PSO.

5.1. Overview of the PSO

The particle swarm optimization (PSO) is an evolutionary computation technique developed by Kennedy and Eberhart in 1995; see [8–11]. Similar to the genetic algorithm (GA), the PSO algorithm is an optimization tool based on population, and the system is initialized with a population of random solutions and can search for optima by the updating of generations. It can be described as follows. At each iteration, each particle can adjust its velocity vector, based on its momentum and the influence of its best position as well as the best position of its neighbors, and then compute a new position where the “particle” is to fly to. Supposing that the dimension for a searching space is 𝐽, the total number of particles is 𝑁, the position of the 𝑖th particle can be expressed as vector 𝐥^𝑖=(𝑙^𝑖₁,𝑙^𝑖₂,…,𝑙^𝑖_𝐽)^𝑇; the best position of the 𝑖th particle being searching until now is denoted by 𝐥^𝑖_𝑙best=(𝑙^𝑖₁,𝑙^𝑖₂,…,𝑙^𝑖_𝐽)^𝑇, and the best position of the total particle swarm being searching until now is denoted by vector 𝐥^𝑔_𝑔best=(𝑙^𝑔₁,𝑙^𝑔₂,…,𝑙^𝑔_𝐽)^𝑇; the velocity of the 𝑖th particle is represented by vector 𝐮^𝑖=(𝑢^𝑖₁,𝑢^𝑖₂,…,𝑢^𝑖_𝐽)^𝑇. So, the velocity and position of the 𝑖th particle can be updated as follows: 𝐮^𝑖(𝑡+1)=𝑊𝐮^𝑖(𝑡)+ℎ₁𝑟₁𝐥^𝑖_𝑙best(𝑡)−𝐥^𝑖(𝑡)+ℎ₂𝑟₂𝐥_𝑔best−𝐥^𝑖(𝑡),𝐥^𝑖(𝑡+1)=𝐥^𝑖(𝑡)+𝐮^𝑖(𝑡+1), (16) where 𝑟₁ and 𝑟₂ are random numbers between 0 and 1; ℎ₁ and ℎ₂ are the acceleration constants with positive values; 𝑊 is the weight function. The following weighting function is usually used in (16): 𝑊=𝑊_start−𝑊_start−𝑊_end𝑡_max×𝑡, (17) where 𝑊_start is the initial weight, 𝑊_end is the final weight, 𝑡_max is the maximum iteration number, and 𝑡 is the current iteration number.

5.2. Fitness Function

Adopting the absolute error (ITAE) as the minimum, we can gain fitness function as 𝐽=^∞₀𝛽𝐳₁(𝑡)𝑑𝑡, (18) with 𝛽 being positive constant.

5.3. Optimization Procedure

The PSO algorithm comprises the following steps.

Step 1. Set 𝑡=0 (𝑡 is the iterative number). Initialization of particle swarm is generated randomly as follows: 𝐥^𝑖(0)=(𝑙^𝑖₁(0),𝑙^𝑖₂(0),…,𝑙^𝑖_𝐽(0))^𝑇=(𝑘₃₁,𝑘₃₂,𝑘₃₃,𝑘₄₁,𝑘₄₂,𝑘₄₃)^𝑇; the swarm has 6-dimension particles which can guarantee 𝐊₃ and 𝐊₄ to be positive matrices. 𝐊₃ and 𝐊₄ are used in the observer mentioned in Section 3. The particle velocities are generated randomly in the range 0~1 and as follows: 𝐮^𝑖=(𝑢^𝑖₁(0),𝑢^𝑖₂(0),…,𝑢^𝑖_𝐽(0))^𝑇.

Step 2. Define objective function values of the particles that are evaluated given by (18). To each particle in the particles’ swarm, set its best position 𝐥^𝑗_𝑙best(0) (the fitness function is the least). Finding the least fitness function values in the initializing particles’ swarm, the best position is 𝐥_𝑔best(0).

Step 3. The position and the velocity of each particle are updated using (16). The fitness function values are calculated for the updated positions of the particles.
If 𝐽(𝐥^𝑗(𝑡))≤𝐽(𝐥^𝑗(𝑡−1)),𝑗=1,2…,𝐽, then 𝐥^𝑗_𝑙best=𝐥^𝑗(𝑡); else 𝐥^𝑗_𝑙best=𝐥^𝑗(𝑡−1).
Setting the least fitness function values in the current particles’ swarm is the 𝑔th particle.
If 𝐽(𝐥^𝑔(𝑡))≤𝐽(𝐥_𝑔best(𝑡−1)),𝑔=1,2,…,𝐽,then 𝐥_𝑔best=𝐥^𝑔(𝑡);else 𝐥_𝑔best=𝐥_𝑔best(𝑡−1).

Step 4. If we attain the maximum iteration, exit, and otherwise go to iterate.

6. Simulation Research

In this section, the performance of the observer tuning by PSO and AR spectrum estimation is verified. A supply ship named as S-175 is used as a case study, and the vessel main particulars are shown in Table 1.

Table 1 

Supply vessel main particulars.

For the purpose of comparison, all the simulations deploy the same parameter settings as follows. The ship controlled by PID controller [1] is maintained in the desired position and heading 𝜂_𝑑=[000^∘]^𝑇. The population size of PSO is 50, the maximum velocity is set as 1, the iteration is 200, and the weights are 𝑊_start=0.9 and 𝑊_end=0.4.

6.1. Example 1. Verify the PSO

In order to demonstrate the efficiency of the proposed PSO, optimization observer by PSO is compared with an observer using genetic algorithmic (GA) [12]. The vessel was exposed to the JONSWAP distributed waves with significant wave height 2.0 m and peak frequency 𝜔_𝑖=0.65 rad/s, 𝜍₁=𝜍₂=𝜍₃=0.1. The parameters are shown in Table 2, and the simulation curve is shown in Figure 2.

Table 2 

Parameters of observer.

(a) Simulation curve of observer using GA

(b) Simulation curve of observer using PSO

(c) The measured position of ship of observer using GA and PSO

Figure 2 

Simulation curve of observer using GA and using PSO.

From Figure 2, we can see that the proposed PSO optimization observer can give better results, with a small position and heading bias.

6.2. Example 2. Verify the Parametrically Adaptive Observer

To demonstrate the performance of parametrically adaptive observer in this paper, the proposed parametrically adaptive observer is compared with the observer (7). The aforementioned two observer gains are optimized using PSO. The JONSWAP wave spectrum with significant wave height to be set at 1.2~3.6 m and peak frequency to be set at 0.8~0.65 rad/s is used to imitate sea conditions. The simulation curve is shown in Figure 3.

(a) Simulation curve of observer without AR estimation

(b) Simulation curve of observer using AR estimation

(c) The measured position of ship of observer without and with AR estimation

(d) The estimated peak frequency

Figure 3 

Simulation curve of observer without and with AR estimation.

From Figures 3(a) and 3(b), the observer with AR spectral estimation can avoid the transients caused by the instability of verying sea states. From Figure 3(c), we can clearly see that the ship using observer with AR estimation has better position and heading. It can be seen from Figure 3(d) that the 𝜔_𝑖 estimated by AR method is very close to predefined peak frequency.

7. Conclusions

This paper proposed a parametrically adaptive observer of ship for the slowly varying seas. The peak frequency of observer can be calculated by AR spectral estimation of the pitch, roll, and heave measurements, some parameters of observer are tuned on-line based on peak frequency, and other parameters of observer are optimized off-line using PSO. The simulation results show that the parametrically adaptive observer can avoid the transients caused by the instability of varying sea states and can improve the accuracy of dynamic positioning systems.

Acknowledgment

The paper was supported by the National Nature Science Foundation of China, no. 60775060.