International Journal
on Marine Navigation
and Safety of Sea Transportation
Volume 5
Number 1
March 2011
111
1 INTRODUCTION
The Auto Track DP modes enable the vessel to fol-
low a predefined path, described by a set of way-
points, with a high degree of accuracy [Kongsberg,
2008]. These modes cover both low-speed and high-
speed operations using different control strategies.
The system can automatically switch between the
strategies depending on the requested speed.
Figure 1. DP MSV (multipurpose supply vessel) in Auto Track
mode. [Kongsberg, 2008]
The positions of the waypoints and the vessel's
heading and speed that are to be used for each track
section are specified by the operator and stored in
waypoint tables. Waypoints can be inserted, modi-
fied and deleted as required. The vessel's heading is
controlled by the following functions:
Present Heading
Set Heading
System Selected Heading
The speed of the vessel along each section of the
track can either be taken from the waypoint table or
specified on-line by the operator using the Set Ves-
sel Speed function.
Depending on the thrusters’ installation and the
vessel design, the maximum speed for a vessel in
Auto Track low speed mode should not exceed ap-
proximately 3 knots since the effect of the lateral
thrusters is reduced considerably at speeds higher
than this.
The figure 1 shows the track a vessel will follow
in Auto Track low speed mode according to the in-
formation contained in the table:
Table 1. Autotrack parameters
___________________________________________________
Waypoint North/East Co-ordinates Speed Heading
No. [m] [m/s] [°]
___________________________________________________
1 1501530 / 503600 0.3 0
2 1501570 / 503680 0.4 330
3 1501650 / 503630 0.4 300
4 1501740 / 503900 0.3 300
___________________________________________________
The operator can select between two alternative
strategies for passing waypoints:
Slowing down at each waypoint before continu-
ing to the next (used when the vessel must remain
on track, even during sharp turns)
Passing the waypoint at a constant speed on a
segment of a circle. The circle's radius can be cal-
culated automatically according to the vessel
speed, the angle of turn and the vessel's turning
characteristics.
Path Following Problem for a DP Ship
Simulation Model
P. Zalewski
Maritime University of Szczecin, Szczecin, Poland
ABSTRACT: For the dynamic positioning (DP) operations the equations of ship’s motion can be simplified
to a 3 degrees-of-freedom (DOF) model. DP station-keeping is a set-point regulation problem, i.e., forcing the
output of the ship to maintain a constant reference position and heading. However since offshore vessels re-
position themselves at different locations, the path following problem has to be solved. The article presents
various solutions of this problem including fuzzy control design which could be utilised in autonomous ship
simulation as well.
112
Both strategies generally require solving of the
path following problem to acquire the desired dy-
namical behaviour of the vessel depending on the re-
sultant thrusters’ set points (geometric and dynamic
problem).
DP vessels have different types of thrusters
tunnel thrusters, azimuth thrusters, main prop and
rudders are the most commend used. Both RPM and
pitch controlled. The pitch follow-up control is per-
formed from the thruster process station. Set points
are received cyclically from the DP controller via
the dual communication network or by use of ana-
logue signals. A PID-controller compares actual
pitch against set value and controls the hydraulic
pitch control valve accordingly.
It would then be convenient to parameterize all
path waypoints in terms of a continuous path and
constraining the vessel to this path [Skjetne, 2005].
The heading of the vessel could be taken as the di-
rection of the tangent vector along the path, or simp-
ly as a constant heading, usually pointed against the
environmental forces like waves and wind (system
selected heading), or set by the DP operator accord-
ing to operational demands (for instance to keep the
transducer in the range of HPR transponders).
The desired dynamic behaviour along this path
would be in first strategy zero speed (fixed position-
ing) at the waypoints, and when moving along the
path from one waypoint to the next the desired path
speed should be commanded online by the operator.
In the second strategy the desired dynamic behav-
iour along the path would be constant speed.
In the works done so far [Skjetne, 2005] the two
methods were used: 1) starting with an already
available tracking controller, and then converting
this into a manoeuvre regulation controller; 2) start-
ing with a parameterized path and a dynamic as-
signment along the path, designing the control, and
tying together the geometric and dynamic objectives
with the final pick of an update law. The second
method seems to be more flexible and has the ad-
vantage that the path variable can be a dynamic state
integrated online in the controller to satisfy the dy-
namic assignment.
These methods will be shortly presented in the
following sections and together with the foundations
of the fuzzy logic controller combining both.
2 CLASSIC TRACKING CONTROLLER IN DP
SYSTEMS
The main functions to be performed in order for a
dynamic positioning system to control a given vessel
position (x, y) and heading (ψ) are [Cadet, 2003]:
Estimate vessel motion
Measure vessel response
Determine error between prediction and meas-
urement
Determine corrective action to be applied
Calculate and allocate appropriate command to
thrusters to achieve desired corrective action
Figure 2 presents block diagram of a DP system.
The kernel of this system is the simplified hydrody-
namic vessel model. This model is a set of equations
of motion that is used to predict the motion of the
vessel when known forces and moment are applied.
In order to separate the wave induced oscillatory
part of the motion from the remaining part of the
motion, the total vessel motion is modelled as the
added outputs of a low-frequency model (LF-model)
and a high-frequency model (HF-model). The HF-
model represents oscillatory wave components in the
vessel motion. The LF-model represents motions in-
duced by wind, thrust and current in surge, sway and
yaw. The low frequency portion of the model is con-
trollable by means of thrusters. The algorithm calcu-
lates values of vessel’s state vector (position, head-
ing and motion variables) by measurements filtration
and then it changes resulting force demand - thrust-
ers allocation to meet position, heading and motion
settings [Zalewski, 2010].
Figure 2. Block diagram of DP system. [Kongsberg, 2008]
113
Thrusters Allocation block is usually ship specific
PID controller responsible for achieving dynamic
and geometric objectives formulated by the Carrot
Computation block. In a Dynamic Positioning appli-
cation a Kalman filter is used to estimate the state of
the vessel (for which a dynamics model has been
developed) based on noisy measurements from ref-
erence systems and sensors.
2.1 Kalman filtering
The estimation problem solved by the Kalman filter
can be expressed as follows [Cadet, 2003]: how to
optimally estimate the state of a vessel with an ap-
proximate knowledge of the vessel dynamics (im-
perfect mathematical model) and with noisy meas-
urements from sensors and position reference
systems? What is the best state estimate one can get
out of all that?
A discrete random dynamic system is described
by two equations in the Kalman filter [Zalewski,
2010]:
state equation (structural model of the process):
kkkk
wxAx +=
11
(1)
measurement equation (measurement model):
kkkk
vxHz +=
(2)
where: x n
th
dimension state vector,
w r
th
dimension state disturbance vector,
z m
th
dimension measurement vector,
vp
th
dimension state disturbance vector
(measurement noise),
A n×n dimension transition matrix,
H m×n dimension measurement matrix,
r ≤ n, p ≤ m.
Besides, it is assumed for disturbance vectors w
and v that they are Gaussian noise of normal distri-
bution, of zero mean vector and are mutually not
correlated. The state equation describes the trend of
the vector concerned, and the measurement model
gives the functional dependence of the measurement
on this vector. The solution to the equation system
(1), (2), taking into consideration the limitations im-
posed on disturbance vectors, can be reached by the
Kalman filter.
The estimation of the state vector in the filter can
be presented by the equations below:
state vector forecast:
+
=
1
ˆˆ
kkk
xAx
(3)
where:
n
k
x
ˆ
forecast or a-priori estimated val-
ue of the state vector at k-moment,
− a-
posteriori estimated value of the state vector,
covariance value of the forecast of state vector:
k
T
kkkk
QAPAP +=
+
1
(4)
where Q is the matrix of the covariance of the state
disturbance (of vector w), index T means matrix
transposing,
filter amplification matrix:
[ ]
1
+=
k
T
kkk
T
kkk
RHPHHPK
(5)
where R is the covariance matrix of measurement
disturbance (of vector v),
estimate of the state vector from filtration
+
k
x
ˆ
after
making measurement z
k
:
[ ]
+
+=
kkkkkk
xHzKxx
ˆˆˆ
(6)
covariance matrix of the estimated state vector:
[ ]
+
=
kkkk
PHKIP
(7)
where I is the identity matrix.
In DP systems controlling vessel’s position with
fixed heading in two dimensions generally the fol-
lowing values are to be estimated: position coordi-
nates (ϕ, λ), projections of the speed vector in rela-
tion to the bottom onto the meridian (N-S axis) and
the parallel (E-W axis) (V
N
, V
E
), acceleration vector
projections onto the meridian and the parallel (a
N
,
a
E
) and the projections of acceleration vector deriva-
tives in relation to the bottom onto the meridian and
the parallel (a’
N
, a’
E
). In this case the state vector
will have the following elements:
[ ]
T
ENENENk
aaaaVVx ',',,,,,,
λϕ
=
(8)
The measured values will be: position coordinates
of the positioning system used (ϕ
PS
, λ
PS
), speed com-
ponents in relation to the meridian and the parallel
from log or DR navigation performed by the posi-
tioning system used (position changes derivatives
V
N
, V
E
), acceleration components in relation to the
meridian and the parallel from an inertial transform-
er MRU (a
N
, a
E
). So the measurement vector will
have the following elements:
[ ]
T
ENENPSPSk
aaVVz ,,,,,
λϕ
=
(9)
Some assumptions like values of variance and
covariance for particular measurements have to be
done, for instance: DGPS system σ
ϕ
= 1.0 m; σ
λ
=
1.0 m; coordinates not correlated; speed components
σ
V
= 0.1 m/s; acceleration components σ
a
= 0.01
m/s
2
.
114
Figure 3. Kalman filter simulation of vessel x co-ordinate esti-
mate.
In the figure 3 the results of the Kalman filter
work referring to one co-ordinate are presented
(simulated in Matlab
®
). Blue are measurements, red
are unknown true values, green are estimated values.
Based on the simplified vessel model, and using
the previous position estimate of the vessel, the pre-
diction step of the Kalman filter gives us a predic-
tion of the vessel position. Based on the forces act-
ing on the vessel, on the vessel model and on the
previous position estimate, this is where the DP sys-
tem thinks the vessel is.
3 PARAMETERIZED PATH TRACKING
CONTROLLER IN DP SYSTEMS
As stated in the section 1, the DP Auto Track in-
volves two tasks called the geometric task and the
dynamic task. While the main concern is to satisfy
the geometric task, the dynamic task, further speci-
fied as a speed assignment, ensures that the system
output follows the path with the desired speed. The
main ideas were published in [Skjetne, 2005].
3.1 Parameterization of path from waypoints –
geometric assignment
One method to generate a path of class C
r
(C
r
means
that path function f is continuously differentiable up
to order r) is to first specify a set of n+1 waypoints,
and then construct a sufficiently differentiable curve
that goes through the waypoints by using splines and
interpolations techniques.
Designating the desired path as
( )
θ
d
p
we get n
curves
( )
θ
id
p
,
, i = 1, 2, … n, between the way-
points. Each of these is expressed as a polynomial in
θ of certain order. Then the expressions for the sub
paths are concatenated at the waypoints to assemble
the full path. To ensure that the overall path is suffi-
ciently differentiable at the way-points, the order of
the polynomials must be sufficiently high. We will
consider the path in two dimension space
.
2
Let
l = {1, 2, … n} be a set of indices identifying each
sub path. A column vector will be stated as col. The
overall desired curve is denoted
( ) ( ) ( )( )
[ ]
,,0,,col nyxp
ddd
=
θθθθ
while individual
segments as
( ) ( ) ( )
( )
,,,col
,,,
liyxp
ididid
=
θθθ
and
( )
}1{,,col += nliyxp
iii
are the waypoints.
Conveniently it is to assume that θ will reach an
integer value at each waypoint, starting with θ = 0 at
WP 1 and θ = i at WP i. The differentiability re-
quirement
( )
r
d
Cp
θ
, means that at the connection
of two sub paths, the following conditions must
hold:
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
θθθθ
θθθθ
θθθθ
θ
θ
θ
θ
θ
θ
θ
θ
θ
θ
θ
θ
θ
θ
θ
θ
θθθθ
rrrr
id
i
id
i
id
i
id
i
id
i
id
i
id
i
id
i
id
i
id
i
id
i
id
i
yyxx
yyxx
yyxx
,
1
1,
1
,
1
1,
1
,
1
1,
1
,
1
1,
1
,
1
1,
1
,
1
1,
1
limlimlimlim
limlimlimlim
limlimlimlim
==
==
==
(10)
where:
r
id
r
id
r
id
r
id
y
y
x
x
rr
θθ
θθ
=
=
,
,
,
,
,
for
}.1{\li
The task solution is to determine the coefficients
},{
,, ijij
ba
of k order polynomials:
( )
( )
ii
k
ikid
ii
k
ikid
bbby
aaax
,0,1,,
,0,1,,
+++=
+++=
θθθ
θθθ
(11)
According to equations (11) for each sub path
there are
)1(2 + k
unknowns ((k+1) unknowns a
i
and (k+1) unknowns b
i
) so there are
)1(2 +kn
un-
known coefficients in total to be determined for the
full path. Two methods are usually used for calculat-
ing these coefficients. After linearization of equa-
tions set (11) we will get
linear equations
cA =
φ
for the full path, which can be solved by a
matrix inversion operation:
.
1
cA
=
φ
However, as
the number n sub paths increases, this soon encoun-
ters numerical problems in the inversion of A. In-
stead, it is possible to calculate the coefficients for
each sub path independently. To ensure the desired
continuity at the connection points the numerical
values which are common for the neighboring sub
paths must be defined.
For a k’th order polynomial
( )
θ
id
x
,
we have that
( )
0
,
=
θ
θ
j
id
x
for
.1+ kj
Hence, we can form equa-
tions from the first k derivatives of
( )
θ
id
x
,
:
C
0
:
n22
equations for
li
:
0 10 20 30 40 50 60 70 80 90 100
8
8.5
9
9.5
10
10.5
11
11.5
12
12.5
Time [ s ]
x [m]
Kalman filter simulation
115
( ) ( )
( ) ( )
1,1,
,,
11
++
==
==
iidiid
iidiid
yiyxix
yiyxix
(12)
C
1
: four equations for the first and last n-waypoint:
( ) ( )
( ) ( )
nnndnnnd
dd
yynyxxnx
yyxxxx
==
==
++ 1,1,
121,121,
00
θθ
θθ
(13)
)1(22 n
equations for intermediate waypoints:
( )
( )
( )
( )
1,...,1
)(
)(
,...,2
)(1
)(1
2,
2,
11,
11,
=
=
=
=
=
=
+
+
+
+
ni
yyiy
xxix
ni
yyiy
xxix
iiid
iiid
iiid
iiid
λ
λ
λ
λ
θ
θ
θ
θ
(14)
where
0>
λ
is a design constant.
5,0=
λ
means
that the slope at WP i is the average of pointing
against WP i-1 and WP i+1; while
5,0<
λ
means
the slope is higher and
5,0>
λ
means the slope is
lower.
C
j
:
n22
equations for
li
, if derivatives of order
2j
are 0:
( ) ( )
( ) ( )
00
0101
,,
,,
==
==
iyix
iyix
jj
jj
idid
idid
θθ
θθ
(15)
which gives
nj 2)1(2 +
equations to solve for
nk 2)1( +
unknowns. The path generation problem
is now set up as n linear, decoupled sets of equations
iii
bA =
φ
;
li
; where the unknown vector
φ
i
is:
)1(2
)
0,...,,0,...,,
}{,}({col
+
==
=
k
R
kjijkjiji
ba
φ
(16)
and A
i
and b
i
are formed correspondingly according
to the above equations.
3.2 Dynamic assignment
The second task is to satisfy a desired dynamic be-
haviour along the path. This can be expressed in
terms of a time assignment, speed assignment, or ac-
celeration assignment along the path [Skjetne,
2005].
A time assignment means to be at specific points
along the path at specific time instants. For a contin-
uous parameterization
( )
θ
d
p
specific values like
21
,
θθ
; etc., must correspond to specific time instants
t
1
,t
2
; etc., dependent through a design function
( )
.
t
v
so that
)(
11
tv
t
=
θ
and
)(
22
tv
t
=
θ
.
A speed assignment is to obtain a desired speed
along the path. If
( )
θ
d
p
is a continuous parameteri-
zation this can be translated into a desired speed for
θ
. This desired speed may depend on the location
along the path given by
θ
, or it may explicitly de-
pend on time. A natural choice is therefore to ex-
press the desired speed for
θ
as a design function
),( tv
s
θ
=
.
An acceleration assignment is to obtain a desired
acceleration along the path. For a continuous param-
eterization
( )
θ
d
p
this can be expressed by a design
function
),,( tv
θθ
α
=
for
θ
, which may depend on
the speed
θ
along the path in addition to
θ
and t.
4 FUZZY CONTROL PATH TRACKING
The design of fuzzy controller solving geometric and
dynamic tasks presented in section 3 can be similar
to the one described in [Zalewski, 2009].
In case of system being a DP ship, steered in Au-
to Track mode in accordance to the designed path,
u(n) will be vector of thrusters setting, y(n) will be a
vector containing six variables defining actual mo-
tion in 3-degrees of freedom: P
xy
- actual position of
selected reference point stored as two variables, v
x
-
actual longitudinal (advance) velocity, v
y
- actual
transverse (lateral) velocity, ω - actual angular (rota-
tion) velocity, ψ actual ship’s heading; and S(n)
will be a vector containing six variables defining re-
quired motion in 3-dof: P
xyr
- required position of
reference point stored as two variables, v
xr
- required
advance velocity, v
yr
- required lateral velocity, ω
r
-
required rotation velocity, ψ
r
- required ship’s head-
ing At time n, y(n) and S(n) are used to compute the
input variables of the fuzzy controller (effect of
thrusters setting on motion): P
xy
- deviation be-
tween required and actual position, v
x
- difference
or deviation between required and actual longitudi-
nal (advance) velocity, v
y
- difference or deviation
between required and actual transverse (lateral) ve-
locity, ω - difference or deviation between required
and actual angular (rotation) velocity, ψ differ-
ence or deviation between required and actual ship’s
heading. So generally the input variables vector can
be designated by:
( ) ( ) ( )
nynSne =
(17)
Input variable scaling factors are used to conven-
iently manipulate the effective fuzzification on the
scaled universes of discourse. The scaled factors
used for e(n) vector in presented research are nor-
malization constants of the five mentioned devia-
tions. Assuming the scaling factors for deviations as
vector K
e
the scaled input vector is: