Chapter 3

Inertial Measurements

3.1 Introduction
Over the years a variety of solutions and control strategies has been proposed to
the DP control problem. A signicant industrial contribution was the application
of Kalman lters and optimal control (Balchen et al. 1976, Grimble et al. 1980,
Balchen et al. 1980). These model based strategies separated the rapid, purely
oscillatory (zero mean) motion caused by the rst order, linear (i.e. proportional to
the instantaneous wave height and with a frequency content equal to the incoming
waves) wave loads, from the more slowly varying forces due to nonlinear wave
eects, so called wave drift forces, wind and ocean currents. The applied thrust was
then calculated from the estimated low-frequency motion thus reducing thruster
modulation and wear and tear.
The objective of the DP system is therefore to counteract constant and the slowly-
varying disturbances due to:

• Higher order wave loads (wave drift).

• Ocean currents.

• Wind forces (not wind gusts).

All these three components have stationary contributions, and some form of in-
tegral action will be required. Furthermore, the vessel’s desired heading angle
should be selected such that the power required to reject the constant force com-
ponents is minimized. This can be achieved automatically by controlling the x-
and y-positions provided that the vessel’s reference point is located at a minimum
distance fore of the centre of gravity (Pinkster and Nienhuis 1986). A more versa-
tile and sophisticated concept called weather optimal positioning control (WOPC)
was introduced by Fossen and Strand (2001).
28 Inertial Measurements

It is not required to estimate the varying disturbances in order to stabilize the

system. Disturbance feed-forward will be an optional feature attenuating devia-
tions caused by varying disturbances. However, not all disturbance components
are measurable. Estimating the dynamic wind forces can be done eciently with a
wind sensor and the knowledge of the ship’s wind coecients, and commercial DP
systems include this “wind feed-forward” option. On the other hand, producing an
estimate of the wave drift forces is non-trivial, and practical measures for handling
these particular dynamic forces have not been explicitly taken until very recently.
Even though PID-like control eciently rejects constant disturbances, it is not
the best option for attenuating varying disturbances in a second order mechanical
system. This stems from the simple fact that a PID-controller, linear or nonlin-
ear, cannot cope with the varying force disturbance due to the phase lag involved.
After all, the control is calculated as a sum of the system’s positions, integrated
position errors and velocities. The resulting behavior will therefore be oscillatory
in contrast to the case where the varying disturbance is completely cancelled out
by an oppositely directed thrust force. Introducing more damping could theoret-
ically reduce the amplitudes of the resulting motion, but due to wave-frequency
residues in the low-frequency velocity estimate, increased damping simultaneously
contributes to thruster modulation. There is therefore a practical limit for the
derivative gain’s magnitude.
In Aalbers et al. (2001) the authors proposed two “wave feed-forward” techniques
to suppress the nonlinear wave eects. They illustrated the performance experi-
mentally with a model ship. The wave drift forces were estimated by measuring
the relative water motion around the waterline of the hull using ten probes. The
estimated drift forces were, together with estimated wind forces, used directly in
the PID feedback loop. For the same thrust power, the position deviations de-
creased. The method is, however, impractical because of the required number of
wave probes and their robustness and expected length of life. It seems that such
an installation could turn out to be expensive and somewhat unreliable.
An alternative method applicable to ships, or in fact to all mechanical systems, is
to measure the accelerations. Active use of measured acceleration for control can
be regarded as manipulating the system’s mass. Negative acceleration feedback
(AFB) increases the mass and positive feedback decreases it. By negative feedback
the system is made virtually heavier as seen from the disturbances, thus their
eects are attenuated and the positioning performance increases. In the opinion of
the author, it seems more attractive to employ measured acceleration in positioning
control of surface vessels to further suppress varying disturbances in general and
slowly varying wave drift in particular rather than implementing a wave feed-
forward approach (Aalbers et al. 2001) because:

• The sensor and engineering cost will be much lower. High-precision commer-
cial inertial measurement units (IMUs) are becoming increasingly aordable
and can be easily interfaced with existing control systems.

• The sensor equipment is intended for aerospace applications and is therefore

extremely solid, reliable and robust.
3.1 Introduction 29

• All kinds of slowly varying disturbances, not only wave drift, will be atten-
uated.

• Acceleration feedback is applicable to all kinds of vessels, in contrast to the

rst force estimation method proposed by Aalbers et al. (2001) which is
primarily applicable to craft with large surface piercing structures such as
tankers.

Active use of measured acceleration is well known in aerospace applications (Blakelock

1991) but is rarely applied in other areas like robotics and ship control simply be-
cause it is in general superuous in stabilizing such systems. If, on the other hand,
the system parameters are uncertain, the tracking performance and robustness may
suer. This is particularly evident in feedback linearization designs. Employing
measured acceleration relaxes the need of an accurate model description because
the right-hand side of Newton’s second law
X
Fk = ma (3.1)
k

divided by the mass (and thereby the model itself) is actually being measured.
Consequently, it seems plausible to expect improved robustness, tracking perfor-
mance and disturbance rejection when AFB is constructively applied.
Using the acceleration in a feedback loop has been referred to as a direct approach
(de Jager 1994) as opposed to indirect where the acceleration signal is used in an
observer to improve the state estimates. de Jager (1994) recommended to consider
using measured acceleration either directly or indirectly, but not in combination,
and he reported increased tracking performance by counteracting uncertainty in
the inertia matrix of a 2D Cartesian manipulator.
For rigid robots, Luo and Saridis (1985) suggested using a decentralized (diagonal)
static linear controller for rigid robots assuming that joint positions, velocities and
accelerations were available. Stability in the sense of Lyapunov and performance
of this controller was later studied by Studenny and Bélanger (1984). In Kosuge
et al. (1989) the authors suggest using low-pass ltered acceleration to improve
the performance and robustness of a feedback linearization design of a two-link
planar robot. The “disturbance” due to inaccurate model parameters introduced
in the linearization was reduced in addition to the environmental disturbances.
Complete disturbance cancellation can only be achieved as the acceleration gain
tends towards innity, and due to unmodeled dynamics, time delays, imperfect
measurements, and other practical limitations, there will be an upper bound on
the acceleration gain.
In low speed ship control, such as DP, a combined approach after de Jager’s de-
nition is almost inevitable due to the heavy notch ltering the observer has to
perform. First order wave loads dominate the measured acceleration signal and
direct application leads to thruster modulation and unbearable wear on the equip-
ment. We let the notch ltered acceleration signal update both the state estimates
and the control law. This particular combined approach proved to be successful.
30 Inertial Measurements

3.1.1 Motivation: AFB in Mass-Damper Systems

Neglecting couplings and considering one degree of freedom (DOF) at a time, a

ship at low speed can be regarded as a linear mass-damper system

mẍ = dx + u + w (3.2)

where x and v = x are the states describing position and velocity,. u is the control
and w is a disturbance. Isolating the velocity equation and dening Tm = m/d we
get the transfer function

v 1 1
h(s) = (s) = (3.3)
w d (1 + Tm s)

The idea is now to incorporate an acceleration term in the control u, that is let u
be the sum
u = ha (s)v + uP ID (3.4)

where uP ID is to be designed later and ha (s) is some dynamic system. The eect
of the negative acceleration term ha (s)v is increased mass, the system’s virtual
mass becomes ma (s) = m + ha (s).
For attenuation of low-frequency disturbances, a suitable ha (s) could be a low-pass
lter with a gain specied as a fraction of the original mass m
m
ha (s) = (3.5)
1 + Tf s

For positive ’s the mass increases to (1 + ) m for frequencies below  f = 1/Tf
thus making the system virtually heavier and less inuenced by varying distur-
bances. For high frequencies the mass is left unaltered because |ha (j)|  0 for
 À f . The resulting transfer function for the velocity dynamics becomes

v v 1 1 + Tf s
(s) = (s) = (3.6)
uP ID w d (1 + Tm s) (1 + Tf s) + Tm s
a e 1 + Tf s
 (3.7)
d (s +  a ) (s +  e )

where
1 1 1+
a = , e = (3.8)
Tm 1 +  Tf

The asymptotic Bode plots from force to velocity (Figure 3.1) with and without
acceleration feedback show that, because the mechanical time constant has been
increased by feedback, the magnitude decreases for frequencies a <  < e .
The “power”, as given by the L2 -induced norm, of the velocity v has been reduced
for disturbances w in the frequency range a <  < e . For frequencies lower
than a , no reduction can be expected.
3.2 Inertial Measurements 31

TF from disturbance to velocity, |v/w(j ω)|

0
Gain [dB]

ω 1/T
m
1/T
f
ω
a e

Virtual mass |m (jω)|

a

m(1+α)
mass [kg]

ω 1/T
m
1/T
f
ω
a e

Frequency ω [rad/s]

Figure 3.1: Top: Bode plot of wv (s) with (solid) and without (dashed) acceleration
feedback. Bottom: The magnitude of virtual mass ma (s) = m + ha (s).

3.2 Inertial Measurements

3.2.1 Angular Rates

When using high-precision gyro measurements it is required to take the Earth’s

rotation into account, the n-frame can no longer be regarded as the inertial frame.
It is customary to express the rotations with respect to the ECI (Earth-centered
inertial) i-frame which by denition does not rotate. The ECEF (Earth-centered
Earth-xed) frame is also located in the center of the Earth but rotates along with
it.
An assumed error free gyro measures the rotation of the b-frame relative to the
inertial ECI frame i.
$bimu = $bib = $bie + $ ben + $ bnb (3.9)

The second component can be ignored for marine applications, while the rst can
be ignored if it is below the gyro noise level which means that for high-precision
gyros $bie must be taken into consideration.
If the n-frame is xed relative to the Earth, or as mentioned slowly-varying as in
32 Inertial Measurements

marine operations, the gyro measurement is reduced to

$ bimu = $ bie + $ bnb (3.10)

3.2.2 Linear Accelerations

Consider now the n-frame as the inertial frame. Suppose the IMU is located at
the distance pbbI from the origin of the body-xed frame. In the inertial frame,
the position of the IMU is then

pnnI = pnnb + Rnb pbbI (3.11)

and its velocity vnI

n
= p nnI is
¡ b ¢
vnI
n
= p nnb + Rnb S($bnb )pbbI = Rnb vnb + S($ bnb )pbbI (3.12)

 n = Rn S( b ). The acceleration an = p̈n is thus

where we used that R b b nb nI nI
³ ´
annI = p̈nnb + Rnb S2 ($ bnb ) + S($ bnb ) pbbI (3.13)
¡ ¢ ³ ´
= Rnb S($bnb )vnb
b
+ v nb
b
+ Rnb S2 ($ bnb ) + S($
 bnb ) pbbI (3.14)

To avoid handling angular accelerations $  bnb , let the accelerometer be located in

the origin of the b-frame, i.e. pbI = 0. Since the acceleration measurements are
b

decomposed in the b-frame, the assumed error free acceleration measurement fimu b

consisting of the actual acceleration anb given by

n

abnb = Rbn annb = v nb

b
+ S($ bnb )vnb
b
(3.15)

and gravity is written

fimu
b
= abnb  Rbn gn (3.16)

Here gn = [0, 0, g]T is the contribution from gravity. Notice the Coriolis eect
S($ bnb )vnb
b
caused by the rotation of the b-frame relative to the inertial n-frame.
Observe that the gravity term could have been expressed using the vessel-parallell
p-frame such that on component form
 
 sin 
gb = Rbn gn = Rbp gp =  sin  cos   g (3.17)
cos  cos 

because gp = gn .
In terms of the -vector, the error-free measured acceleration can be written as

fimu
b
=  1 + S( 2 ) 1  gb =  1 +  2 ×  1  gb (3.18)
3.3 Compensators 33

3.3 Compensators
The error-free measured acceleration (3.16) demonstrates that yacc b
6=  1 = v nb
b
.
In order to construct a “measured”  1 some kind of compensator has to be imple-
mented:
y 1 = fimu
b
 S($̂ bnb )v̂nb
b
+ R(ˆ
pb )gp   1 (3.19)
Even small roll and pitch angles will lead to gravity components in the acceleration
measurements along the surge and sway axes. Those components will, because
gravity is the dominating force acting upon the vessel, dominate. Consequently,
accurate measurements of surge and sway acceleration require good roll and pitch
measurements. An integrated navigation system estimates the position pnnb , the
orientation  nb and the velocities vnb
b
and $ bnb . The compensator (3.19) can thus
be realized using such systems.
An integrated navigation system is, however, not strictly required for gravity com-
pensation. Below two dierent types of g-compensation are proposed and dis-
cussed, one static and one dynamic. The static compensator uses only accelera-
tion measurements to remove the gravity forces, while the dynamic approach is
an attitude observer also utilizing gyro measurements. However, neither the static
nor the dynamic g-compensator is able to cancel the Coriolis component in

abnb =  1 +  2 ×  1 (3.20)

which means that unless the navigation system is able to accurately estimate
 1 = vnb
b
, this Coriolis eect cannot be removed and therefore isolating  1 is
generally speaking impossible and care must be taken when using the gravity
compensated acceleration measurement in the control design.
Still, in positioning operations at sea, the need for such integrated systems can be
relaxed as explained below.

3.3.1 Static Low-Speed Gravity Compensator

If we can assume that the velocities vnb

b
and $ bnb are small, the error-free measured
acceleration (3.16) can be approximated by

fimu
b
 v nb
b
 Rbp gp (3.21)

