# Multiple Model Adaptive Nonlinear Observer of Dynamic Positioning Ship

Authors: Yehai Xie; Xiaogong Lin; Xinqian Bian; Dawei Zhao
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Considering the filtering problem of dynamic positioning (DP) ship for the slowly varying sea state, a multiple model adaptive observer (MMAO) for dynamic positioning ship is presented. The MMAO consists of a bank of nonlinear subobserver and a dynamic weighting signal generator, in which each sub-observer is designed based on different peak frequency of wave spectrum model. To improve the performance of the observer, subobserver using the measurement of position, velocity, and acceleration is used to update the estimated velocity of ship. The observer parameters are optimized using particle swarm optimization (PSO). Finally, the method is verified effective by the computer simulation.


---

## Body

## 1. Introduction

Dynamic positioning (DP) 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 introduced the nonlinear passive observe [4]. 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 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 improved the nonlinear passive observer with recursively adaptive WF filtering [5]; 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 1st order wave components. Secondly, the observer tuning can be quite difficult. Torsetnes et al. give a gain scheduled observer [6] the gain scheduled wave filtering is achieved by measuring the slowly varying wave model parameters online and therefore requires a priori knowledge of the sea state. In [6] the observer gains were parameterized by the wave peak frequencies and spectral analysis techniques were used to estimate the wave spectrum in surge, sway, and yaw from position and heading measurements. This approach is sensitive to measurement noise and may have latency problems because it requires that the samples acquired be buffered to estimate the Power Density Spectrum of the measurement time series. In [7, 8], a multiple model adaptive Kalman filter (MMAK) was designed; a bank of Kalman filters is designed for a finite number of parameter values, each corresponding to a different peak frequency of the assumed wave spectrum model. The Kalman filter is based on a linear model, and ship motion is a nonlinear process; if a linearization ship model is used to design wave filter, the performance of filter will be reduced.In this paper, inspired by previous pioneering work on DP, a multiple model adaptive nonlinear observer is proposed. In this MMAO, a bank of nonlinear observer relies on measurements of the vessel’s position, velocity, and acceleration and is designed for a finite number of parameter values, each corresponding to a different peak frequency of the assumed wave spectrum model. The main emphasis of the paper is on the nonlinear observer and the use of a multiple model scheme for adaptive observer; however, for the sake of completeness, in the numerical simulations an adaptive backstepping controller is used to control the position of the vessel. All observer parameters are optimized using particle swarm optimization. Numerical simulations are given to demonstrate the effectiveness of the proposed method.The organization of this paper is as follows. Section2 gives ship model. In Section 3, the model of observer measuring position, velocity, and acceleration is given. The multiple model adaptive observer is showed in Section 4. Section 5 briefly explains PSO technique, in finding optimal settings for observer. In Section 6, the numerical simulations are given to verify the effectiveness of the proposed method. Summary and conclusion are given in Section 7.


## 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 Figure1. The Earth-fixed reference frame is denoted as the OE-XEYEZE, and the body frame is O-XYZ (see [1]).Figure 1
Earth-fixed, reference-parallel, and body-fixed frame.The state vector is defined byη=[x,y,ψ]T, (x,y)∈ℜ2 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. v=[u,v,r]T is the velocity vector in body frame O-XYZ.If only surge, sway, and yaw (3-DOF) are considered, the kinematics and the state vectors are reduced to(1)η˙=R(ψ)v;R(ψ) is defined as follows:
(2)R(ψ)=[cosψ-sinψ0sinψcosψ0001].


### 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 lst-order wave loads. The nonlinear LF equation of motion is driven by 2nd-order mean and slowly varying wave, current, and wind loads as well as 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)Mv˙+Dv=τ+RT(ψ)b;τ=[τx,τy,τψ]T is the control vector consisting of forces and moments produced by the thruster system; b∈ℝ3×1 if environmental disturbances; M∈ℝ3×3 is system inertia matrix including added mass; D∈ℝ3×3 is 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)p˙w=Apwpw+Epwwpw,ηw=Cwpw;ηw=[xw,yw,ψw]T is the position and orientation measurement vector and pw∈ℝ6, wpw∈ℝ3 is a zero-mean Gaussian white noise vector. The system matrices Apw, Epw, Cw are given by
(5)Apw=[03×3I3×3Aw21Aw22],Epw=[03×3Σ],Cw=[03×3I3×3]T,
where
(6)Aw21=diag[-wo12-wo22-wo32],Σ=diag[k1k2k3],Aw22=diag[-2ξ1wo1-2ξ2wo2-2ξ3wo3].The relative damping ratioξi will typically be in the range [0.05–0.1]; woi is the wave frequency.The bias model may also be modelled as follows:(7)b˙=Ebwb.


## 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 Figure1. The Earth-fixed reference frame is denoted as the OE-XEYEZE, and the body frame is O-XYZ (see [1]).Figure 1
Earth-fixed, reference-parallel, and body-fixed frame.The state vector is defined byη=[x,y,ψ]T, (x,y)∈ℜ2 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. v=[u,v,r]T is the velocity vector in body frame O-XYZ.If only surge, sway, and yaw (3-DOF) are considered, the kinematics and the state vectors are reduced to(1)η˙=R(ψ)v;R(ψ) is defined as follows:
(2)R(ψ)=[cosψ-sinψ0sinψcosψ0001].


## 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 lst-order wave loads. The nonlinear LF equation of motion is driven by 2nd-order mean and slowly varying wave, current, and wind loads as well as 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)Mv˙+Dv=τ+RT(ψ)b;τ=[τx,τy,τψ]T is the control vector consisting of forces and moments produced by the thruster system; b∈ℝ3×1 if environmental disturbances; M∈ℝ3×3 is system inertia matrix including added mass; D∈ℝ3×3 is 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)p˙w=Apwpw+Epwwpw,ηw=Cwpw;ηw=[xw,yw,ψw]T is the position and orientation measurement vector and pw∈ℝ6, wpw∈ℝ3 is a zero-mean Gaussian white noise vector. The system matrices Apw, Epw, Cw are given by
(5)Apw=[03×3I3×3Aw21Aw22],Epw=[03×3Σ],Cw=[03×3I3×3]T,
where
(6)Aw21=diag[-wo12-wo22-wo32],Σ=diag[k1k2k3],Aw22=diag[-2ξ1wo1-2ξ2wo2-2ξ3wo3].The relative damping ratioξi will typically be in the range [0.05–0.1]; woi is the wave frequency.The bias model may also be modelled as follows:(7)b˙=Ebwb.


## 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)Mv˙+Dv=τ+RT(ψ)b;τ=[τx,τy,τψ]T is the control vector consisting of forces and moments produced by the thruster system; b∈ℝ3×1 if environmental disturbances; M∈ℝ3×3 is system inertia matrix including added mass; D∈ℝ3×3 is 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)p˙w=Apwpw+Epwwpw,ηw=Cwpw;ηw=[xw,yw,ψw]T is the position and orientation measurement vector and pw∈ℝ6, wpw∈ℝ3 is a zero-mean Gaussian white noise vector. The system matrices Apw, Epw, Cw are given by
(5)Apw=[03×3I3×3Aw21Aw22],Epw=[03×3Σ],Cw=[03×3I3×3]T,
where
(6)Aw21=diag[-wo12-wo22-wo32],Σ=diag[k1k2k3],Aw22=diag[-2ξ1wo1-2ξ2wo2-2ξ3wo3].The relative damping ratioξi will typically be in the range [0.05–0.1]; woi is the wave frequency.The bias model may also be modelled as follows:(7)b˙=Ebwb.


## 3. Observer Based Measured Position, Velocity, and Acceleration

In [9], a nonlinear observer based model was proposed. It was based on the nonlinear passive observer developed in [4] but extended with the option to utilize both velocity and acceleration measurements in addition to the position and heading measurements.


### 3.1. Ship Model Contains Positions, Velocities, and Accelerations