If the vessel is at complete rest

fimu
b
= Rbp gp (3.22)

from which can solve for the roll and pitch angles as follows
μ ¶
fy
 = arctan , fz > 0 (3.23)
fz
 
f
 = arctan  q
x
(3.24)
fz2 + fy2
34 Inertial Measurements

where we used the individual components of the measured signal fimu

b
= [fx , fy , fz ]T .
Provided the inclinations are small
fy fx
  (3.25)
g g
and consequently, the roll and pitch errors  and  exhibit a similar dependence
on the measurement errors fx and fy , that is
fy fx
     (3.26)
g g
A 1 mg accelerometer error thus gives a static roll and pitch accuracy of 1 mrad
 0.06 deg.
A successful g-compensation relies on good roll and pitch estimates (measure-
ments). Suppose that the velocities are small in the sense that gravity forces dom-
inate S($ bnb )vnb
b
such that (3.23)-(3.24) are reasonable approximations for roll and
pitch. The accelerometers are contaminated with an error acc = [fx , fy , fz ]T
due to various eects. The error-free measured acceleration fimu b
and the actual
measured accelerations are

fimu
b
= ab  Rbp gp (3.27)
yacc
b
= fimu
b
+ acc (3.28)

From the latter we can estimate ab using the estimated roll and pitch angles
ˆpb = [,
 ˆ ˆ, 0]T
âb = yacc
b
+ R̂bp gp (3.29)
The error ãb = ab  âb is thus
³ ´
ãb = fimu
b
+ Rbp gp  fimu
b
+ acc + R̂bp gp
³ ³ ´´
 ˆpb ) gp  acc