The by assumption uncorrelated wave induced positions, velocities, and accelerations, respectively, are given by(8)p˙w=Apwpw+Epwwpw,v˙w=Avwvw+Evwwvw,a˙w=Aawaw+Eawwaw,
where pw∈ℝ6, vw∈ℝ6, aw∈ℝ4 describe the first order wave-induced positions, velocities, and accelerations, respectively. Apw∈ℝ6×6, Avw∈ℝ6×6, Aaw∈ℝ4×4 are assumed by Hurwitz and describe the first order wave induced motion. For position, velocity, and acceleration measurements, j=p,v,a, a cascade of second order linear systems
(9)Ajw=[0I-Ωj-Λj],Cjw=[0I],
with
(10)Λj=diag(2ςj,1ωj,1,…,2ςj,nyiωj,nyj),nyp=nyv=3Ωj=diag(ωj,12,…,ωj,nyj2),nya=2,
where ωj,k>0 is the resonance frequency and ζj,k>0 is the relative damping factor which determines the width of the spectrum.In the observer, the position and heading measurements are always required. The number of velocity measurements and acceleration measurements utilized are 3 and 2, respectively. We gety1∈ℝ3, y2∈ℝ3, y3∈ℝ2 and define the measurements:
(11)y=[y1y2y3]=[η+Cpwpwγ2v+Cvwvwγ3v˙+Cawaw],
where γ2=I3×3 and γ3=[100;010]T are the projections extracting the measured velocities and accelerations, respectively, from the actual three DOF velocity and accelerations vectors.Equations (10) and (11) can be rephrased in a standard form as
(12)z˙=TpT(ψy)ApTp(ψy)z+Bqτq+Eqw,y=Cpz+Dpτq.The transformation matrixTp(ψy)∈ℝ25×25 is given by
(13)Tp(ψy)=diag(RT(ψy),…,RT(ψy),I13).The statez∈ℝ25 is given by
(14)z=[pwTηTbTvwTawTvT]T.The system matricesAp∈ℝ25×23, Eq∈ℝ25×12, and Bq∈ℝ3×25, Dp∈ℝ3×8, Cp∈ℝ8×25 are given by
(15)Bq=[036033033036034M-1]T,Ap=[Apw063063066062063036033033036032I33036033-Tb-1036032033066063063Avw062063046043043046Aaw043036033M-1036032-M-1D],Eq=[Epw063063063033033033033033Eb033033063063Evw063043043043Eaw033033033033],Dp=[033033M-Tγ3T]T,Cp=[CpwI033036034033036033033Cvw034γ2026023γ3M-1RT026Caw-γ3M-1D].


### 3.2. Observer Design

By duplicating the model dynamics (14) and introducing a low-pass filter in order to achieve a certain roll-off effect, the following model based observer is proposed:(16)a^˙f=Tf-1(-a^f+y~3),z^˙=TpT(ψy)ApTp(ψy)z^+Bqτq+Kq(ψy)y~+Kpfa^f,y^=Cp(ψy)z^+Dpτq,
where y~ is estimation error. To reduce the number of interconnections, the observer gain matrices are chosen to be
(17)Kpf=[036033033036034KpaT]T,Kq(ψy)=[K11063062K21033032K31033032063K42032043043K53K61RT(ψy)K62032],
where K11∈ℝ6×3, K42∈ℝ6×3, K53∈ℝ4×2, K31∈ℝ3×3, K62∈ℝ3×3, K21∈ℝ3×3, K61∈ℝ3×3are the observer gain matrices. af∈ℝ2×2 is low-pass filter; Kpa∈ℝ3×2, Τf∈ℝ2×2 are the filter constants. K11, K42, K53, and K31 are given the following structure:
(18)K11=[diag{k11,k12,k13}diag{k14,k15,k16}],K42=[diag{k41,k42,k43}diag{k44,k45,k46}],K53=diag{k51,k52k54,k55},K31=diag{k37,k38,k39}.In order to ensure passivity and to relate the observer gains of (18) to the dominating wave response frequencies, it is proposed that
(19)ki=-2(ζni-ζi)ωciωii=11,12,13,41,42,43,51,52,ki=2ωi(ζni-ζi)ωciωii=14,15,16,44,45,46,54,55,ki=ωcii=37,38,39,
where ωci>ωi is the filter cut-off frequency and ωi=ωp is peak frequency. ζni>ζi is a tuning parameter to be set between 0.1 and 1.0. K21, K62, and K61 should be sufficiently high to ensure proper bias estimation.


## 3.1. Ship Model Contains Positions, Velocities, and Accelerations

The by assumption uncorrelated wave induced positions, velocities, and accelerations, respectively, are given by(8)p˙w=Apwpw+Epwwpw,v˙w=Avwvw+Evwwvw,a˙w=Aawaw+Eawwaw,
where pw∈ℝ6, vw∈ℝ6, aw∈ℝ4 describe the first order wave-induced positions, velocities, and accelerations, respectively. Apw∈ℝ6×6, Avw∈ℝ6×6, Aaw∈ℝ4×4 are assumed by Hurwitz and describe the first order wave induced motion. For position, velocity, and acceleration measurements, j=p,v,a, a cascade of second order linear systems
(9)Ajw=[0I-Ωj-Λj],Cjw=[0I],
with
(10)Λj=diag(2ςj,1ωj,1,…,2ςj,nyiωj,nyj),nyp=nyv=3Ωj=diag(ωj,12,…,ωj,nyj2),nya=2,
where ωj,k>0 is the resonance frequency and ζj,k>0 is the relative damping factor which determines the width of the spectrum.In the observer, the position and heading measurements are always required. The number of velocity measurements and acceleration measurements utilized are 3 and 2, respectively. We gety1∈ℝ3, y2∈ℝ3, y3∈ℝ2 and define the measurements:
(11)y=[y1y2y3]=[η+Cpwpwγ2v+Cvwvwγ3v˙+Cawaw],
where γ2=I3×3 and γ3=[100;010]T are the projections extracting the measured velocities and accelerations, respectively, from the actual three DOF velocity and accelerations vectors.Equations (10) and (11) can be rephrased in a standard form as
(12)z˙=TpT(ψy)ApTp(ψy)z+Bqτq+Eqw,y=Cpz+Dpτq.The transformation matrixTp(ψy)∈ℝ25×25 is given by
(13)Tp(ψy)=diag(RT(ψy),…,RT(ψy),I13).The statez∈ℝ25 is given by
(14)z=[pwTηTbTvwTawTvT]T.The system matricesAp∈ℝ25×23, Eq∈ℝ25×12, and Bq∈ℝ3×25, Dp∈ℝ3×8, Cp∈ℝ8×25 are given by
(15)Bq=[036033033036034M-1]T,Ap=[Apw063063066062063036033033036032I33036033-Tb-1036032033066063063Avw062063046043043046Aaw043036033M-1036032-M-1D],Eq=[Epw063063063033033033033033Eb033033063063Evw063043043043Eaw033033033033],Dp=[033033M-Tγ3T]T,Cp=[CpwI033036034033036033033Cvw034γ2026023γ3M-1RT026Caw-γ3M-1D].


## 3.2. Observer Design

By duplicating the model dynamics (14) and introducing a low-pass filter in order to achieve a certain roll-off effect, the following model based observer is proposed:(16)a^˙f=Tf-1(-a^f+y~3),z^˙=TpT(ψy)ApTp(ψy)z^+Bqτq+Kq(ψy)y~+Kpfa^f,y^=Cp(ψy)z^+Dpτq,
where y~ is estimation error. To reduce the number of interconnections, the observer gain matrices are chosen to be
(17)Kpf=[036033033036034KpaT]T,Kq(ψy)=[K11063062K21033032K31033032063K42032043043K53K61RT(ψy)K62032],
where K11∈ℝ6×3, K42∈ℝ6×3, K53∈ℝ4×2, K31∈ℝ3×3, K62∈ℝ3×3, K21∈ℝ3×3, K61∈ℝ3×3are the observer gain matrices. af∈ℝ2×2 is low-pass filter; Kpa∈ℝ3×2, Τf∈ℝ2×2 are the filter constants. K11, K42, K53, and K31 are given the following structure:
(18)K11=[diag{k11,k12,k13}diag{k14,k15,k16}],K42=[diag{k41,k42,k43}diag{k44,k45,k46}],K53=diag{k51,k52k54,k55},K31=diag{k37,k38,k39}.In order to ensure passivity and to relate the observer gains of (18) to the dominating wave response frequencies, it is proposed that
(19)ki=-2(ζni-ζi)ωciωii=11,12,13,41,42,43,51,52,ki=2ωi(ζni-ζi)ωciωii=14,15,16,44,45,46,54,55,ki=ωcii=37,38,39,
where ωci>ωi is the filter cut-off frequency and ωi=ωp is peak frequency. ζni>ζi is a tuning parameter to be set between 0.1 and 1.0. K21, K62, and K61 should be sufficiently high to ensure proper bias estimation.


## 4. Multiple Model Adaptive Observer For DP