I  S(pb )  I  S(
˜pb  acc
= S(gp ) (3.30)

The linear approximations of the rotation matrices are valid for small inclinations.
Taylor expansions of the roll and pitch components of the attitude error  ˜pb =
ˆpb are
pb  

˜  fy

g
˜ fx
  
g
Consequently, by inserting these errors into (3.30), we see that a g-compensator
not only cancels out the gravity components but also removes the accelerometer
error acc on the x and y-axis
£ ¤T
ãb   0 0 fz (3.31)
3.3 Compensators 35

3.3.2 Dynamic Gravity Compensation

In order to measure dynamic roll and pitch accurately, a lter that integrates gyro
and accelerometer measurements must be designed. This is often referred to as
a vertical reference unit (VRU) and such can be implemented using a Kalman
lter or an observer. Here the observer developed in Vik and Fossen (2001) is
summarized and used as basis in a dynamic gravity compensation (DGC) scheme.
The “low-speed assumption” above that was required for using (3.23)-(3.24) as
roll-pitch measurements is no longer applicable for vessels exposed to incoming
waves. However, this requirement can be relaxed when employing angular velocity
measurements.
The observer is written (Vik and Fossen 2001):
1 h i 1
q̂ = Tq̂ (q̂)R(q̃) $ bimu + b̂gyro + K1 %˜q sgn(˜
 q )  (q̂)$ nin (3.32)
2 2
b̂ 1 1
gyro = Tgyro b̂gyro + K2 %̃q sgn(˜ q ) (3.33)
2
where q̂  H, is the four element unit quaternion
n o
H= q | qT q =1, q = [q , %Tq ]T , q  R, %q  R3 (3.34)

$ imu  R3 is the angular velocity vector measured by the gyros, b̂gyro  R3 is the
gyro bias, and $ nin = $ nie + $ nen where $nie is the earth rate vector and $ nen is the
angular velocity due to the movement of the ship over the Earth. Computation of
$ nie requires knowledge of true north. Tgyro  R3×3 is a diagonal time constant
matrix, and K1  R3×3 and K2  R3×3 are diagonal matrices. Finally,
· ¸
%̂Tq
Tq̂ (q̂) = (3.35)
ˆq I + S(%̂q )
· ¸
%̂Tq
(q̂) = (3.36)
ˆq I  S(%̂q )
%̃q and ˜q are components of the quaternion error q̃, which is computed from the es-
timated quaternion and a measurement quaternion derived from the accelerometer
based attitude measurements (3.23)-(3.24). Since the accelerometer attitude mea-
surement only needs to prevent the integrated gyro signal from drifting, the gains
are usually chosen very small. The low gain means that horizontal accelerations
have little inuence on the roll and pitch measurements. Thus, these measurements
can be used to compensate for the g-vector in the surge and sway acceleration mea-
surements. Moreover, since the accelerometers are used for attitude computation,
accelerometer bias will not inuence the surge and sway measurements.
DGC is then carried out after (3.29), that is
âb = yacc
b
+ R̂bp gp

where R̂bp now reects it is being calculated based on the VRU’s attitude estimates
q̂.
36 Inertial Measurements

Figure 3.2 shows the measured and g-compensated acceleration in a static test
using the Litton LN-200.
Surge acceleration [m/s ]
2

0.04

0.02

-0.02

-0.04
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
Sway acceleration [m/s ]
2

0.3

0.2

0.1

-0.1
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
Heave acceleration [m/s ]
2

-5

-10
500 1000 1500 2000 2500 3000 3500 4000 4500 5000
Time [s]

Figure 3.2: Measured and g-compensated accelerations. The noise level is about
120 μg (1), and the oset after compensation is a few μg on all three axes.

3.4 Positioning Control

In ship positioning, oscillatory motion caused by rst order wave loads domi-
nate the vessel’s dynamic behavior. By separating the wave motion from the
low-frequency dynamics, the body-xed velocities and linear accelerations can be
written as

 bnb =  bLF +  bW F (3.37)

$ bnb = $ bLF + $ bW F (3.38)

where the subscripts now identify the low-frequency and and wave frequency parts
instead of relative motion. The error-free acceleration measurement (3.16) can be
written as
b ¡ b ¢
fimu
b
= v LF
b
+ v W
b
F + S($ LF + $ W F ) vLF + vW F  Rp g
b b b p

= v LF
b
+ S($bLF )vLF
b
 Rbp gp + zbW F (3.39)
3.5 Conclusions 37

where zbW F is a signal dominated by the linear wave induced motion

¡ b ¢
zW F = S($bLF )vW
b b
F + S($ W F ) vLF + vW F
b
(3.40)

A reasonable assumption motivated by the small low-frequency velocities occurring

in a positioning operation is:

A1 The Coriolis terms from the low-frequency motion can be neglected, that is

S($ bLF )vLF

b
0 (3.41)

Hence, the measured acceleration can be approximated by

yacc
b
 v LF
b
 Rbp gp + zbW F + acc (3.42)

which is the sum of the low-frequency acceleration, gravity contributions and some
signal zbW F of a signicantly higher frequency content than v LF
b
. Therefore, using
a dynamic gravity compensator to remove g and acc together with a wave
b

lter (or observer) to remove the wave frequency components zbW F as depicted
in Figure 3.3, it is at least in a positioning operation possible to isolate v LF
b
=
 1 . Consequently, using  1 in a low-speed control design is possible also without

Figure 3.3: Realization of a g-compensation system and observer for reconstruction

of the low-frequency body-xed acceleration v LF
b
and angular velocity $ bnb .

utilizing an integrated navigation system.

3.5 Conclusions
The active use of negative acceleration feedback has been briey introduced by
regarding it as a change of the system’s mass. In order to utilize measured lin-
ear accelerations on surface vessels, some kind of dynamic compensation or pre-
processing has to be performed in order to relate the measurement to the b-frame.
By separating the low freqency motion from the wave induced motion (of relative
high freqency content), it was demonstrated that a particular VRU with gravity
compensation serves this purpose well as long as the vessel itself operates at low
speed.
38 Inertial Measurements
Chapter 4

Observer Design

4.1 Introduction
An observer lters available measurements to provide online estimates of the mea-
sured and unmeasured states within a system. In ship control, the most commonly
needed states are the low-frequency (LF) parts of the positions, the heading, the
velocities and stationary (or slowly varying) disturbances due to wind, ocean cur-
rent and nonlinear wave eects. Based on those estimated states, the controller
calculates its thrust demand. The three main objectives for the observer are:

• State estimation: To produce the state estimates from which the controller
calculates its desired propeller thrust forces and moments.

• Wave ltering: The estimator should attenuate the fast, oscillatory motion
due to rst-order wave loads. It is useless trying to compensate this sinu-
soidal (zero mean) behavior, and doing so will only lead to excessive thruster
system wear.

• Handling dead-reckoning: In the case of temporary sensor failure, the ob-

server must for a period of time be able to adequately predict the motion of
the ship such that the positioning operation can continue. This is a system
redundancy requirement which is mandatory with the classication societies
(Det Norske Veritas 1990).

The third objective implies that some kind of model based lter must be imple-
mented.
Classic solutions to the DP problem of surface vessels are output-feedback de-
signs using a state-estimator to lter out 1st-order wave induced motion from the
LF positions while reconstructing LF velocities (Balchen et al. 1976, Balchen et
al. 1980, Sælid et al. 1983, Grimble et al. 1980, Fung and Grimble 1983, Sørensen et
al. 1996). All these were realized using linear stochastic theory (Kalman Filters),
40 Observer Design

but also H -solutions been proposed (Katebi et al. 1997). Unfortunately, the lin-
earization of the nonlinear kinematics implies that the results are only valid locally.
However, if the nonlinearities satisfy a global Lipschitz-condition, a modication
(Reif et al. 1999) of the extended Kalman lter ensures global exponential stability.
Another approach with comparable performance is to utilize the model structure
and let the observer “linearize” itself about the measured compass heading. As
opposed to traditional extended Kalman-lters, the on-line explicit linearization
is avoided, and global stability properties are more easily established since the
nonlinear kinematics can be treated as a known time-varying block. Examples
are the passivation designs (Fossen and Strand 1999, Strand and Fossen 1999),
further extensions to higher order monotonic damping terms (Aamo et al. 2001),
and non-dissipative linear damping terms (Lindegaard and Fossen 2001b).
As discussed previously, using the measured accelerations, we are able to better
keep up with unmodeled disturbances like slowly varying wave forces which must
be counteracted by the control system. Slowly varying wave induced forces is a
phenomenon well known from nonlinear hydrodynamic theory (Faltinsen 1990),
yet they are dicult to express in a form suited for controller design. However,
feeding the measured accelerations uncritically into the closed loop system is not
recommended due to the high-frequency, large amplitude oscillations caused by
1st-order wave loads. Therefore, some kind of notch ltering of the measured
accelerations is required in order to remove the wave frequency components.
In this chapter we focus on extending the proven observer structure from Fossen
and Strand (1999). More specically, the contributions are:

• Identication of structural conditions unifying linear and nonlinear observer

design for surface vessel at low speed: If some structural constraints are satis-
ed, the nonlinear kinematics can be disregarded in the stability analysis and
linear design tools may be applied. While previous nonlinear observer de-
signs either assumed a passive, and thus stable, vessel-xed dynamics (Fossen
and Grøvlen 1998, Fossen and Strand 1999, Strand and Fossen 1999, Loría et
al. 2000, Aamo et al. 2001) or prescribed using non-diagonal observer gain
matrices to handle the possibly unstable sway-yaw dynamics (Robertsson
and Johansson 1998), we show that one can obtain uniform global exponen-
tial stability of the observer errors for all kinds of vessels using xed, diagonal
observer gains.
• Optional inclusion of velocity and acceleration measurements.

Common for all previously mentioned designs is that they are derived under the
assumption that only the positions and compass heading were available for feed-
back. Today high performance inertial measurement units (IMU) are becoming
increasingly aordable, and integrated navigation systems (INS) integrating IMU
and GPS reproduce not only positions but also velocities and linear accelerations
with great accuracy to a reasonable price. This development in sensor technology
is reected in the proposed designs.
Two dierent observers will be analyzed:
4.2 Common Model Description 41

1. The rst observer (Lindegaard and Fossen 2001a, Lindegaard and Fossen
2001b) contains one single wave model generating accelerations, velocities
and positions. It requires position measurements, and the structure allows
the inclusion of velocity and acceleration. Without acceleration measure-
ments the resulting error dynamics is shown to be UGES provided that the
wave model and selected observer gains satisfy the structural properties.
With accelerations, a specic bound on the yaw rate  must be imposed and
uniform semi-global exponential stability (USGES) can be guaranteed.

2. In the second observer (Lindegaard et al. 2002) the wave models for ac-
celeration, velocity and position and treated as separate phenomena. This
facilitates the tuning procedure signicantly, and a particular tuning proce-
dure based on pole placement is proposed. A similar structural constraint
as for the rst observer must also be assumed in this case, and USGES of
the error dynamics is shown.

4.2 Common Model Description

The two observers are quite similar, the main dierence between the two being in
the implementation of the wave model used to lter out rst order wave motion.
This section discusses the LF model description used in both designs.
We consider the dynamics of a vessel in three degrees of freedom, the horizontal
plane, and we choose to express the model in the Earth- and body-xed coordinate
frames. The body-xed frame coincides with the principal axes of the vessel and it
is rotated an angle y with respect to the Earth-xed frame. This transformation
of coordinates is represented by the orthogonal rotation matrix
 
cos y  sin y 0
R( y ) =  sin y cos y 0  (4.1)
0 0 1

4 d

and its time-derivative is R() =dt (R()) = SR()
 where the skew-symmetric
matrix S = ST is given by
 
0 1 0
S= 1 0 0  (4.2)
0 0 0

Let  = [x, y, ]T be the LF position vector where x and y are the North and East
positions respectively, and being the LF heading.  = [u, v, r]T contains the LF
body-xed velocities, i.e. surge, sway, and yaw. The LF ship model is assumed to
satisfy:

A1 The orientation angle between the Earth-xed and body-xed frame is the
42 Observer Design

measured heading y such that:

 = R( y )
b = T1
b b + Eb wb (4.3)
M = GRT ( y )  D +  + RT ( y )b

Here   R3 is the applied thruster force, M = MT > 0 is the sum of rigid

body mass and hydrodynamic added mass, D  R3×3 contains linear damp-
ing coecients and G  R3×3 describes the mooring forces. The bias forces
b  R3 are modelled as Markov processes with a positive semi-denite diag-
onal matrix Tb  R3×3 of time constants. wb  R3 is a bounded disturbance
signal, and Eb  R3×3 is a gain factor.

In the following we will frequently utilize a commutation property between the

Earth-xed parameters and the rotation R().

Property 4.1 A matrix A  R3×3 is said to commute with the rotation R() if

AR() = R()A (4.4)

Examples of matrices A satisfying Property 4.1 are linear combinations A =

a1 R() + a2 I + a3 kT k for scalars ai ,  and k = [0, 0, 1]T , the axis of rotation. Also
note that since R() is orthogonal, that is RT () = R1 (), Property 4.1 implies
that
A = RT ()AR() = R()ART () (4.5)
Furthermore, if A is nonsingular, A1 commutes with R() too. That is
A is nonsingular
AR() = R()A  A1 R() = R()A1 (4.6)

4.3 Observer With Consistent Wave Model

In Section 4.3.1 below we describe the vessel model and some of its properties with
attention given to how the individual measurements t into the model framework.
Section 4.3.2 concerns the lter design and stability analysis. An important part of
this section is the discussion regarding structural properties: The structure of gain
matrices updating the Earth-xed error dynamics can not be selected arbitrarily.
Section 4.3.3 presents some conclusions and remarks.

4.3.1 Complete Ship and Environment Models

The LF model is described by (4.3). Now, let xw  R3nw describe the rst order
wave-induced motion where nw denotes the number of states used to describe
the wave frequency motion in each degree of freedom (DOF). Here we let the
4.3 Observer With Consistent Wave Model 43

WF model be expressed entirely in the Earth-xed frame, and in accordance with

Fossen and Strand (1999), we employ a linear wave model on the form

x w = Aw xw + Ew ww (4.7)

where Aw is assumed Hurwitz, Ew  R3nw ×3 is a gain matrix, and ww  R3 is

a zero-mean, bounded disturbance input. In a Kalman lter setting ww should
be a white, Gaussian process. We are ready to impose two additional model
assumptions:

A2a For D = {dij } i, j = 1, · · · , 3, the elements d11 , d22 > 0.

A3a The bias time constant matrix Tb and each 3 × 3 sub-block of Aw satises
Property 4.1.

Note in Assumption A2a that there are no restrictions on neither d23 , d32 nor d33 .
A2a is thus less restrictive than assuming D + DT > 0 (Fossen and Strand 1999)
and its interpretation is that separate surge and sway motions are dissipative.
It does, however, include cases with potentially unstable sway/yaw dynamics
(Robertsson and Johansson 1998). The last assumption, A3a, implies that the
mean wave motion period, relative damping, and bias time constants in the North
and East directions are identical. It should be emphasized that this is not as re-
strictive as it may sound since the dominating frequency of the rst order wave
induced motions will be approximately the same in surge and sway.

First Order Wave Motion

A suitable linear representation of the oscillatory motion caused by the 1st-order

wave loads can be approximated by a set of three de-coupled linear transfer func-
tions
kwi s
h(s) = 2 (4.8)
(s + 2 i 0i s + 20i )
2

in each of the 3 DOF. The motivation for using a function of fourth order instead
of other approximations, is that this choice ensures that the transfer functions
between the excitation and positions, velocities as well as accelerations, will be
strictly proper.
A minimal realization with xw  R12 can in state-space be described by

x w = Aw xw + Ew ww (4.9)
   
0 I 0 0 0
  0 I 0
Aw =  , Ew = 
0 0 0 I  0 
0 0   diag (kw )
where  and  are diagonal matrices holding the the wave motion resonance
frequencies  oi and relative damping factors i , i = 1, · · · , 3 for the North, East
44 Observer Design

and heading respectively like this

 = diag(201 , 202 ,  203 ) (4.10)

 = diag(2 1  01 , 2 2  02 , 2 3  03 ) (4.11)

and kw = [kw1 , kw2 , kw3 ]T is a gain vector. Assumption A3 requires 01 = 02
and 1 = 2 . The Earth-xed wave induced position, velocity and acceleration
can be extracted from xw as follows

pw = Cp xw , vw = Cv xw , aw = Ca xw (4.12)

where Cp , Cv , Ca  R3×12
£ ¤
Cp = I 0 0 0 (4.13)
£ ¤
Cv = 0 I 0 0 (4.14)
£ ¤
Ca =   0 I (4.15)
¡ T ¢ ¡ T¢
d
Since dt R vw = dt d
R vw + RT v w =   y RT Svw + RT aw , the experienced
wave induced velocities and accelerations will in the body frame be given as

vw
b
= RT ( = RT ( y )Cv xw
y )vw (4.16)
³ ´
abw = RT ( y ) Ca   y SCv xw (4.17)

The acceleration term depending on measured rotation rate ry =  y can be re-

garded as a Coriolis-like term.
Collect the Earth-xed states in x1  R6+3nw and the body-xed in x2  R3
£ T ¤T
x1 = xw  T bT (4.18)
x2 =  (4.19)

and dene the block diagonal transformation matrix

¡ ¢ ¡ ¢
T( y ) = Diag(RT y , · · · , RT y , I) (4.20)

On compact form using Assumption A3 we get with x = [xT1 , xT2 ]T and w =

T
[ww , wbT ]T
x = TT ( y )AT( y )x + B + Ew (4.21)
where  
Aw 0 0 0
0 0 0 I
A=  (4.22)
0 0 T1
b 0
0 M1 G M1 M1 D
   
0 Ew 0
0 0 0
B= , E= (4.23)
0  0 Eb 
M1 0 0
4.3 Observer With Consistent Wave Model 45

Measurements

We intend to cover all combinations of position, velocity and acceleration measure-

ments. The positions are usually given in an Earth-xed reference frame, while
velocities and accelerations are given in a body-xed coordinate system.
There might be cases where not all kinds of measurements are available, either due
to sensor failure, or simply because that particular vessel was not equipped with
that kind of instrument. Denote the measurements y. We have that y  Rny where
3  ny  8 depending on the conguration. Dene 2 and 3 as the projections
extracting the measured velocities and accelerations respectively from the actual
three DOF velocity and accelerations vectors. It is possible to measure all three
velocities pretty accurately, which means that quite often 2 = I. While linear
accelerations are easy to measure, the angular acceleration is not. Therefore, most
likely only the accelerations in surge and sway are available and hence ny3 = 2
and · ¸
1 0 0
3 = (4.24)
0 1 0

Let y1  R3 contain the Earth-xed positions and compass heading, y2  Rny2

the vessel-xed velocities and y3  Rny3 accelerations. Then,
y1 =  +  w =  + Cp xw (4.25)
¡ ¢
y2 = 2  + R Cv xw
T
(4.26)
¡ ¢
y3 = 3  + awb
³ ´ ³ ´
= 3 M1 GRT  + RT b  D + 3 RT Ca   y SCv xw (4.27)

Compactly written
y = Cy ( y,
 y )x + Dy  (4.28)
where

Cp I
Cy ( y,
 y) =  2 Cv RT ( y ) 0
1
3 (Ca   y SCv )RT ( y ) 3 M GR (
T
y)

0 0
0 2 
3 M1 RT ( y ) 3 M 1
D
£ ¤T
Dy = 0 0 MT T3 (4.30)
When only positions are available, ny2 = ny3 = 0, the model is reduced to the
traditional DP observer problem.
In the stability analysis below it will be convenient to make an assumption on how
the velocity and acceleration feedback is congured:

A4a Let i = Ti i  R3×3 , i = 2, 3. Valid congurations are those which allow
i to commute with R(), that is R()i = i R() for all   R.
46 Observer Design

Objective

For the model (4.21) with output (4.28)-(4.30), under Assumptions A1, A2a, and
A3a, we seek a deterministic, model based observer that is exponentially stable
for all possible sensor combinations satisfying Assumption A4a.

4.3.2 Observer Design

By copying the system dynamics (4.21), the following observer is proposed:

x̂ = TT ( y )AT( y )x̂ + K( y )ỹ (4.31)

The estimated output is

ŷ = Cy ( y,
 y )x̂ + Dy  (4.32)

and hence when the estimation error is x̃ = x  x̂,

ỹ = y  ŷ = Cy ( y,
 y )x̃ (4.33)

Although it somewhat restricts the exibility, we suggest not to update the Earth-
xed estimates from the acceleration error ỹ3 at this stage, because this would in-
troduce transmission zeros. Therefore, this particular observer gain matrix K( y )
with constant K1i  R12×3 , K2i  R3×3 , K3i  R3×3 and K4i  R3×3 is suggested
 
K11 K12 R( y )T2 0
K21 K22 R( y )T2 0
K( y ) =   (4.34)
K31 K32 R( y )T2 0
K41 R ( y )
T
K42 2T
K43 3T

where the following assumption is made:

A5 Each 3 × 3 sub-block of the gain matrices Kji , 1  j  3, i = 1, 2 commute

with R().

This implies that the gains in North and East must be identical. The body-xed
gain matrices K4i , 1  i  3 can, however, be selected freely.

Error Dynamics

Since y is measured and the constant parameter matrix A is assumed known, we

obtain:
x̃ = TT ( y )AT( y )x̃  K( y )Cy ( y ,  y )x̃ (4.35)
which can be written:
x̃ = TT ( 
y )Ao ( y )T( y )x̃ (4.36)
4.3 Observer With Consistent Wave Model 47

Assumptions A3a, A4a, and A5 are sucient requirements for this. Moreover, it
can be shown that the resulting Ao can be written as

Ao (  y ) = A0 +  y A1 (4.37)
· ¸ · ¸
A0,11 A0,12 0 0
A0 = A1 = (4.38)
A0,21 A0,22 A1,21 0

where A0,11  R18×18 , A0,12  R18×3 , A0,21  R3×18 , A0,22  R3×3 and A1,21 
R3×18 . Denote K̄43 = I  K43 3 such that K43 3 = I  K̄43 . Then:
 
Aw  K11 Cp  K12 2 Cv K11 0
A0,11 =  K21 Cp  K22 2 Cv K21 0  (4.39)
K31 Cp  K32 2 Cv K31 T1 b
 
K12 2
A0,12 =  I  K22 2  (4.40)
K32 2
 ¡ ¡ ¢ ¢T T
 K41 Cp + K42 2 Cv + I  K̄43 Ca
¡ ¢T
A0,21 =  K41 + K̄43 M1 G  (4.41)
¡ ¢T
 K̄43 M1

£ ¤
A0,22 = K42 2  K̄43 M1 D (4.42)
¡ ¢
A1,21 = I  K̄43 SCv (4.43)

Stability Analysis

The form of the observer error dynamics (4.36) is very attractive because the known
transformation T( y ) can be eliminated from the analysis when Assumptions A3-
A5 are employed. Although the eigenvalues of TT Ao (  y )T are identical to the
ones of Ao (  y ) since TT (s) = T1 (s) for all s, Re( i (TT Ao T)) < 0 y if and
only if Ao is Hurwitz. In general, an eigenvalue analysis of a linear time-varying
system will not be sucient to prove stability (Khalil 1996). We have to nd a
Lyapunov function candidate to conclude on that.
The idea is to analyze the error-dynamics in the vessel-xed coordinate system and
selecting a quadratic Lyapunov function candidate V = zT Pz where the P-matrix
also satises some structural constraints. The following lemma will be useful in
that respect (Lindegaard and Fossen 2001b).

Lemma 4.1 Linear time-varying systems on the form

 1 = Ā11 1 + H()Ā12 2 (4.44a)

 2 = Ā21 HT ()1 + Ā22 2 (4.44b)
48 Observer Design

1  Rn1 , 2  Rn2 interconnected by a rotation H : R  Rn1 ×n1 and where

 : R0  R is a known signal and
Ā11 = HĀ11 HT (4.45)
are uniformly globally exponentially stable (UGES) if there for a Q = Q > 0 T

exists a structurally constrained P = PT > 0

· ¸
P11 P12  T H = HT HP
P11 H  11
P= ,  12 = 0 (4.46)
P12 P22
T
H HP
T

such that
PĀ + ĀT P  Q (4.47)
where Ā is the system matrix
· ¸
Ā11 Ā12
Ā = (4.48)
Ā21 Ā22

The system matrix Ā must be Hurwitz, otherwise no such P can be found.

Proof. Dene  = [T1 , T2 ]T . Because of the structural constrains (4.45) on Ā11
we may write
 = TT ()ĀT() (4.49)
where · ¸
HT () 0
T() = (4.50)
0 I
Dene z = T(). Since HT () is a rotation, |z| = ||. Now, abusing the notation

slightly, T() d
= dt (T()), we get
³ ´

z = T() + T() = T() + ĀT() 
³ ´
= 
T()T T
() + Ā z (4.51)

such that dierentiating V = zT Pz along the trajectories yields

V = zT Pz + zT Pz
¡ ¢ ³ ´
= zT PĀ + ĀT P z + zT PTT  T + TT
 TP z
³ ´
 zT Qz + zT PTT  T + TT TP z (4.52)

The structural constraints on P imply that the last term is zero and thus UGES
is proven.
Even though this lemma indeed provides sucient conditions for the elimination
of the kinematics term and the dependence on the varying signal y , we still have
to deal with the time-varying signal  y . Physically  y describes the yaw rate
of the vessel, and intuitively this quantity will be bounded even when exposed
to incoming waves provided that the applied control is appropriate. Therefore,
if a set of simultaneous Lyapunov inequalities are satised at the minimum and
maximum of  y , the error dynamics (4.36) will be ULES. This is summarized in
the following theorem.
4.3 Observer With Consistent Wave Model 49

Theorem 4.1 Consider the observer (4.3)-(4.34) and let the Earth-xed ob-
server gains be selected according to Assumption A5. Assume that:

Ao () = A0 + A1 (4.53)

is Hurwitz for  = 0. If there exists a P = PT > 0

· ¸
P11 P12
P= (4.54)
PT12 P22

where there are structural constraints on P11 and P12

 T H = HT HP
P11 H  11 (4.55)
T 
H HP12 = 0 (4.56)

and an > 0 such that for  = mint (  y ) and ¯ = maxt (  y ) the simultaneous
Lyapunov inequalities are satised

PAo () + ATo ()P   I (4.57)

PAo (¯) + ATo (¯)P   I (4.58)

the error-dynamics (4.36) is uniformly semi-globally exponentially stable (ULES).

Proof. Ao (0) = A0 being Hurwitz is a straightforward requirement, likewise is the

Lyapunov inequalities (4.57)-(4.58) sucient for ensuring that for all   [, ¯] R
since ¡ ¢
f () = xT PAo () + ATo ()P x (4.59)
¯
is linear in  and thus convex such that if f (), f ()   , f ()   for any  
[, ¯].
The rest of the proof consists of verifying that the error-dynamics can be expressed
as (4.44) such that Lemma 4.1 can be employed.

Remark 4.1 For congurations where only position and/or velocity feedback are
used, A1 = 0 and the assumption of  y being bounded is removed. The problem is
thus reduced to nding a suitable P = PT > 0 such that PA0 + AT0 P   I. In
those cases, the observer is UGES.

Remark 4.2 This approach to handling the varying  y is conservative in the

sense that it guarantees exponential stability for arbitrarily fast variations in  y ,
i.e. as long as  y is bounded there is no bound on | ¨ y |.

The simultaneous Lyapunov inequalities can be represented as an LMI feasibility

problem and hence solved using standard software packages: Find a P = PT > 0
in accordance with the structural requirements such that
· ¸
PAo () + ATo ()P 0
<0 (4.60)
0 PAo (¯
) + ATo (¯)P
50 Observer Design

4.3.3 Concluding Remarks

An existing nonlinear model-based observer with wave ltering capabilities for

surface vessels has been extended to optionally include velocity and acceleration
measurements. If the environmental model and some of the gain matrices sat-
isfy certain structural properties, global exponential stability of the lter error
dynamics can be concluded using a quadratic Lyapunov function with structural
constraints. Once these model and gain constraints are satised, stability is de-
termined by an eigenvalue analysis of the observer error’s system matrix Ao ().
Without feedback from acceleration, the observer error is globally exponentially
stable when Ao (0) = A0 is Hurwitz. With acceleration feedback, on the other
hand, the observer error will be exponentially stable semi-globally.

4.4 Simplied Observer

In the previous section we introduced a model based observer with wave ltering
capabilities for surface vessels at low speed. Although this observer was the rst
model based integrated design that could incorporate partial velocity and accel-
eration measurements, there are two reasons why this design is unsuited from a
practical point of view. At the core of these problems is the suggested model of
rst order wave induced motion:

1. The tuning procedure is much more complicated when velocity and/or accel-
eration measurements are included compared to the pole placement strategy
used for position measurements. This is due to the fact that in the general
case the observer gains enter non-anely in the expressions describing the
eigenvalue of the observer error-dynamics. One solution is to solve an alge-
braic Riccati equation (Kalman gains or H -ltering techniques) either a
priori or on-line, but having complete control of the notch-eects is almost
impossible. As a consequence, it is very likely that the time spent tuning
the DP system during sea-trials will increase.

2. A common wave model for all state derivatives could be fatal for the stability
of the combined wave motion model if the individual measurements are out
of synchronization with respect to each other. This occurs e.g. if the time-
delays from the sensor system components are dierent.

Here we propose a model based observer where the wave models for position,
velocity, and acceleration measurements are considered separately. The main idea
is that wave induced acceleration is “uncorrelated” with the induced velocity, an
assumption that is motivated more from engineering experience rather than from
physics. Global exponential stability of the error dynamics may still be guaranteed
using a structured P matrix. Still, we chose to trade global results in order to relax
Assumption A2a and the need of a positive denite Tb . With a bounded yaw rate
the stability results will be valid only in a semi-global sense.
4.4 Simplied Observer 51

4.4.1 Complete Ship and Environment Model

As for the previously proposed observer, the LF vessel and bias models are de-
scribed by (4.3). The by assumption uncorrelated wave induced positions, veloci-
ties, and accelerations respectively are given by

p w = Apw pw + Epw wpw (4.61)

v w = Avw vw + Evw wvw (4.62)
a w = Aaw aw + Eaw waw (4.63)

where the order of each wave model number is arbitrary, but it is recommended
to keep the order fairly low. Second or fourth order linear models are sucient.
Let mp , mv , ma denote the order of the position, velocity and acceleration wave
models respectively. Then pw  R3mp ,vw  Rny2 ·mv , aw  Rny3 ·ma describe
the rst order wave-induced positions, velocities and accelerations respectively.
Apw  R3mp ×3mp , Avw  Rny2 ·mv ×ny2 ·mv , Aaw  Rny3 ·ma ×ny3 ·ma are assumed
Hurwitz and describes the rst order wave induced motion. The wave and bias
models are driven by disturbances of appropriate dimensions.
In order to make use of the commutation properties, we have to assume

A2b The bias time constant matrix Tb and each 3 × 3 sub-block of Apw satisfy
Property 4.1.

Now, collect all the Earth-xed states in x1  R6+3mp and stack the body-xed
ones into x2  R3+ny2 ·mv +ny3 ·ma
£ T ¤T
x1 = pw T bT (4.64)
£ T ¤
T T
x2 = vw aw 
T
(4.65)

Let n denote the dimension of x = [xT1 , xT2 ]T and dene the block diagonal trans-
formation matrix T : R  Rn×n

T( y) = Diag(RT ( y ), · · · , RT ( y ), I3+ny2 ·mv +ny3 ·ma ) (4.66)

On compact form using Assumption A2b and w = [wpw

T
, wbT , wvw
T
, waw
T T
] we get

x = TT ( y )AT( y )x + B + Ew (4.67)

where the parameters A have been separated from the rotation R( y ). The model
parameters in (4.67) are
 
Apw 0 0 0 0 0
0 0 0 0 0 I
0 0 T1 0 0 0
A= b (4.68)
0 0 0 Avw 0 0
 0 0 0 0 Aaw 0 
0 M1 G M1 0 0 M1 D
52 Observer Design
   
0 Epw 0 0 0
0 0 0 0 0
0 0 Eb 0 0
B= , E= (4.69)
0 0 0 Evw 0
 0   0 0 0 Eaw 
M1 0 0 0 0

Measurements

Position and heading measurements are always required, and the number of ve-
locity and acceleration measurements available are denoted 0  ny2  3 and
0  ny3  3, respectively. Let y1  R3 , y2  Rny2 and y3  Rny3 describe the
position, velocity and acceleration measurement vectors. We dene the measure-
ments as

y1 =  + Cpw pw (4.70)
y2 = 2  + Cvw vw (4.71)
y3 = 3  + Caw aw (4.72)

where, 2 and 3 are projections isolating the components of the LF-model that
are actually measured. Written compactly,

y = Cy ( y )x + Dy  (4.73)

where

Cpw I
Cy ( y) =  2 Cv RT ( y) 0
0 3 M1 GRT ( y)

0 0 0 0
0 2 Cvw 0 
3 M1 RT ( y ) 3 M1 D 0 Caw

£ ¤T
Dy = 0 0 MT T3 (4.74)

Physically, however, it should be pointed out that the LF linear accelerations that
are being measured is not  as claimed in (4.72) since 
¨ 6=  when the Earth-xed
frame is considered as being the inertial frame. More specically, considering the
LF dynamics
y3,LF = ¨ =  y SR( y ) + R( y ) (4.75)

which means that (4.72) is approximately correct for small angular rates. For large
angular rates, however, an auxiliary pre-processor should be used to compensate
for the Coriolis eect  y SR( y ). The need for an external processing unit will
in fact always be there as discussed Section 3.3.
4.4 Simplied Observer 53

4.4.2 Observer Design

By duplicating the model dynamics and introducing a low-pass lter in order to

achieve a certain roll-o eect, see (4.101), the following model based observer is
proposed
a f = T1
f (af + ỹ3 ) (4.76)

x̂ = TT ( y )AT( y )x̂ + B + K( y )ỹ + Kf af (4.77)

and its estimated output is

ŷ = Cy ( y )x̂ + Dy  (4.78)

and hence when the estimation error is x̃ = xx̂, the output error is ỹ = Cy ( y )x̃.

A pragmatic selection of observer gain matrices K( y ) and Kf reducing intercon-

nections is  
K11 0 0
K21 0 0
K31 0 0
K( y ) = (4.79)
0 K42 0
 0 0 K53 
K61 RT ( y ) K62 0
£ ¤T
Kf = 0 0 0 0 0 KTa (4.80)
In order to apply the concept of commutating matrices, we have to impose the
following requirement on some of the gains:

A3b Each and every 3 × 3 block of K11 , K21 and K31 commute with the rotation
R( y ) (Property 4.1).

A schematic drawing of this observer, without bias estimation, is given in Figure

4.1.

Stability Analysis

When Assumption A2b and A3b are satised, the rotations can be separated
from the parameters. The observer error-dynamics can hence be rewritten on the
compact form

x̃ = TT ( y )Ao T( y )x̃ + Kf af + Ee w (4.81)

a f = Tf af + T1
1
f C3 T( y )x̃ (4.82)
· ¸
A11 A12
Ao = (4.83)
A21 A22
54 Observer Design

Figure 4.1: Observer with rst order cut-o lter for acceleration. Bias estimation
is not shown.
 
Apw  K11 Cpw K11 0
A11 =  K21 Cpw K21 0  (4.84)
1
K31 Cpw K31 Tb
 
0 0 0
A12 =  0 0 I  (4.85)
0 0 0
 
0 0 0
A21 = 0 0 0  (4.86)
K61 Cpw M G  K61 M1
1
 
Avw  K42 Cvw 0 K42 2
A22 =  0 Aaw  K53 Caw K53 3  (4.87)
K62 Cvw 0 M1 D  K62 2
£ ¤
C3 = 0 3 M1 G 3 M1 3 M1 D 0 Caw (4.88)

Stacking (4.81)-(4.82) together into z  Rnz , more specically z = [x̃T , aTf ]T , we

then get
z = TTz ( y )Az Tz ( y )z + Ez w (4.89)
where Tz = Diag(T, Iny3 ) and
· ¸ · ¸
Ao Kf Ee
Az = , Ez = (4.90)
T1
f C3 T1
f 0

We now state a robustness-like theorem for the stability of this lter. The limiting
factor is the yaw rate  y = ry , and we could just as well repeat using a P-matrix of
4.4 Simplied Observer 55

a certain structure that commutes with Tz ( y ) in a quadratic Lyapunov function

provided that Assumption A2b and A3b hold. Because the upper bound on |ry (t)|
is likely to be larger than the physical limit, we assign an arbitrary P. The
advantage of an arbitrarily selected P is that there are no longer any restrictions
on the selection of cross-terms. As a consequence we could have let kTb k  ,
then the bias is modelled as an open integrator and we obtain true integral action.
Before we state the theorem we need to introduce a skew-symmetric matrix Sz
that appears when the rotation Tz ( y ) is dierentiated.

4 d ¡ ¢
z =
T Tz ( y) =  y Sz Tz ( y) =  y Tz ( y )Sz (4.91)
dt
In our case Sz = Diag(ST , ..., ST , 0nz 12xnz 12 ) where S is given by (4.2).

Theorem 4.2 The observer error dynamics (4.89) is exponentially stable for small
|ry (t)| < rmax (ULES) if and only if Az is Hurwitz. Suppose an rmax > 0 is explic-
itly given, then (4.89) is uniformly globally exponentially stable (UGES) if there
exists a matrix P = PT > 0 such that the following two LMIs are feasible for some
>0 ¡ ¢
PAz + ATz P + I  rmax ¡PSz + STz P¢
(4.92)
PAz + ATz P + I  rmax PSz + STz P

Notice that stability can be characterized without dealing with the rotations
Tz ( y ). However, there is a bound on the rotation rate rmax making the observer
USGES due to observability.
Proof. Use the non-singular rotation Tz ( y) as a mapping  = Tz z. Then,

  z z + Tz z
= T
=  y Sz Tz z + Tz TTz Az Tz z
³ ´
= Az +  y Sz  (4.93)

Consider the Lyapunov function V = T P

¡ ¢ ¡ ¢
V = T PAz + ATz P  +  y T PSz + STz P  (4.94)

Since V is linear in  y for xed P and  it is also convex and it suces to verify
that V < 0 at the boundaries of  y , namely ±rmax since rmax   y  rmax by
assumption. We therefore have to make sure that
¯
¯
V 1 = V ¯ <0 (4.95)
 =rmax

¯ y
¯
V 2 = V ¯ <0 (4.96)
 y =rmax

Inserting the LMIs from (4.92) we get for k = 1, 2

V k   I (4.97)
56 Observer Design

For small enough rmax , there will always exist an > 0 if and only if Az is
Hurwitz.
This result could also be proven by the circle criterion, but the maximum allowable
rmax is likely to be smaller due to the required SPR-property.
Notice that stability can be proven also for the limiting case T1
b = 0. The
passive design (Fossen and Strand 1999) and our previous version (Lindegaard
and Fossen 2001a), on the other hand, require T1
b > 0.

4.4.3 Observer Tuning

In this section we suggest models for the rst order wave loads and then we suggest
tuning rules that based on those models generate the desired frequency response
between the measurements and the LF estimates.
For position, velocity, and acceleration measurements, i = p, v, a, a cascade of
second order linear systems
· ¸
0 I £ ¤
Aiw = , Ciw = 0 I (4.98)
i i

can be used to represent the wave induced motion whereby we obtain the desired
wave ltering capability. Treat each DOF separated from the others by setting

i = diag( 2i,1 , ...,  2i,nyi ) (4.99)

i = diag(2 i,1  i,1 , ..., 2 i,nyi i,nyi ) (4.100)

where  i,k > 0 is the resonance frequency and i,k > 0 is the relative damping
factor which determines the width of the spectrum.
Depending on the number pi = mi /2, where i = p, v, a of second order models
in cascade, the desired transfer function between any measurement and the LF
estimate is
³ ´pi
s2 + 2 i,k i,k s + 2i,k
hdi (s) = c,k ³ ´pi i = p, v, a (4.101)
s2 + 2 i,k i,k i,k s +  2i,k ( c,k + s)

which is a notch-lter, with center frequency at  i,k , the wave model resonance,
and notch “width” given by  i,k 1, in series with a low-pass lter that guarantees
a certain roll-o for frequencies larger than  c,k . In order to achieve good perfor-
mance, the roll-o frequency  c,k should be larger than the resonance frequency
of the notch-lter, that is c,k  i,k .

Wave Model Gains

We apply a same pole-placement technique to nd observer gains for position

and velocity innovation. If we were dealing with second order wave models, the
4.4 Simplied Observer 57

tuning rules from Fossen and Strand (1999) apply. However, due to the increasing
power of the wave frequency components in the velocity and acceleration signal,
we suggest using a fourth order wave model, at least for acceleration, in order to
achieve satisfactory wave ltering capabilities. We were in fact unable to get good
results for acceleration using second order models. Below, we therefore present
the extension of Fossen and Strand (1999) to fourth order models.
Consider the surge dynamics being updated from position measurements. We aim
to nd the elements to put in K11 and K31 in order to create the desired notch
and roll-o hdi (s) in (4.101)
ˆi
(s) = hdi (s) (4.102)
y1,i
This can be obtained one degree of freedom at the time by dening p,i =  2p,i > 0
and p,i = 2 p,i p,i > 0 and letting
 
0 1 0 0
p,i  p,i 0 1
Apw,i =   (4.103)
0 0 0 1
0 0 p,i  p,i
£ ¤
C̄pw,i = 1 0 0 0 (4.104)

and selecting observer gains according to

k11,i = L1
i ci (4.105)
k31,i =  c,i (4.106)

where  
1 0 0 0
2 p,i 1 0 0
Li =  (4.107)
p,i + 2p,i p,i 0 1 
p,i 1 1 0
 
¡ 2 ¢p,i ( p,i  1)
2 2
p,i  p,i  1 + 2 p,i ( ¡p,i  1) ¢ c,i
ci =  (4.108)
p,i p,i ( p,i  1) + 2p,i  2p,i  1 c,i 
2 p,i ( p,i  1) c,i
ensures that the specied notch-eect and roll-o is indeed acquired. The gains
K31 and K61 , the gains from position innovation which update the bias and LF
velocity, can be selected freely as long as Az remains Hurwitz.
The very same approach can be applied to assign values to K42 , K62 in order to
obtain a notch-eect for the velocity measurements.

Acceleration Gains

The acceleration part of the lter possesses another feature as well. Measuring the
acceleration could be regarded as an alternative to using a model based observer
58 Observer Design

because the model actually estimates the acceleration while an accelerometer mea-
sures it. The gain Kf serves as a weight factor determining how much emphasis
we should put on the model. When Kf = 0 we choose not to utilize acceleration
feedback to update the lter at all and when Kf = 1, the LF model description is
completely disregarded for low frequencies.
The low-pass lter between acceleration innovation ỹ3 and ˆ  takes care of the
roll-o. The lter constants Tf should therefore be selected as
T1
f = diag( c,1 , ...,  c,ny3 ) (4.109)

Next, to obtain the desired notch-ltering around the resonance frequency, select
· ¸
0
K53 = (4.110)
diag( a,1 , ...,  a,ny3 )

4.4.4 Experiments

The experiment was carried out with “Cybership II”, and the dynamic compen-
sator scheme presented in Chapter 2 was implemented.
Based on the principle of certainty of equivalence, an observer-feedback PID-like
tracking controller on the form

 = 
ˆ  d (4.111)
 = Ki RT ( ˆ )  Kp RT ( ˆ ) (ˆ
  d )
   d)