The sea state may undergo large variations; therefore, the observer in charge of reconstructing the LF motion should adapt to the sea state itself. In this paper, a multiple model adaptive observer (MMAO) is proposed for dynamic positioning ship under varying seas. Figure2 shows the architecture of a MMAO system. For some early references on Multiple-Model Adaptive Estimator see [7, 8, 10, 11].Figure 2
The MMAO architecture.The structure of MMAO, in Figure2, consists of: (i) the dynamic weighting signal generator (DWSG) and (ii) a bank of observer (22), where each local observer is designed based on one of the representative parameters. The state estimate is generated by a probabilistically weighted sum of the local state-estimates produced by the bank of observer. It is assumed that a linear time-invariant plant G is driven by white noise and a known deterministic input signal and that it generates measurements that are corrupted by white measurement noise. If the plant has an uncertain real-parameter vector, say ωp, one can imagine that it is “close” to one of the elements of a finite discrete representative parameter set, Ω≔{ωp1,ωp2,…,ωpN}. One can then design a bank of standard observer, where each observer uses one of the discrete parameters ωpi, in its implementation, and i∈{1,…,N}. It turns out that if indeed the true plant parameter is one of the discrete values, then the conditional probability density of the state is the sum of Gaussian densities.Consider a finite set of candidate peak frequency valuesΩ≔{ωp1,ωp2,…,ωpN} indexed by i∈{1,…,N}. The state estimate of the MMAO is given by
(20)z^(t)≔∑i=1Npi(t)z^i(t),
where z^(t) is the estimate of the state z(t) (at time t) and pi(t) is the conditional probability that ωp=ωpi, given the measurements record. In (20), each z^i(t) corresponds to a “local” state estimate generated by the ith steady state observer (16).In the proposed MMAO, the dynamic weightspi(t)∈ℝ, i∈{1,…,N} satisfy
(21)p˙i(t)=-λ(1-βi(t)e-mi(t)∑j=1Npj(t)βj(t)e-mj(t))pi(t),
where λ is a positive constant, βi(t) is a signal assumed to satisfy the condition c1≤βi(t)≤c2 for some positive constants c1, c2, and mi(·) is a continuous function called an error measuring functionthat maps the measurable signals of the plant and the states of the ith local observer to a nonnegative real value. The βi(t) and mi(·) are given by
(22)βi(t)≔1detSi,mi(t)≔12∥y(t)-y^i(t)∥Si-12,
where Si is a uniformly bounded positive definite weighted matrix and ∥z∥S=(zTSz)1/2. The matrices Si are important to scale the energy of estimation error signals y(t)-y^i(t) in order to make them comparable. In what follows, we refer to (21) as the dynamic weighting signal generator (DWSG).We impose the constraint that the initial conditionspi(0) be chosen such that pi(0)∈(0,1) and ∑i=1Npi(0)=1. The parameters Q,R and the functions βi,mi are tuning parameters/functions of the MMAO chosen by the designer based on the system being modeled.


## 5. Optimizing the Observer

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 matricesK31, K61, K62, Ka, and Tf of observer are optimised by particle swarm optimization (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 [12–15]). 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 that the “particle” is to fly to. Supposing the dimension for a searching space is J, the total number of particles is N, and the position of the ith particle can be expressed as vector1 li=(l1i,l2i,…,lJi)T; the best position of the ith particle being searched until now is denoted as llbesti=(l1il2i⋯lJi)T, and the best position of the total particle swarm being searched until now is denoted as vector lgbestg=(l1gl2g⋯lJg)T; the velocity of the ith particle is represented as vector ui=(u1iu2i⋯uJi)T. So, the velocity and position of the ith particle can be updated as follows:
(23)ui(t+1)=Wui(t)+h1r1(llbesti(t)-li(t))+h2r2(lgbest-li(t)),li(t+1)=li(t)+ui(t+1),
where r1 and r2 are a random number between 0 and 1; h1 and h2 are the acceleration constants with positive values; W is the weight function. The following weighting function is usually used in (23):
(24)W=Wstart-Wstart-Wendtmax×t,
where Wstart is the initial weight, Wend is the final weight, tmax is the maximum iteration number, and t is the current iteration number.


### 5.2. Fitness Function

Adopting the absolute error (ITAE) as the minimum, we can gain fitness function as follows:(25)J=∫0∞β|y~i(t)|dt,
with β being positive constant and y~i(t) being estimation error of the ith local observer.


### 5.3. Optimization Procedure