Kd (ˆ (4.112)

was used to keep the boat on the position  d = [0.3, 0, 0]T ,  d = 0. The controller
and the thrust allocation algorithm is described and analyzed in Lindegaard and
Fossen (2003).
From t  20 seconds and onwards, the model ship was exposed to JONSWAP-
distributed irregular head waves. The peak period and signicant wave height
were set to Ts = 0.75 and H1/3 = 0.02 meters respectively.
Time series plots of the measured positions (dotted) and their respective LF esti-
mates are reproduced in Figure 4.2 together with the observer’s surge bias estimate.
Notice that the surge bias converges towards the controller’s I-term, that is the
mean of applied surge propeller force  1 . A large wave slammed into the vessel at
t  115 generating a temporary drift o in East and heading because the vessel
had a small oset angle at the time of the impact. The slow oscillations are due
to nonlinear wave eects and not to the rst order induced motion. Figure 4.3
shows g-compensated measurements of the surge and sway accelerations. Here,
the wave frequency motion (rst order wave loads) dominate the picture. But as
the empirical transfer functions (Figure 4.4) of the measured signals and the state
derivatives show, for low frequencies the estimated LF-accelerations are excellent,
because they follow the measured signals at frequencies below f = 0.1 Hz. As
required, frequency components around the wave frequency peak f = 1/Ts = 1.33
Hz have been successfully attenuated.
4.4 Simplied Observer 59

North pos. East pos.

-0.24 0.06

-0.26
0.04
-0.28

-0.3 0.02
[m]

[m]
-0.32 0
-0.34
-0.02
-0.36

-0.38 -0.04
0 50 100 150 0 50 100 150
Time [s] Time [s]
Heading ψ Surge/North bias and -mean( τ )
y 1
6 0.05

4
0
Force [N]

2
[deg]

-0.05
0

-0.1
-2

-4 -0.15
0 50 100 150 0 50 100 150
Time [s] Time [s]

Figure 4.2: Top left, measured (dotted) and LF estimated North position. Top
right, measured (dotted) and LF estimated East position. Bottom left, measured
(dotted) and LF estimated heading. Bottom right, estimated bias and mean of
applied thrust thrust  1 .

4.4.5 Concluding Remarks

A simple model based state estimator for surface vessels with wave ltering capa-
bilities has been proposed and analyzed along with an intuitive tuning procedure.
For bounded yaw rate, the observer error dynamics was shown to be exponentially
stable. Inertial measurements, that is linear accelerations and yaw rate, were in-
cluded in the lter to improve performance. Due to the acceleration measurement
ambiguity, a g-compensation system had to be utilized in order to remove gravity
components from the linear acceleration terms.
Experimental results with a model ship performing a DP operation as it was
exposed to incoming irregular waves illustrated the performance of the lter. Em-
pirically calculated frequency responses between available measurements and esti-
mated low frequency positions, velocities and accelerations documented that the
desired notch ltering of rst order wave induced motion was achieved.
60 Observer Design

Surge accelerations
0.75

0.5

0.25
[m/s ]
2

-0.25

-0.5

-0.75
0 20 40 60 80 100 120 140
Time [s]

Sway accelerations
0.15

0.1

0.05
[m/s ]
2

-0.05

-0.1

-0.15
0 20 40 60 80 100 120 140
Time [s]

Figure 4.3: g-compensated surge and sway accelerations. Measurements (dotted)

and LF estimates (solid).
4.4 Simplied Observer 61

TF x- and y-position TF heading ψ (solid) and yaw rate (dotted)

20 20
0 0

-40 -40
Gain [dB]

Gain [dB]

-80 -80

-120 -120

-2 -1 0 1 -2 -1 0 1
10 10 10 10 10 10 10 10
Freq. [Hz] Freq. [Hz]

TF surge- and sway-acceleration TF position to uhat and vhat

20 20
0 0

-40 -40
Gain [dB]

Gain [dB]

-80 -80

-120 -120

-2 -1 0 1 -2 -1 0 1
10 10 10 10 10 10 10 10
Freq. [Hz] Freq. [Hz]

Figure 4.4: Top left, TF from measured North (solid) and East (dotted) position
to their respective LF estimates. Top right, TF from measured heading (solid) and
yaw rate (dotted) to LF estimates. Bottom left, TF from measured surge (solid)
and sway (acceleration) to LF estimates. Bottom left, TF from North position to
surge velocity (solid) and from East to sway (dotted).
62 Observer Design
Chapter 5

Controller Design

5.1 Introduction
Problems of motion control can be classied into three groups (Encarnação and
Pascoal 2001):

• Point stabilization: The objective is to stabilize and keep a vehicle at a

specied point and orientation.
• Trajectory tracking: The task of making a vehicle track a reference trajectory
parameterized in time.
• Path following: The vehicle is required to converge to and follow a desired
path without an implicit speed assignment.

Trajectory tracking is the most commonly implemented control approach in com-

mercial DP systems today. This chapter focuses on both linear and nonlinear
control strategies for low speed tracking and positioning control (point stabiliza-
tion) of fully or overactuated surface vessels in the horizontal plane. Whenever
the desired velocities are zero, the desired trajectory collapses into a single point
and the tracking controller becomes a positioning controller. There are, however,
operations where starting time and encountered delays are of minor concern. An
intuitive example is way-point tracking, the task of following a path of specied
way-points in the horizontal plane. Furthermore, more typical DP operations like
pipe lying and dredging are also likely to benet from a control strategy with less
attention to time constraints.
Path following is typically applied to underactuated ships in transit, and many
recent developments have been reported: In Zhang et al. (2000) the authors ad-
dressed the problem of following straight lines, and to achieve this a new output
combining the cross-track error and the heading angle was dened. Asymptotic
stability to the path was shown by the use of sliding mode control. Another rede-
nition emulating that of an experienced helmsman was proposed by Pettersen and
64 Controller Design

Lefeber (2001); the desired heading was dened as a function of the cross-track
error. A more exible alternative is the introduction of a Serret-Frenet frame to
represent the cross-track and heading error (Encarnação et al. 2000, Encarnação
and Pascoal 2000, Skjetne and Fossen 2001). The advantage is that the path can
be regarded as a more general smooth curve in the plane rather than consist of
straight lines.
Underactuated trajectory tracking and stabilization has received a lot of attention,
see Do et al. (2002) and references therein. A natural consequence of the ship being
underactuated is that geometric restrictions apply on the reference trajectory.
More specically, the desired yaw velocity must be persistently excitative.
A method combining trajectory tracking and path following was proposed in Hind-
man and Hauser (1996). Suppose a desired path  d : R0  Rn is given and that
it is a continuous function of the path variable  0. The basic idea is to determine
 by projecting the current state of the vehicle onto the reference trajectory. This
returns the appropriate trajectory “time” given the current state of the system.
Under some geometrical path conditions, it is shown that feeding  d () instead
of  d (t) into an already existing tracking controller guarantees convergence to the
path. This procedure, however, requires a path specication for the full state, and
it is applicable to feedback linearizable systems. An output maneuvering extension
to Hindman and Hauser (1996) was proposed in Encarnação and Pascoal (2001).
By employing backstepping (Krstić et al. 1995), the need for time derivatives of
 d , the full path specication, was relaxed. This approach is well suited for me-
chanical systems such as ships, but it is less tted for systems of relative degree
higher than two due to the need of higher order derivatives of the path variable .
Recently, Skjetne et al. (2003) proposed a more general robust maneuvering design
for systems on strict feedback form which tackles the relative degree restriction.
Trajectory tracking is geometrically speaking relatively simple compared to path
following, and the above described methods unifying tracking and path following
promise increased exibility of the controllers derived in this chapter. We extend
existing DP designs (Strand 1999, Berge 1999) and provide ideas and suggestions
to improve and facilitate the design and implementation of observer feedback po-
sitioning control systems. Ships suited for traditional dynamic positioning opera-
tions are usually overactuated. This means that the desired thrust  can almost
always be satised. In fact, due to propeller rate and thrust magnitude constraints,
there are indeed practical limits, but it is assumed that these issues can be disre-
garded or at least addressed elsewhere. The control laws will be derived assuming
full state feedback. Later substituting the state variables with their respective low
frequency estimates, the so called principle of certainty of equivalence, we show
that the combination of an observer and controller stabilizes the entire closed loop
system. Results from the study of nonlinear cascaded systems will be used in the
analysis.
These are the main objectives of this chapter:

1. Derive a simple yet exible framework for low speed trajectory tracking
of fully actuated surface vessels. First, a method separating the nonlinear
5.2 Trajectory Generation 65

kinematics from the otherwise linear closed-loop dynamics will be proposed.

This type of design, which has many similarities to more conventional PID-
compatible schemes, is attractive because one can apply any linear design
tool to shape the error-dynamics. Some examples will be given. The method
called commutating design will be described in Section 5.3. A generalization
to more nonlinear ship models (higher velocities) is given in Section 5.5.

2. To rene the already derived framework in order to incorporate acceleration

terms in the controller (Section 5.4).

5.2 Trajectory Generation

The desired position or trajectory of a DP vessel is usually given with respect to the
so called center of rotation (COR), about which any rotation of the vessel should
be conducted. The COR is a freely selectable xed point inside or outside the
physical boundaries of the ship, and its coordinates will be assumed given relative
to the body-xed b-frame. This section derives some basic kinematic relations
between the desired trajectory given in Earth-xed coordinates (n-frame) and
the transformation into the b-frame which is the reference frame where the DP
controller operates.

Figure 5.1: Denition of the Earth-xed position of center of rotation rnnr and its
relation to the position of the vessel rnnb .

Let rbbr = [xbbr , ybr

b T
] denote the position of the COR relative the the b-frame.
The desired Earth-xed position and heading of the COR is given by  dr (t) =
66 Controller Design

[(d rnnr (t))T , d (t)]T where again d rnnr = [d xnnr ,d ynr

n T
] is the distance in the hori-
zontal plane from the n-frame to the desired position of the COR. Observe that
the desired heading d (t) of the COR is identical to the desired heading of the
ship. According to Figure 5.1 the desired position of the b-frame is given by

rnb
d n
= d rnnr  R̄( b
d )rbr (5.1)

where the rotation R̄( d) is dened as

· ¸
cos  sin
R̄( d) =
d d
(5.2)
sin d cos d

Frequently, it is dersirable to let the vessel rotate about some other point than
CG. For instance when deploying a device on the sea bed using a crane, the ship
should rotate about the crane head rather than CG. Consider Figure 5.2 where
a 180 deg change of heading instructs the DP controller to move the vessel on a
circular arc while turning rather than revolving about the center of the b-frame.

Figure 5.2: Change of heading 180 deg about the COR: The ship is moved from
rnnb (t0 ) towards rnnb (tf ) while turning.

The Earth-xed velocity and acceleration of the COR are denoted  dr (t) = [(d r nnr (t))T , rd (t)]T
and 
¨dr (t) = [(d r̈nnr (t))T , rd (t)]T respectively. Desired velocity and acceleration of
5.3 Commutating Control 67

the b-frame decomposed in the n-frame are found by dierentiating d rnnb .

b
r nb
d n
= r nr
d n
 rd R̄( d )S̄rbr (5.3)
 rd2 R̄( 2 b
b
r̈nb
d n
= r̈nr
d n
 rd R̄( d )S̄rbr d )S̄ rbr

+ rd2 R̄(
b
= r̈nr
d n
 rd R̄( d )S̄rbr
b
d )rbr (5.4)

Here · ¸
0 1
S̄ = (5.5)
1 0

The desired speed and acceleration along the trajectory should not exceed neither
the physical nor the imposed limitations of the ship. Considering these limita-
tions in an Earth-xed setting is dicult, but the control will be regained when
expressing the trajectory in a reference parallell frame, the d-frame: When the
vessel tracks the desired trajectory perfectly, the velocities in the d-frame will be
exactly those of the vessel itself (b-frame). The desired velocity is
4
vnb
d d
= R̄T ( d) r nb
d n
b (5.6)
= R̄T ( d) r nr
d n
 rd S̄rbr

and for the acceleration we have

d ¡d d ¢
v nb
d d
= vnb
dt
b
= R̄T ( d )d r̈nnr  rd R̄T ( d )S̄ r
 nr
d n
 rd S̄rbr (5.7)

The translation and rotation of the d-frame can thus be regarded as a virtual ship,
yet the motion does not need to be that of a numerical ship model.
To summarize, suppose the COR is located at rbbr = [xbbr , ybr b
, 0]T , and once a
smooth, Earth-xed reference trajectory for the COR is given,  dr : R0  R3 ,
 dr : R0  R3 , and 
¨ dr : R0  R3 , the desired trajectory for the orgin of the
vessel will be given by

d =  dr  R( d )rbbr (5.8)
d = RT ( d ) dr  rd Srbbr (5.9)
 d = RT ( d )¨
 dr  rd RT ( d )S
 dr  rd Srbbr (5.10)

Observe that this trajectory satises  d = R( d ) d . The corresponding velocities

 d (t) = [ud , vd , rd ]T and accelerations  d (t) = [u d , v d , rd ]T are thus decomposed
in the reference parallel d-frame.

5.3 Commutating Control

From a practical point of view it is important that the tuning procedure for the
DP controller is intuitive. There is usually little time available for detailed tuning,
68 Controller Design

and still the adjustments made during a short sea trial have to perform well in
weather conditions ranging from calm sea to the more extreme. During the limited
time available personnel with limited theoretic background make decisions inu-
encing the vessel’s positioning performance for many years to come. Computer
simulations do provide an acceptable initial tuning, but some modications are
almost always needed. If a controller gain adjustment does result in a predictable
behavior of the vessel, it is likely that the engineer is condent that the tuning
is acceptable, the overall procedure takes less time, and the customer eventually
gets satised.
Here we elaborate on nonlinear PID-like tracking controllers simply because such
controllers are readily interpreted and analyzed. Low speed vessel tracking and
positioning control can be accomplished with relatively simple means, such as
PID-control, and in a critical situation where understanding of physics is needed
experimenting with other approaches may not be advisable. The price we pay
for concentrating on such straightforward controllers is that stability can only be
guaranteed for bounded yaw rates. However, the resulting upper bound rmax for
a well-behaving controller usually exceeds the physical limitations for the ship.
A denite advantage is that we can incorporate static or dynamic feedback from
acceleration directly within the derived framework. The derived controllers are
linear in the sense that their respective terms are bounded linearly in the error
variables. They are nonlinear in the sense that the kinematics is included.

5.3.1 Full State Feedback PID Tracking Control

The motivation for studying PID or even PD in detail is to show how dierent
linear design techniques can be applied in DP control. Similarly to the observer
design, it is possible to separate the kinematics from the design procedure and
later include these rotations again when implementing the controller. Finding the
controller’s gains themselves is a task which may be performed without considering
the kinematics at all.
The objective is to track a given smooth trajectory (d (t),  d (t),  d (t)). Let the
position error  e be given in the Earth-xed n-frame and the velocity and accel-
eration error in the b-frame, that is

e =   d (5.11)
e =   RT ( e ) d (5.12)
 e =    e ST RT ( e ) d  RT ( e )
d (5.13)

Observe that in DP  d =  d = 0 such that

e =   d
e = 
 e = 
5.3 Commutating Control 69

The low-speed ship model considered is

 = R( ) (5.14)
M + DL  =  + RT ( )w (5.15)

where w = wb + R( )wsv is an disturbance consisting of a constant bias wb and a

slowly varying term wsv . The former will be attenuated be integral action. During
the design we will, however, let w = 0 in order to simplify the stability analysis.
A controller with respective proportional, integral, and derivative gains Kp , Ki , Kd 
R3×3 can be formulated as

 = e (5.16)
 = Ki RT ( )  Kp RT ( ) e  Kd  e +  r (5.17)

where e =  d . Notice that the rst three terms form the PID feedback
control while the latter term  r is the reference feed-forward given by
³ ´
 r = DL RT ( e ) d + M  e ST RT ( e ) d + RT ( e ) d (5.18)

Also note that the gains Ki and Kp have been put to the left of RT ( ) making
them body-xed gains. This makes more sense to an operator than having them on
the right of RT ( ) (Loría et al. 2000) because the compass heading (t) should not
inuence the convergence rates. The drawback of doing this is that energy-based
stability proofs can no longer be applied directly.
By inserting  into (5.15) we get
³ ´
M    e ST RT ( e ) d + RT ( e ) d =Ki RT ( )Kp RT ( ) e (Kd + DL )  e

Consequently, the velocity error dynamics can be written

M e = Ki RT ( )  Kp RT ( ) e  (Kd + DL )  e (5.19)

By collecting the states xe = [T ,  Te ,  Te ]T and simultaneously considering the

slowly varying disturbance wsv only, we may now rewrite the closed-loop dynamics
on the compact form

x e = TT ( )Ac T( )xe + Bwsv (5.20)

where T( ) is dened as the block-diagonal T( ) = Diag(RT ( ), RT ( ), I) and

 
0 I 0
Ac = A  BK =  0 0 I  (5.21)
M1 Ki M1 Kp M1 (DL + Kd )
£ ¤T
B = 0 0 MT (5.22)
£ ¤
K = Ki Kp Kd (5.23)

Due to the pair (B, A) being controllable, we have complete control over the
eigenvalues of Ac by assigning appropriate gains Kp , Ki , and Kd . Similarly to
70 Controller Design

the observer design, the eigenvalues of TT ( )Ac T( ) are constant for all and
equal to the ones of Ac because TT ( ) = T1 ( ). The applied control  is using
these denitions written as

 = KT( )xe +  r (5.24)

where the gain K  R3×9 is given by (5.23) and  r by (5.18).

We have deliberately left out a Coriolis term in the controller (5.24) compensating
for the time-derivative of R( ) that would have allowed us to establish global sta-
bility properties. The limiting factor is thus the rotation rate  = r. Nevertheless,
uniform local exponential stability (ULES) can be veried for bounded r(t), that
is |r(t)|  rmax . Notice also that limt xe (t) = 0 does not imply that  0 but
rather  d even though occurs in error dynamics x e = TT ( )Ac T( )xe .

Theorem 5.1 Consider the system (5.4)-(5.5) controlled by (5.24). Suppose

|r(t)|  rmax and wsv = 0 for all t t0 . The origin xe = 0 of (5.20) is uniformly
locally exponentially stable provided rmax > 0 is suciently small and if and only
if Kp , Kd , Ki  R3×3 are chosen such that Ac as dened by (5.2) is Hurwitz. If
rmax is larger than any physical upper limit for |r(t)|, (5.20) is said to be uniformly
globally exponentially stable.

Proof. The necessity of Ac being Hurwitz is well known. For proving suciency,
 ) we hereby mean d (T( )). Then,
dene z = T( )xe . By T( dt

z = Tx  T z + Ac Txe = (Ac + rST ) z

 e + Tx e = TT (5.25)

where ST = Diag(ST , ST , 0). If and only if Ac is Hurwitz there exists a P = PT >

0 such that
PAc + ATc P = Q (5.26)
for a given Q = QT > 0.
Consider the radially unbounded storage function V (xe , ) = xTe TT ( )PT( )xe =
zT Pz. Dierentiated along the trajectories we get

V = z T Pz + zT Pz
¡ ¡ ¢¢
= zT ATc P + PAc + r PST + STT P z
2
 zT Qz + 2rmax max (P) |z|
2
= ( min (Q)  2rmax max (P)) |xe | (5.27)

which is negative denite provided that rmax is small enough.

If it is known a priori that r(t)  L2 L , convergence to zero can be shown
no matter how large supt0 r(t) actually is. On the other hand, if we do not in
advance know that r(t)  0, an ecient but admittedly conservative method for
estimating rmax > 0 is solving a generalized eigenvalue problem (Boyd et al. 1994)
as explained in the following Corollary:
5.3 Commutating Control 71

Corollary 5.1 For a given Ac , a yaw rate bound rmax > 0 that guarantees expo-
nential stability of (5.20) can be found be solving the following generalized eigen-
value problem in the decision variables P and

minimize
P = PT > 0, >0 ¡ ¢ (5.28)
subject to PST + STT P <  ¡ATc P + PAc ¢

PST  STT P <  ATc P + PAc

where rmax = 1/ .

Proof. From the proof of Theorem 5.1 we note that if

¡ ¡ ¢¢
g(r) = max ATc P + PAc + r PST + STT P < 0

for all t t0 then (5.20) is exponentially stable. As g(r) is linear in r it suces

to show g(rmax ) < 0 and g(rmax ) < 0, and the constraints in (5.28) therefore
ensures g(r) < 0 for all |r(t)|  rmax .
Observe that Corollary 5.1 yields a conservative rmax because it allows innitely
fast changes in r(t), that is unbounded |r(t)|,
 as long as |r(t)| < rmax . Still, nding
this rmax given any Hurwitz Ac is a relatively simple task using available software
packages like Matlab’s LMI Toolbox (Gahinet et al. 1995).

5.3.2 LMI Control Strategies

Yet another feature of the LMIs given by Corollary 5.1 is that it may be combined
by LMI based control synthesis methods in order to provide state-feedback con-
trollers rendering the closed loop globally exponentially stable for any specied
rmax . This section briey explores how LMI synthesis provides a state feedback
K resulting in a globally exponentially stable closed loop dynamics for any given
rmax .
First, we briey summarize three common and widely applicable kinds of LMI
control strategies that may be used separately or in combination, those being H ,
H2 , and pole-clustering. When applied together, we say that the resulting control
is multiobjective. Using the more general description of a linear plant

x = Ax + Bu + Ew
H z = C x + Dw w + Du u (5.29)

z2 = C2 x + D2u u

where z  Rn and z2  Rn2 are the outputs used in the H and H2 cost
criteria, respectively. Let A  Rn×n and B  Rn×m be the nominal model,
and without any further discussion on model scaling, we assume that E  Rn×p ,
C  Rn ×n , Dw  Rn ×p , Du  Rn ×m , C2  Rn2 ×n , and D2u  Rn ×m
are properly scaled. The control objective is to compute a state-feedback controller

u = Kx (5.30)
72 Controller Design

that fullls certain specications on the closed-loop behavior. These closed loop
design criteria can be cast as convex optimization problems satisfying certain LMIs.
The following summary is taken from a variety of sources (Boyd et al. 1994, Chilali
and Gahinet 1996, Scherer et al. 1997).

H -Control

Let Hwz (s) denote the closed loop realization between w and z . Minimization
of the H -gain  from w to z can be cast as an LMI optimization problem in
the matrix variables X  Rn×n and Y  Rm×n (with X = XT > 0) and Y
while minimizing  > 0 in the following LMI
 
AX + X AT + BY + YT BT E X CT + YT DTu
 ET  2 I DTw  < 0 (5.31)
C X + Du Y Dw I
This inequality is known as the bounded real lemma. Provided that (5.31) is
feasible, it is guaranteed that the H gain is below , that is kHwz (s)k < ,
when applying the state-feedback matrix
K = YX1
 (5.32)
A prerequisite for this method to complete successfully is that the pair (A, B)
is controllable and (A, C ) is observable. This means that for the augmented
integral control to work, the augmented state  must be reected in the output
matrix C . The H performance is convenient to enforce robustness to model
uncertainty, and it is a direct measure for the L2 gain from disturbance w L2 to
the respective output z L2 .

H2 -Control

Let Hwz2 (s) denote the closed loop realization between w and z2 . Then, kHwz2 (s)k2 <
if there exist X2 = XT2 < 0, Y  Rm×n , and Z  Rn×n2 such that the following
LMIs are feasible
· ¸
AX2 + X2 AT + BY + YT BT E
<0
ET · I ¸
X2 X2 C2 T
(5.33)
<0
C2 X2 Z
Trace(Z) < 2
Consequently, the state-feedback u = Kx, where K is dened similarly to (5.32)
as K = YX1 2 , guarantees that the H2 -gain from w to z is below .

Pole Clustering

Assigning closed loop poles of a linear system can be seen as a tool for specifying a
minimum decay rate. This technique may also be used to ensure a minimum closed
5.3 Commutating Control 73

loop damping factor or a maximum bandwidth in order to avoid fast dynamics and
high frequency gain in the controller.
In order to describe convex regions in the complex plane in which the poles are
supposed to be put, we follow the denitions of Chilali and Gahinet (1996):

Denition 5.1 A subset D of the complex plane C is called an LMI region if there
exist a symmetric matrix   RnD ×nD and a matrix   RnD ×nD such that
D = {z  C : FD < 0} (5.34)
with
FD (z) =  + z + z̄ T < 0 (5.35)
where z̄ is the complex conjugate of z.

An LMI region is convex and symmetric about the real axis. Furthermore, LMI
regions are invariant under set intersection: The intersection of two LMI regions
D1 and D2 is also an LMI region with the characteristic function
FD1 D2 (z) = Diag(FD1 (z), FD2 (z)) (5.36)
As a consequence, an arbitrary region consisting of the intersection of conic curves,
vertical strips, and/or horizontal strips can be expressed in terms of LMI regions.

Denition 5.2 A matrix A is said to be D-stable if all its eigenvalues lie in D.

We may now summarize this for control synthesis purposes by the following theo-
rem (Chilali and Gahinet 1996):

Theorem 5.2 Let   RnD ×nD be a real symmetric matrix and   RnD ×nD .
Then Acl = A  BK has all its eigenvalues in the LMI region (5.35) if and only
if a real, symmetric, positive denite X  Rn×n and a real Y  Rm×n exist such
that the LMI
  X +   V +  T  VT < 0 (5.37)
where  is the Kronecker product and
V = AX + BY (5.38)
K = YX1 (5.39)
is feasible.

Observe that this is a generalization of Lyapunov stability for linear systems, that
is the eigenvalues of Ac having negative real part (Hurwitz), because the left half
complex plane is described by  = 0 and = 1. Then, Theorem 5.2 simply states
that
Ac X + XATc < 0 (5.40)

For more examples of various LMI regions and how to construct them, please refer
to Chilali and Gahinet (1996).
74 Controller Design

Stability For a Prescribed rmax

Let us now recast the results from Corollary 5.1 into a set of LMIs that can be
included with other LMIs in a state feedback synthesis procedure. For the sake of
clarity, this result is formulated as a Theorem.

Theorem 5.3 Suppose the system (5.4)-(5.5) is to track a suciently smooth

trajectory  d (t) = [xd (t), yd (t), d (t)]T by using the controller (5.6)-(5.7). If
there for a specied rmax > 0 exist matrices X = XT > 0 and Y  Rm×n such
that
³ ´
AX + XAT + BY + YT BT + rmax ST X + XSTT < 0
³ ´ (5.41)
AX + XA + BY + Y B  rmax ST X + XST
T T T T
< 0

is feasible, then the state feedback gain K = YX1 guarantees uniform (global)
exponential stability of xe .

Proof. Revisiting Theorem 5.1 we remember that the exponential stability of xe

is equivalent to the existence of a P = PT > 0 such that for all r(t) < rmax t 0
¡ ¢
(A  BK)T P + P (A  BK) + r PST + STT P < 0 (5.42)

Pre- and post-multiplying with P1 and substituting Y = KP1 we get

³ ´
AX + BY + XAT + YT BT + r ST X + XSTT < 0 (5.43)

Negative deniteness is ensured by (5.41) since this inequality is convex in r.

As for the other LMI criteria above, in a control synthesis procedure we end up
searching for a X = XT > 0 and a Y. Consequently, (5.41) can be combined
with any suitable LMI design criterion, that being H , H2 , pole clustering or
combinations of those, to obtain globally exponentially stable tracking controllers
for all rmax .

5.4 Acceleration Feedback

Having established the PID control framework we are now ready to expand the
state feedback control synthesis above with additional acceleration feedback. The
proposed controller and the certainty equivalence realization are particular contri-
butions of this thesis.

5.4.1 Dynamic Acceleration Feedback

First, let us consider the velocity and acceleration dynamics. Suppose that mea-
sured acceleration ya  Rnya where 1  nya  3 is given in the b-frame and let
5.4 Acceleration Feedback 75

  Rnya ×3 be the projection extracting those available accelerations ya = .


Let the applied control  be
³ ³ ´´
a f = Af af + Bf ya    e ST RT ( e ) d + RT ( e )
d
(5.44)
 =  PID  Ka af
where af  Rna is the ltered acceleration error described by the Af  Rna ×na
and Bf  Rna ×nya matrices.  PID is to be determined shortly. The velocity error
dynamics can be written

M e = D e +  PID  Ka af + w (5.45)

³ ´
At low frequencies af  ya    e ST RT ( e ) d + RT ( e ) d =  e which
means that the acceleration feedback term Ka af can be regarded as a change in
the system’s mass. Taking the position error dynamics into account, we get
 e = R( ) e
(5.46)
(M + Ka )  e = D e +  PID + w

Next, we want to nd a PID-like control law  PID for the system (5.46). As in
Section 5.3 integral action is obtained by integrating the position deviation and
assign gains Ki  R3×3 , Kp  R3×3 , Kd  R3×3 . Applying the very same PID-
controller as in (5.16)-(5.17), that is

 = e
(5.47)
 PID = Ki RT ( )  Kp RT ( ) e  Kd  e +  r

where the reference feed-forward  r is dened by (5.18). Collecting the states

into xe = [T ,  Te ,  Te , aTf ]T  R9+na , we can by dening the block-diagonal T :
R  R(9+na )×(9+na ) as follows
T() = Diag(RT (), RT (), I3+na ) (5.48)
1
and letting BM = Bf M express the complete error-dynamics on the compact
form
x e = TT ( )Ac T( )xe + Ew (5.49)
where
 
0 I 0 0
0 0 I 0
Ac =  (5.50)
M1 Ki M1 Kp M1 (D + Kd ) M1 Ka 
BM Ki BM Kp BM (D + Kd ) Af  BM Ka
£ ¤T
E= 0 0 I BTM (5.51)
The system parameters and controller gains have been isolated in the matrix Ac ,
and due to controllability the eigenvalues of Ac are freely assignable.
Theorem 5.1 and Corollary 5.1 can now be employed to establish stability prop-
erties for (5.49).
76 Controller Design

5.4.2 Controller Tuning

One of the main objectives of this thesis was to illustrate that constructive use
of measured acceleration could improve performance of DP systems compared to
PID or PD designs irrespective of their design philosophy. A reasonable way to do
this is to assign more or less identical poles to the resulting closed-loop systems.
The method is described in this section.
Step 1: Acceleration Feedback.
A rst-order low-pass lter was used to remove high-frequency noise components
from the two acceleration signals available, surge- and sway-acceleration. Thus,
we let Ka  R3×2 and

Af = Bf = diag(1/Tf , 1/Tf ) (5.52)

· ¸
1 0 0
 = (5.53)
0 1 0

where the lter constant Tf was selected so small

μ ¶
m11 m22
Tf ¿ min , (5.54)
d11 d22

that the acceleration feedback term could be regarded as a direct manipulation of

the mass (5.46), that is Ma = M + Ka .
Step 2: PID-control.
When the rotations are disregarded, we are left with nding a state-feedback
control u = Kx to the system

x e = Axe + Bu (5.55)

where    
0 I 0 0
A= 0 0 I  , B= 0  (5.56)
1
0 0 Ma D M1
a

When the gain matrix K  R3×9 is partitioned as follows

£ ¤
K = Ki Kp Kd (5.57)

the closed loop dynamics can be written as, remember that   R3 is the position
vector, Z t
¨ + 2 + 2  + K̄i
 (s)ds = w (5.58)
0
which we recognize as a second order dynamic system with integral action.

2 = M1
a (D + Kd ) (5.59)
2 = M1
a Kp (5.60)
K̄i = M1
a Ki (5.61)
5.4 Acceleration Feedback 77

If K̄i = 0 the matrices ,   R3×3 determine the natural frequency and rela-
tive damping respectively. In this study we de-coupled the individual degrees of
freedom by the following selection of gains

 = diag(1 , 2 ,  3 ) (5.62)
 = diag( 1 , 2 , 3 ) (5.63)
K̄i = diag(ki1 , ki2 , ki3 ) (5.64)

Consequently, for a constant selection ,  and K̄i , acceleration feedback control

as dened in Step 1 does not inuence the system’s tracking capabilities.

5.4.3 Output Feedback

This section discusses the extension of the state-feedback controller (5.44), (5.47)
to output-feedback by using an observer with wave ltering capabilities. Global
asymptotic stability of the complete system is established using Theorem A.2.
This analysis is valid for any commutating design following the algorithm outlined
in Section 5.3. For a more specic treatment of a simpler PID-controller please
see Lindegaard and Fossen (2003).

Observer review

The main objective of the implemented observer is to reconstruct the system’s LF

states and accelerations
ˆT , ˆ
T
T , 
xo = [ˆ  ]T (5.65)
based on measured positions and other types of sensor data, that is in particu-
lar (partial) velocity and acceleration feedback. Note the dierence between ˆ ,

the estimated acceleration, and ˆ  , the dierential equation used to update the
estimated velocity ˆ.
In Lindegaard et al. (2002) an observer for low-speed ship applications was pro-
posed and discussed. Based on the same stability arguments as the state-feedback
controller in Section four, that is for bounded yaw rate, the observer error

x̃o = xo  [ T ,  T ,  T ]T (5.66)

was shown to converge exponentially to zero

|x̃o (t)|  k |x̃o (0)| exp(t) , t 0 (5.67)

for some k,  > 0.

Observer-Feedback Control

Based on the principle of certainty equivalence, an observer-feedback controller is

realized substituting the actual states in the state-feedback controller (5.44) and
78 Controller Design

(5.47) by their estimated values  ˆ , and ˆ

ˆ,   . A new set of error variables following
the denitions (5.11)-(5.13) are thus


ˆe = 
ˆ  d (5.68)

ˆe ˆ  R ( ˆ e ) d
=  T
(5.69)
³ ´
ˆ
 e = ˆ
  r̂e ST RT ( ˆ e ) d + RT ( ˆ e ) d
³ ´
= ˆ
  (r̂  rd ) ST RT ( ˆ  d ) d + RT ( ˆ  d )
d (5.70)

The proposed observer-feedback controller with dynamic acceleration feedback for

(5.14)-(5.15) is

 = 
ˆe
a f = Af af + Bf ˆ  e (5.71)
ˆ = Ki RT ( ˆ )  Kp RT ( ˆ )ˆ
 e  Kd 
ˆ e  Ka af + ˆ r

where the reference feed-forward is given by

³ ´
ˆ r = DL RT ( ˆ e ) d + M (r̂  rd ) ST RT ( ˆ e ) d + RT ( ˆ e ) d (5.72)

Notice that in the rotations we have substituted with ˆ and consequently

state estimates appear non-anely in the control ˆ . Due to the linear bound
of kR( )  Ik and the inherent linear characteristics of the system, using (5.71)
instead of (5.44) and (5.47) does not compromise the asymptotic stability estab-
lished for state-feedback control. This can be summarized as follows

Theorem 5.4 Consider the system (5.4)-(5.5) for which there exists an ob-
server whose errors x̃o (t) converge asymptotically to zero (5.67). The observer-
feedback control (5.7)-(5.72) will guarantee UGAS of xe = [T ,  Te ,  Te , aTf ]T pro-
vided that

. The gains Ki , Kp , Kd , Ka are selected such that Ac as dened by (5.50) is

Hurwitz.

2. Maximum yaw rate rmax calculated by Corollary 5. does not exceed the phys-
ical bound, |r(t)|  rmax .

The formal proof is given in Appendix B.1.

5.5 Nonlinear Control

The main objective of this section is to identify the lacking terms in the commu-
tation based controller that if included would have guaranteed global exponential
5.5 Nonlinear Control 79

stability of the tracking error. A second objective is to show the similarities be-
tween controllers formulated in a reference parallel frame and our two-frame based
design. Using integrator backstepping (Krstić et al. 1995), we demonstrate that,
under some restrictions, reference parallel control (Strand 1999) is equal to body
xed control (Berge 1999). Furthermore, the derived tracking controller

• Considers Coriolis, centripetal, and nonlinear damping terms, while it does

not directly cancel any of these forces.

• Handles smooth, time-varying trajectories.

• Implements “true” integral action in the sense that the position deviation
updates the controller’s integral action.

5.5.1 Model Description

Consider a 3 DOF model with Coriolis C() = CT () and nonlinear damping
D() = DL + DN ()

 = R( )
(5.73)
M = C()  D() +  + w
A smooth reference trajectory is assumed given in the reference parallel d-frame
according to (5.8)-(5.10).

5.5.2 State Feedback Backstepping Control

Backstepping is a constructive design procedure for the control of nonlinear sys-

tems on feedback form. The original state variables are transformed into a new
set of variables, the errors z, for which a stabilizing controller is derived stepwise
together with a block diagonal Lyapunov function. There is, however, no dominat-
ing general method for obtaining controller integral action, yet this can be imple-
mented by introducing an additional ”step” (Aarset et al. 1998, Strand 1999) or
parameter adaptation (Godhavn et al. 1997, Berge 1999, Fossen et al. 2001). The
rst alternative complicates the derived controller (backstepping controllers tend
to be comparatively complicated and an additional step increases the complexity
even more). In a parameter adaptation setting, the ”integrator” will be updated
by a combination of the states. For example, in mechanical systems this means
that the integral action is updated by a sum of the position and velocity errors
in contrast to conventional linear control where only position error is used. The
advantage of such designs is that the integral gain can be chosen arbitrarily large
due to a relative degree of one between the controller and the constructed output
updating this integrator. On the other hand, the close connection between a linear
controller and a backstepping controller becomes less obvious.
We suggest yet another method for obtaining integral action in the tracking of
ships: The idea is to augment the position error with an extra integrator such that
80 Controller Design

two steps will be sucient. As a result, ”true” integral action, in the sense of it is
being updated by the position error alone, is achieved At the same time controller
complexity is reduced, at least when expressed in the z-variables, compared to a
three step method.
The objective is to follow the suciently smooth reference trajectory  d ,  d ,  d :
R0  R3 where  d = R( d ) d . Let the position error e be given in the d-frame

e = RT ( d ) (   d ) = RT ( d ) e (5.74)

such that

e =   d SRT ( d ) e + RT ( d ) (R( )  R( d ) d )
=   d Se + R( e ) e (5.75)

where the heading error is e =  d . A kind of integral action performed in

the Earth-xed frame can be augmented as follows:

 =  + CT11 R( d )e (5.76)

The matrix C11  R3×n is a projection used to isolate the components of e that
is subject to integral action.   Rn ×n should be diagonal and contain the
inverse of some large time constants. The larger these constants are the closer we
get to true integral action in the sense that  becomes an open integrator. For the
purpose of a simpler analysis we let  be non-zero as this allows us to establish
exponential stability more easily.
£ ¤
Theorem 5.5 below conrms that the tracking error eT ,  Te where  e =  
RT ( e ) d is UGES. It also provides conditions on how the individual gain matri-
ces should be selected such that reference parallel (Strand 1999) control becomes
identical to body-xed control (Berge 1999).

Theorem 5.5 Assume that there exist symmetric and positive denite 1 
Rn ×n , 2  R3×3 where 2 commutes with R( ) and C11 1 is selected such
that
R( )C11 1 = C11 1 CT11 R( )C11 (5.77)
Applying the controller

 =  + CT11 R( d )e
(5.78)
 = Ki ()RT ( )  Kp (, )RT ( e )e  Kd  e +  r

where the gains are

³ ´
Ki () = C() + D() + C2 + rMST 1
2 C11 1 C11
T

M12 C11 1 C11

T
(5.79)
Kp (, ) = 2 + (C() + D() + C2 ) C12
+M12 R ( )C11 1 C11 R( ) + rMC12 S
T T T
(5.80)
Kd = C2 + MC12 (5.81)
5.5 Nonlinear Control 81

and the reference feed-forward is

d ¡ T ¢
 r = (C() + D()) RT ( e ) d +MR ( e ) d (5.82)
dt
h iT
yields uniform global exponential stability of xe = T ,  Te ,  Te = 0 whenever
C12  R3×3 , C2  R3×3 , and   Rn ×n are selected such that

1  + T 1 > 0 (5.83)
2 C12 + CT12 2 > 0 (5.84)
D() + DT () + C2 + CT2 > 0 (5.85)

The proof is given in Appendix B.2.

If all positions are subject to integral action C11 = I, the gains are reduced to
1 1
Ki () = (C() + D() + C2 ) 1
2 1  M2 1  + rM2 1 S (5.86)
T

Kp () = 2 + (C() + D() + C2 ) C12 + M1

2 1 + rMC12 S
T
(5.87)
Kd = C2 + MC12 (5.88)

Considering now the low speed model (5.15) and C11 = I, the following gains
should be used
1 1 T
Ki (r) = (DL + C2 ) 1
2 1  M2 1  + rM2 S 1 (5.89)
1
Kp (r) = 2 + (DL + C2 ) C12 + M2 1 + rMC12 ST (5.90)
Kd = C2 + MC12 (5.91)

and the same reference feed-forward as for the commutating designs (5.18). From
(5.89)-(5.90) the missing yaw rate dependent terms needed for establishing UGES
are clearly visible. Since RT ( e )e = RT ( ) e and R( d )e =  e , we immediately
recognize that (5.78) is the very same controller as (5.17) apart from that the gains
are selected according to a dierent strategy.
It is also worth emphasizing that the yaw rate dependendent gains, those are
M1 2 1 S and MC12 S in (5.89) and (5.90) respectively, can be kept small
T T

while simultaneoulsy increasing the sum of the constant terms in Kp and Kd . In

other words, in accordance with the commutating controller (5.24), this leads to
higher bounds of rmax if the yaw rate components in Ki (r) and Kp (r) were to be
neglected.

5.5.3 Output Feedback

Even though the state feedback controller derived by backstepping guarantees

global exponential stability of the tracking error, it cannot be used directly in
a certainty equivalence kind of output feedback control setting. In fact, the re-
striction experienced with the commutation based designs, the bounded yaw rate
82 Controller Design

|r(t)|  rmax , is indeed relaxed, but the price paid is signicant: The error in the
compass heading estimate ˜ in combination with the nonlinear terms of the pro-
portional and integral gains excludes employing the results on cascaded systems
in establishing asymptotic stability of the closed loop system. The interconnec-
tion term between the observer errors 2 and the tracking errors 1 is no longer
linearly bounded in the states of 1 .
In Aarset et al. (1998) the authors used observer backstepping to avoid using
˜ = 0, but this contradicts the separation principle in the sense that the resulting
controller tuning depends on the observer. Another approach is pursued for more
general Euler-Lagrange systems (Loría and Panteley 1999, Aamo et al. 2001) where
is assumed measured and available for feedback. Neglecting the linearly wave
induced motion and thus using y (the measured heading angle) in the derived
controller (5.78), this procedure would be globally asymptotically stable here as
well. In the following we are, however, going to analyze a certainty equivalence
design to point out the ˜ -dependency and then suppose ˜ = 0 to eliminate the
conicting terms.

Prerequisites

First we summarize some of the properties of the Coriolis and damping matrices
needed in the forthcoming analysis.

• The Coriolis and centripetal matrix is linear in its argument, that is for
a,b  R3
C(a + b) = C(a) + C(b) (5.92)
As a consequence, it is bounded linearly in the argument as well: There exist
positive scalars cm ,cM > 0 such that
cm |a|  kC(a)k  cM |a| (5.93)

• The damping is linear plus quadratic throughout the entire velocity range
D() = DL + DN (). Then, for any a  R3 there are bounds dm and dM
such that
dm |a|  kDN (a)k  dM |a| (5.94)
Assume furthermore that the error in DN dened as
Derr (a, b) = DN (a  b)  DN (a) (5.95)
is bounded linearly in b. More specically for some derr > 0 the following
holds
kDerr (a, b)k  derr |b| (5.96)

Remark 1 In the scalar case it is known that

|a  b|  |a|  |b| a, b  R
which means that (5.96) at least covers any diagonal DN .
5.5 Nonlinear Control 83

Dene

N() = C() + D() (5.97)

Nerr (, 
˜) = Derr (, 
˜ )  C(˜
) (5.98)

then
 ) = N() + Nerr (, 
N(ˆ ˜) (5.99)
From the assumptions above we see that

|Nerr (, 
˜)|  (derr + cM ) |˜
| (5.100)

Observe when disregarding Coriolis and quadratic damping that N() = DL and
thus Nerr (, 
˜ ) = 0.

Observer-Feedback Control

The proposed observer feedback controller is obtained by substituting the state

variables with their estimated counterparts (certainty of equivalence), and without
loss of generality, integral action is applied in all three DOFs.
Let position and velocity error be dened by (5.68) and (5.69) respectively. The
position error decomposed in the reference parallel frame is thus

ê = RT ( e
d )ˆ = e  RT ( 
d )˜ (5.101)

The controller reads

 =  + R( d )ê
(5.102)
ˆ  )RT ( ˆ )  Kp (ˆ
= Ki (ˆ  )RT ( ˆ e )ê  Kd 
ˆe + ˆ r

where the reference feed-forward is

³ ´
 )RT ( ˆ e ) d + M r̂e ST RT ( ˆ e ) d + RT ( ˆ e ) d