The PSO algorithm comprises the following steps.Step 1.
Sett=0 (t is the iterative number). Initialization particles swarm is generated randomly as follows: li(0)=(l1i(0),l2i(0),…,lJi(0))T=(k31,k32,k33,k41,k42,k43)T; the swarm has 6-dimension particles which can guarantee K31, K61, K62, Ka, and Tf to be positive matrices. K31, K61, K62, Ka, and Tf are used in the observer mentioned in Section 3. The particle velocities are generated randomly in the range 0~1 as follows: ui=(u1i(0)u2i(0)⋯uJi(0))T.Step 2.
Define objective function values of the particles are evaluated given by (25). To each particle in the particles swarm, set its best position llbestj(0)   (the fitness function is the least). Finding the least fitness function values in the initializing particles swarm, the best position is lgbest(0).Step 3.
The position and the velocity of each particle are updated using (23). The fitness function values are calculated for the updated positions of the particles.
IfJ(lj(t))≤J(lj(t-1)), j=1,2,…,J, then llbestj=lj(t), else llbestj=lj(t-1).
Setting the least fitness function values in the current particles swarm is thegth particle.
IfJ(lg(t))≤J(lgbest(t-1)), g=1,2,…,J, then lgbest=lg(t), else lgbest=lgbest(t-1).Step 4.
If the maximum iteration is attained, exit; otherwise, go to iterate.


## 5.1. Overview of the PSO

The particle swarm optimization (PSO) is an evolutionary computation technique developed by Kennedy and Eberhart in 1995, (see [12–15]). 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 that the “particle” is to fly to. Supposing the dimension for a searching space is J, the total number of particles is N, and the position of the ith particle can be expressed as vector1 li=(l1i,l2i,…,lJi)T; the best position of the ith particle being searched until now is denoted as llbesti=(l1il2i⋯lJi)T, and the best position of the total particle swarm being searched until now is denoted as vector lgbestg=(l1gl2g⋯lJg)T; the velocity of the ith particle is represented as vector ui=(u1iu2i⋯uJi)T. So, the velocity and position of the ith particle can be updated as follows:
(23)ui(t+1)=Wui(t)+h1r1(llbesti(t)-li(t))+h2r2(lgbest-li(t)),li(t+1)=li(t)+ui(t+1),
where r1 and r2 are a random number between 0 and 1; h1 and h2 are the acceleration constants with positive values; W is the weight function. The following weighting function is usually used in (23):
(24)W=Wstart-Wstart-Wendtmax×t,
where Wstart is the initial weight, Wend is the final weight, tmax is the maximum iteration number, and t is the current iteration number.


## 5.2. Fitness Function

Adopting the absolute error (ITAE) as the minimum, we can gain fitness function as follows:(25)J=∫0∞β|y~i(t)|dt,
with β being positive constant and y~i(t) being estimation error of the ith local observer.


## 5.3. Optimization Procedure

The PSO algorithm comprises the following steps.Step 1.
Sett=0 (t is the iterative number). Initialization particles swarm is generated randomly as follows: li(0)=(l1i(0),l2i(0),…,lJi(0))T=(k31,k32,k33,k41,k42,k43)T; the swarm has 6-dimension particles which can guarantee K31, K61, K62, Ka, and Tf to be positive matrices. K31, K61, K62, Ka, and Tf are used in the observer mentioned in Section 3. The particle velocities are generated randomly in the range 0~1 as follows: ui=(u1i(0)u2i(0)⋯uJi(0))T.Step 2.
Define objective function values of the particles are evaluated given by (25). To each particle in the particles swarm, set its best position llbestj(0)   (the fitness function is the least). Finding the least fitness function values in the initializing particles swarm, the best position is lgbest(0).Step 3.
The position and the velocity of each particle are updated using (23). The fitness function values are calculated for the updated positions of the particles.
IfJ(lj(t))≤J(lj(t-1)), j=1,2,…,J, then llbestj=lj(t), else llbestj=lj(t-1).
Setting the least fitness function values in the current particles swarm is thegth particle.
IfJ(lg(t))≤J(lgbest(t-1)), g=1,2,…,J, then lgbest=lg(t), else lgbest=lgbest(t-1).Step 4.
If the maximum iteration is attained, exit; otherwise, go to iterate.


## 6. Simulation Research

In this section, numerical results are used to demonstrate the proposed MMAO. The simulations are performed using the Marine Systems Simulator (MSS) developed by the Norwegian University of Science and Technology. The MSS incorporates high fidelity models, denoted as process plant model or simulation model in [16], at all levels (plants and actuators). It captures hydrodynamic effects, generalized coriolis and centripetal forces, nonlinear damping and current forces, and generalized restoring forces. It is composed of different modules such as environmental module, vessel dynamics module, thruster and shaft module, and vessel control module. A supply ship named as S-175 is used as case study (Table 1).Table 1
Supply vessel main particulars.
Over length
175.00 m
Breadth
25.40 m
Design draught
9.50 m
M
107  ×  [2.64 0 0; 0 3.34 1.49; 0 1.49 652.09]
D
105  ×  [0.22 0 0; 0 2.22 −17.75; 0 −17.75 715.06]In these simulations, the different environment conditions from calm to high seas are simulated using the spectrum of the Joint North Sea Wave Project (JONSWAP) [2]. The ship is maintained in the desired position and heading ηd=[000°]T.For each candidate value, a nonlinear observer is developed based on the model described with (22) and a MMAO is derived with the dynamic weights given by (28). A nonlinear adaptive backstepping controller is designed that usesη^ and v^ provided by MMAO to control the position of the ship, see [17]. Because the emphasis of this paper is not on control, but rather on filtering, we eschew the details of controller design.Three sets of computational simulations were conducted to test the performance of the proposed MMAO. The first set compares PSO algorithm with the GA algorithm (GA) of [18]. The first set compares the PSO and GA. The second set compares MMAO with the single observer.


### 6.1. Verify the PSO

In order to demonstrate the efficiency of the proposed PSO, optimization observer by PSO is compared with an observer using GA. The vessel was exposed to the JONSWAP distributed waves with significant wave height 3.0 m and peak frequencyωi=0.65 rad/s, ς1=ς2=ς3=0.1. The population size of PSO is 50, the maximum velocity is set as 1, the iteration is 200, and the weight Wstart=0.9, Wend=0.4. The parameters are shown in Table 2; the simulation curve is shown in Figures 3 and 4.Table 2
Gains of observer using GA and PSO.
Matrix
GA
PSO
K
31
0.1
×
(
4
5
4
)
×
diag
(
1.265,1.543,1.881
)
0.1
×
(
4
5
4
)
×
diag
(
1.221,1.531,1.486
)
K
61
(
4
5
4
)
×
diag
(
1.278,1.197,1.472
)
(
4
5
4
)
×
diag
(
1.577,1.493,1.216
)
K
62
diag
(
1.631,1.727,1.735
)
diag
(
1.973,1.425,1.834
)
K
p
a
diag
(
1.498,1.864,1.399
)
×
γ
3
diag
(
1.288,1.606,1.339
)
×
γ
3
K
f
diag
(
1.887,1.635
)
diag
(
1.425,1.715
)Figure 3
The estimated position of ship by observer using GA and PSO.Figure 4
The real position of ship using GA and PSO.From Figure3, it can be seen that the observer using GA and PSO can estimate the low-frequency position of ship effectively. From Figure 4, we can see that the proposed PSO optimization observer can give better results, with a small position and heading bias.


### 6.2. Comparison of the MMAO and Single Observer

In order to demonstrate the good performance of the proposed MMAO, the simulation results are compared with the single observer. The JONSWAP wave spectrum with significant wave height is set between at 0 and 6 m. A set of four candidate values of the peak frequency were selected as{0.6,0.68,0.79,0.93} rad/s. The simulation curves are shown in Figures 5–9.Figure 5
The estimated position of ship using single observer.Figure 6
Evolution of dynamic weights of MMAO.Figure 7
The estimated position of ship using MMAO.Figure 8
The estimated velocity of ship using MMAO.Figure 9
Comparison of the real position of ship using single observer and MMAO.Figure6 shows the dynamic weights in the MMAO and the dynamic weights adaptively track the changes in the sea state. In this simulation, in the first 1500 seconds, the dominant wave frequency is 0.93 (rad/s). During this interval of time, p2, p3, p4 go to zero; in the next 1500 seconds the dominant wave frequency is 0.79 (rad/s). In the next 1500 seconds, the dominant wave frequency is 0.68 (rad/s) and in the last 1500 seconds, the dominant wave frequency is 0.60 (rad/s). Figures 7 and 8 show the time evolution of the low frequency estimation of the position and velocity of the ship. Figure 9 demonstrates that the ship using MMAO has better disturbance rejection property than single observer.


### 6.3. Comparison of MMAK and MMAO