ˆ r = N(ˆ (5.103)

and the gains are

1 1
Ki (ˆ  ) + C2 ) 1
 ) = (N(ˆ 2 1  M2 1  + r̂M2 1 S (5.104)
T

1
Kp (ˆ
 ) = 2 + (N(ˆ ) + C2 ) C12 + M2 1 + r̂MC12 ST
(5.105)
Kd = C2 + MC12 (5.106)

Theorem 5.6 Consider a suciently smooth reference trajectory d ,  d ,  d : R0 

R3 where  d = R( d ) d , and assume there exists an observer with asymptoti-
cally converging states xo . The controller (5.02)-(5.06) guarantees UGAS of the
£ ¤T
tracking error eT ,  Te provided that ˜ = 0 and the gains are selected according
to Theorem 5.5.

The proof is given in Appendix B.3.

84 Controller Design

Remark 2 Admittedly, by imposing bounds on the velocity vector  and the es-
timated heading error ˜ it would be possible to establish asymptotic stability by a
completion of the squares in the derivative of the Lyapunov function derived under
state-feedback. This, however, contradicts the objective of obtaining UGAS under
observer feedback. Instead we settled for ˜ = 0 ( is perfectly measured).

Separating the linear control terms from the nonlinear ones can be done by intro-
ducing the following denitions:
1
Gi = C2 1
2 1  M2 1  (5.107)
Xi1 = C2 1
2 1 (5.108)
Xi2 = M12 1 S
T
(5.109)
Gp = 2 + C2 C12 + M1
2 1 (5.110)
Xp = MC12 ST (5.111)

The gains are thus given as

 ) = Gi + N(ˆ
Ki (ˆ  )Xi1 + r̂Xi2 (5.112)
Kp (ˆ
 ) = Gp + N(ˆ
 )C12 + r̂Xp (5.113)
Kd = C2 + MC12 (5.114)

where the Gi  R3×3 and Gp  R3×3 are the constant (linear control) gains. The
terms requiring ˜ = 0 are the nonlinear factors of Ki (ˆ  ) and Kp (ˆ
 ). It should
be noted that even though ˆ r indeed contains N(ˆ ), for this term alone it is not
necessary to assume ˜ = 0 because | d | is known to be bounded.