In this subsection the proposed MMAO is compared with MMAK designed by [6]. The simulation results are shown in Figures 10 and 11.Figure 10
Evolution of dynamic weights of MMAK.Figure 11
Comparison of the real position of ship using MMAK and MMAO.Figure10 shows that the dynamic weights of the MMAK adaptively track the changes in the sea state. From Figure 5, we can draw the conclusion that the MMAK can switch to the best subobserver according to the varying seas. From Figure 10, it can be seen that the vessel using MMAO has a small position bias, and we can see that the MMAO has a better performance than MMAK.


## 6.1. Verify the PSO

In order to demonstrate the efficiency of the proposed PSO, optimization observer by PSO is compared with an observer using GA. The vessel was exposed to the JONSWAP distributed waves with significant wave height 3.0 m and peak frequencyωi=0.65 rad/s, ς1=ς2=ς3=0.1. The population size of PSO is 50, the maximum velocity is set as 1, the iteration is 200, and the weight Wstart=0.9, Wend=0.4. The parameters are shown in Table 2; the simulation curve is shown in Figures 3 and 4.Table 2
Gains of observer using GA and PSO.
Matrix
GA
PSO
K
31
0.1
×
(
4
5
4
)
×
diag
(
1.265,1.543,1.881
)
0.1
×
(
4
5
4
)
×
diag
(
1.221,1.531,1.486
)
K
61
(
4
5
4
)
×
diag
(
1.278,1.197,1.472
)
(
4
5
4
)
×
diag
(
1.577,1.493,1.216
)
K
62
diag
(
1.631,1.727,1.735
)
diag
(
1.973,1.425,1.834
)
K
p
a
diag
(
1.498,1.864,1.399
)
×
γ
3
diag
(
1.288,1.606,1.339
)
×
γ
3
K
f
diag
(
1.887,1.635
)
diag
(
1.425,1.715
)Figure 3
The estimated position of ship by observer using GA and PSO.Figure 4
The real position of ship using GA and PSO.From Figure3, it can be seen that the observer using GA and PSO can estimate the low-frequency position of ship effectively. From Figure 4, we can see that the proposed PSO optimization observer can give better results, with a small position and heading bias.


## 6.2. Comparison of the MMAO and Single Observer

In order to demonstrate the good performance of the proposed MMAO, the simulation results are compared with the single observer. The JONSWAP wave spectrum with significant wave height is set between at 0 and 6 m. A set of four candidate values of the peak frequency were selected as{0.6,0.68,0.79,0.93} rad/s. The simulation curves are shown in Figures 5–9.Figure 5
The estimated position of ship using single observer.Figure 6
Evolution of dynamic weights of MMAO.Figure 7
The estimated position of ship using MMAO.Figure 8
The estimated velocity of ship using MMAO.Figure 9
Comparison of the real position of ship using single observer and MMAO.Figure6 shows the dynamic weights in the MMAO and the dynamic weights adaptively track the changes in the sea state. In this simulation, in the first 1500 seconds, the dominant wave frequency is 0.93 (rad/s). During this interval of time, p2, p3, p4 go to zero; in the next 1500 seconds the dominant wave frequency is 0.79 (rad/s). In the next 1500 seconds, the dominant wave frequency is 0.68 (rad/s) and in the last 1500 seconds, the dominant wave frequency is 0.60 (rad/s). Figures 7 and 8 show the time evolution of the low frequency estimation of the position and velocity of the ship. Figure 9 demonstrates that the ship using MMAO has better disturbance rejection property than single observer.


## 6.3. Comparison of MMAK and MMAO

In this subsection the proposed MMAO is compared with MMAK designed by [6]. The simulation results are shown in Figures 10 and 11.Figure 10
Evolution of dynamic weights of MMAK.Figure 11
Comparison of the real position of ship using MMAK and MMAO.Figure10 shows that the dynamic weights of the MMAK adaptively track the changes in the sea state. From Figure 5, we can draw the conclusion that the MMAK can switch to the best subobserver according to the varying seas. From Figure 10, it can be seen that the vessel using MMAO has a small position bias, and we can see that the MMAO has a better performance than MMAK.


## 7. Conclusion

In order to improve the performance of observer under varying seas, a multiple model adaptive nonlinear observer is proposed for dynamic positioning ship. The MMAO consists a bank of nonlinear observer that relies on measurements of the vessel’s position, velocity, and acceleration. PSO is proposed to tune the gains of the observer. The computer simulation in this paper illustrates that the MMAO can switch to the best observer according to the sea state, and the MMAO has a better performance than single observer. The PSO can improve the performance of observer and the control efficiency of the DP.



---

*Source: 893081-2013-09-16.xml*