Topic Dynamic Positoning

Topics in dynamic positioning
System identification, GNSS and MRU lever

arm estimation, and hybrid integral action
Svenn Are Tutturen
Marine Technology
Submission date: August 2014
Supervisor:
Roger Skjetne, IMT
Co-supervisor:
ivind Kre Kjerstad, IMT
Norwegian University of Science and Technology

Department of Marine Technology
NTNU Trondheim
PROJECT DESCRIPTION SHEET

Name of the candidate:
Svenn Are Tutturen
Field of study:
Marine control engineering
Thesis title (Norwegian):
Tema innen dynamisk posisjonering: Systemidentifikasjon, estimering

av GNSS og MRU mlearmer, og hybrid integralregulering.
Thesis title (English):
Topics in dynamic positioning: System identification, GNSS and MRU

lever arm estimation, and hybrid integral action.
Background
The low-speed dynamic model of a DP vessel is well known, both high fidelity process models and lowfidelity design models. However, the methods and procedures for identification of the model parameters
is still a field of development. For specialized vessels such parametric studies are traditionally
conducted through model-scale experiments in a model basin and through hydrodynamic computer
programs. For full-scale testing, on the other hand, recent technological development of sensor
technologies, data acquisition, and communication make more accurate dynamic data available for
analysis and estimation of vessel parameters.
The ship R/V Gunnerus, owned by NTNU, is planned in the spring of 2015 to be retrofitted with new
podded propulsion. The old propellers and rudders will then be replaced by new prototype podded
permanent magnet (PM) azimuth thrusters. As part of this replacement, a study on the efficiency of the
thruster system, low-speed responsiveness, and maneuverability of the ship is conducted before and
after the retrofit to document changes and improvements to the ships dynamic behavior. This study is
conducted through several sea trials for the ship before and after the retrofit.
The scope of this project is to study thruster dynamics and low-speed ship models for R/V Gunnerus,
and to analyze the corresponding sea trials data collected through the different tests on August 14-15,
2013. The objective here is to draw conclusions on methods for and results from system identification of
the thruster dynamics and low-speed dynamics of a DP vessel based on numerical studies for R/V
Gunnerus. In addition, the objective is to develop a method for estimation of the lever arms for the posref and IMU sensors in the vessel based on persistently exciting motions of the DP vessel.
Work description
1) Perform a literature review, providing relevant references, on:
Modeling of thruster dynamics and low-speed dynamics for DP vessels.
Methods for system identification, especially step tests and system identification software.
Write a list with abbreviations and definitions, and a section explaining particularly relevant terms
and concepts related to DP systems and system identification methods in an alphabetic order.
2) Propose relevant dynamic models for low-speed thruster and rudder dynamics, considering
particularly:
advanced model(s) (incl. the relationships between thrust force, RPM, and power),
simplified model(s) that are appropriate for system identification based on available
measurements, and
thruster configuration and net thruster force/moment (thrust allocation).
3) Propose relevant dynamic models for low-speed and wave motions of a DP vessel, considering
particularly:
high-fidelity models incl. environmental loads, and
simplified low-fidelity models (3 DOF DP; forward speed dynamics; steering dynamics)
appropriate for system identification based on available measurements.
NTNU
Faculty of Engineering Science and Technology

4) Analyze the results from Test 006 Thrust agility test and the settling times for the produced
thrust in each commanded thrust direction. Generate a polar plot with response time in the radial
direction and plot the corresponding envelopes for thrust response and speed response based on the
different tests carried out.
5) Analyze the results from Test 4-corner DP maneuver. Based on the KM log, report in a table the
duration for each of the 4 maneuvers, and the accumulated thrust force and power for each thruster
during each of the 4 maneuvers and in total. Repeat the analysis for the same tests conducted later
on 06.11.2013 and compare and discuss the results.
6) For each of the subtests conducting pure DP (zero speed) in the test reports of 15.08.2013 and
08.11.2013, make a table and report the following values (based on a 5 min. interval of the test):
Average wind; Observed waves; Ocean current; Positioning and heading accuracy (Std); Average
thruster force and accumulated power (energy).
7) Based on step test responses for Gunnerus, identify relevant thrust [] and propeller speed
[] time constants for the thrusters.
8) Based on step test responses in surge for Gunnerus, and assuming the thruster steps respond
instantaneously, identify relevant surge speed [/] time constants for the ship.
9) Based on step test responses in sway for Gunnerus, and assuming the thruster steps respond
instantaneously, identify relevant sway speed [/] time constants for the ship.
10) Based on step test responses in yaw for Gunnerus, and assuming the thruster steps respond
instantaneously, identify relevant yaw rate [/] time constants for the ship.
11) For the estimation of pos-ref lever arms, show the following:
The regressor matrix in an adaptive setup satisfies a PE requirement given some assumptions of
motion of the vessel. Conclude with a proposition.
The lever arms and the rotation point 0 are uniformly completely observable through the
measurements, given some assumptions of motion of the vessel. Conclude with a proposition.
The physical interpretation of the rotation point 0 .
Then do the following:
Propose an observer-based estimation design and conclude with a theorem.
Propose an alternative estimation design (adaptive, RLS, etc.) and conclude with a theorem.
Verify the methods using relevant sea-trial tests for Gunnerus as a case study.
Discuss the numerical results of the two estimation methods when compared to the real
measured values on the ship.
12) Propose a setup for also including the MRUs into the lever arm identification system:
Derive and discuss the observability properties of the system.
Propose a numerical estimation design based on measurements during vessel motion.
Verify the method using Gunnerus sea-trial as a case study and discuss the results.
Tentatively:
13) Derive a hybrid PID-type DP control law, including a discrete resetting of the integral action to
better compensate fast variations in the bias forces. Prove stability of the resulting closed-loop
hybrid dynamical system, and perform simulations to verify the effectiveness of the algorithm
compared to a conventional design.
NTNU
Faculty of Engineering Science and Technology

Guidelines
The scope of work may prove to be larger than initially anticipated. By the approval from the supervisor,
described topics may be deleted or reduced in extent without consequences with regard to grading.
The candidate shall present his personal contribution to the resolution of problems within the scope of work.
Theories and conclusions should be based on mathematical derivations and logic reasoning identifying the
various steps in the deduction.
The report shall be organized in a rational manner to give a clear exposition of results, assessments, and
conclusions. The text should be brief and to the point, with a clear language. The report shall be written in
English (preferably US) and contain the following elements: Abstract, acknowledgements, table of
contents, main body, conclusions with recommendations for further work, list of symbols and acronyms,
references, and (optionally) appendices. All figures, tables, and equations shall be numerated. The original
contribution of the candidate and material taken from other sources shall be clearly identified. Work from
other sources shall be properly acknowledged using quotations and a Harvard citation style (e.g. natbib
Latex package). The work is expected to be conducted in an honest and ethical manner, without any sort of
plagiarism and misconduct. Such practice is taken very seriously by the university and will have
consequences. NTNU can use the results freely in research and teaching by proper referencing, unless
otherwise agreed upon.
The thesis shall be submitted with two printed and electronic copies, to 1) the main supervisor and 2) the
external examiner, each copy signed by the candidate. The final revised version of this work description
must be included. The report must appear in a bound volume or a binder according to the NTNU standard
template. Computer code and a PDF version of the report shall be included electronically.
Start date:
1 February, 2014
Supervisor:
Co-advisor(s):
Roger Skjetne
ivind K. Kjerstad
Due date:
As specified by the administration.
Trondheim, __________________
_______________________________
Roger Skjetne
Supervisor
Summary
This thesis consists of three main parts. The first part of the master thesis looks at the
identification of thruster dynamics and low speed ship dynamics. The relevant parameters
identified are time constants and time delays in the system. Simple step tests are used for
the identification. Different models for identification are suggested, both for uncoupled
surge, sway, and yaw dynamics. Other test results, such as agility plots, DP 4 corner tests,
and pure DP tests (stationkeeping) are reported. All the results are to be compared to
similar tests performed after R/V Gunnerus has a retrofit of the thruster system.
The second part discusses another problem, and that is the topic of numerically estimating the body frame position of the GNSS and MRU sensors. For the GNSS position an
Luenberger observer design and an adaptive scheme are proposed and analyzed. The estimation designs are tested using numerical simulations and experimental data from the
Gunnerus sea trials. A similar Luenberger observer is proposed for the MRU positions,
and experimental data from the sea trials are used to test the observer.
The third part discusses a hybrid augmentation of integral action. The motivation is a DP
system, where typically the integral action is tuned very low to avoid oscillations due to
the integral action. When there is a sudden load change, such as a ice load that hits the
vessel, or if a mooring wire snaps, then a hybrid update augmentation could be useful, to
speed up the convergence of the integral action. The update law is a linear update law
based on the error in the states (the velocity for the DP system). The augmentation can
significantly improve performance, especially for very large disturbance changes.
Sammendrag
Denne avhandlingen bestar av tre hoveddeler. Den frste delen av masteravhandlingen ser
pa identifisering av thrusterdynamikk og lavhastighets skipsdynamikk. Relevante parametre som blir identifisert er tidskonstanter og tidsforsinkelser i systemet. Enkle sprangresponser er brukt i identifikasjonen. Forskjellige modeller for identifisering er foreslatt,
bade for jag-, svai- og gir-dynamikk. Ander testresultater, deriblant en grafisk fremstilling
av thrusterrespons, en DP-manver i lav hastighet, samt en DP-manver hvor skipet holder
seg i ro er rapportert. Alle resultatene vil bli sammenlignet med lignende tester etter FF
Gunnerus har fatt ombygd sitt thrustersystem.
Del to av avhandlingen omhandler en annen problemstilling, og det er a numerisk estimere
armene til GNSS- og MRU-sensorer. For GNSS-problemet sa er et Luenberger estimatordesign presentert, og en adaptiv lsning er ogsa foreslatt og analysert. De foreslatte
estimatorene er testet ved a bruke numeriske simuleringer som benytter eksperimentell
data fra sjprvene med FF Gunnerus. Et lignende Luenberger-oppsett er foreslatt for
MRU-armene, og eksperimentell data fra sjprvene er brukt til teste estimatoren.
Den tredje delen diskuterer en hybrid augmentering til integraleffekt i en kontroller. Motivasjonen er dynamisk posisjonering, hvor integraleffekten typisk er stilt inn ganske lavt
for a forhindre at integraleffekten induserer svingninger i systemet. Nar det intreffer en
plutselig lastforandring, for eksempel fra en islast, eller at en forankringsline ryker, sa kan
det vre nyttig med en hybrid augmentasjon av integraleffekten, for a ke hastigheten
til integraleffekten sin konvergens. Oppdateringsloven er proposjonal med strrelsen til
feilen i tilstandene i systemet (for det foreslatte DP-systemet brukes feilen i hastighet for a
bestemme denne strrelsen). Denne hybride augmenteringen kan ke ytelsen til systemet,
spesielt ved store og plutselige lastendringer.
ii
Acknowledgement
This thesis was completed between January and July 2014 alongside a course for the integrated PhD studies as the last part of my masters degree in Marine Technology at the
Norwegian University for Science and Technology (NTNU).
The two first parts of the thesis are a continuation of the project thesis from last semester.
This is a stand alone thesis, and the project thesis is not needed for the understanding of
this work. The third and last part of the Thesis is an extension of a paper written for the
PhD course MR8500 - Advanced Topics in Marine Control Systems.
I would like to thank the crew at R/V Gunnerus for two pleasant cruises (sea trials), and
MARINTEK for all the data collection. I would like to thank Seatex for the MRUs at R/V
Gunnerus, and allowing us to use the data from these, and I would like to thank Professor
Sverre Steen at the Department of Marine Technology, NTNU, that have led the sea trials.
I would especially like to thank my supervisor Professor Roger Skjetne at the Department
of Marine Technology, NTNU. He has been a great motivation for the study of control
theory, especially the nonlinear part. Roger has also helped guide me through the master
thesis, and I am grateful. I would also like to thank Professor Andrew Teel at UCSB that
taught MR8500 (at NTNU) in hybrid control theory, and he has been a great inspiration
for working on hybrid control. Lastly, I would like to thank my co-advisor PhD candidate
ivind K. Kjerstad at the Department of Marine Technology for helpful discussions.
iii
Table of Contents
Summary
Sammendrag
ii
Acknowledgements
iii
List of Tables
ix
List of Figures
xi
3
iv
Introduction
1.1 Background . . . . . . . . . . .
1.1.1 R/V Gunnerus sea trials
1.2 Main contributions . . . . . . .
1.3 Thesis structure . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
1
2
3
3
Modeling
2.1 Vessel models . . . . . . . . . . . . . .
2.1.1 High fidelity model . . . . . . .
Kinematics . . . . . . . . . . .
Kinetics . . . . . . . . . . . . .
2.1.2 Uncoupled surge dynamics . . .
2.1.3 Sway-yaw subsystem . . . . . .
First order approximation . . .
2.1.4 Low speed models . . . . . . .
DP (3 DOF) - Linearized model
2.2 Thruster dynamics . . . . . . . . . . .
2.2.1 Simplified model . . . . . . . .
2.2.2 Thrust allocation . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
5
5
6
7
7
8
8
8
9
10
11
Model Identification
.
.
.
.
.
.
.
.
.
.
.
.
13
3.1
3.2
3.3
3.4
3.5
3.6
4
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
System identification by step tests . . . . . . . . . . . . . . . . . . . . .
3.2.1 First order system . . . . . . . . . . . . . . . . . . . . . . . . . .
2nd order system . . . . . . . . . . . . . . . . . . . . . . . . . .
First order plus dead time approximation (FOPDT) . . . . . . . .
Approximating higher order systems with first order time delay
functions . . . . . . . . . . . . . . . . . . . . . . . . .
Proposed identification models for the step tests . . . . . . . . . . . . . .
3.3.1 Thrust and shaft speed time constants . . . . . . . . . . . . . . .
3.3.2 Speed time constants . . . . . . . . . . . . . . . . . . . . . . . .
Proposed model for surge step tests . . . . . . . . . . . . . . . .
Proposed model for the sway step responses . . . . . . . . . . . .
Proposed model for the yaw step direction test responses . . . . .
Uavailable data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Thrust data . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Thruster mapping for the main thrusters . . . . . . . . . . . . . .
Tunnel thruster . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Power data . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Identification results . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.1 Thrust agility tests . . . . . . . . . . . . . . . . . . . . . . . . .
Agility tests with no yawing moment applied . . . . . . . . . . .
Agility test with applied yawing moment . . . . . . . . . . . . .
3.5.2 DP 4 corner tests . . . . . . . . . . . . . . . . . . . . . . . . . .
DP 4 corner - calm sea . . . . . . . . . . . . . . . . . . . . . . .
DP 4 corner - at sea . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.3 Pure DP tests . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Results - calm sea . . . . . . . . . . . . . . . . . . . . . . . . . .
Results - at sea . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.4 Step tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
About how tests are read . . . . . . . . . . . . . . . . . . . . . .
Surge step tests . . . . . . . . . . . . . . . . . . . . . . . . . . .
Sway step tests . . . . . . . . . . . . . . . . . . . . . . . . . . .
Yaw step tests . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Yaw rate time constant . . . . . . . . . . . . . . . . . . .
Rudder time constant . . . . . . . . . . . . . . . . . . . .
Conclusions and further work . . . . . . . . . . . . . . . . . . . . . . . .
Lever arm estimation for GNSS and MRU sensors

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Scope . . . . . . . . . . . . . . . . . . . . . . .
4.2 GNSS . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Luenberger observer design . . . . . . . . . . .
Observability assessment . . . . . . . . . . . . .
Luenberger observer design and stability analysis
4.2.2 Adaptive observer design . . . . . . . . . . . . .
Adaptive observer design and stability analysis .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
13
13
13
14
15
15
15
16
16
17
17
17
17
18
18
19
19
19
19
20
22
23
23
24
27
27
29
33
34
34
37
41
43
44
45
46
46
47
48
48
49
52
53
54
v
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
55
57
58
59
60
61
61
61
62
63
63
65
67
67
68
69
71
72
73
73
Hybrid integral action for DP

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . .
5.1.2 Literature review . . . . . . . . . . . . . . . . . . . .
5.1.3 Scope . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Mathematical modelling . . . . . . . . . . . . . . . . . . . .
5.2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . .
About the approach . . . . . . . . . . . . . . . . . . .
5.2.2 First order linear systems . . . . . . . . . . . . . . . .
Flow set . . . . . . . . . . . . . . . . . . . . . . . . .
Flow map . . . . . . . . . . . . . . . . . . . . . . . .
Jump map . . . . . . . . . . . . . . . . . . . . . . . .
5.2.3 DP system . . . . . . . . . . . . . . . . . . . . . . .
Flow set . . . . . . . . . . . . . . . . . . . . . . . . .
Flow map . . . . . . . . . . . . . . . . . . . . . . . .
Jump map . . . . . . . . . . . . . . . . . . . . . . . .
5.3 Case studies . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 First order linear system . . . . . . . . . . . . . . . .
Results first order system - with knowledge of x2 . . .
Results first order system - without knowledge of x2 .
5.3.2 Example - DP . . . . . . . . . . . . . . . . . . . . . .
Results DP - with knowledge of x3 . . . . . . . . . .
Results DP - without knowledge of x3 - sampling of x2
5.4 Conclusions and further work . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
75
75
75
75
76
76
76
77
77
78
78
79
81
82
82
82
85
85
87
88
89
90
92
93
4.3
4.4
4.5
5
vi
Percistence of excitation . . . . . . . . . . .
MRU . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Luenberger observer design . . . . . . . . .
Observability assessment . . . . . . . . . . .
4.3.2 Observer equations and stability . . . . . . .
Case studies . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Case studies - GNSS (GPS) lever arms . . .
Gunnerus data . . . . . . . . . . . .
Installed GPS antennas . . . . . . . .
Case study 1 - Simulated GPS data . . . . .
Luenberger observer results . . . . .
Adaptive observer results . . . . . . .
Case study 2 - Experimental GPS data . . . .
Luenberger observer results . . . . .
Adaptive observer results . . . . . . .
Simulated vs recorded GPS data . . . . . . .
4.4.2 Case study: MRU lever arms - Gunnerus data
Installed MRUs . . . . . . . . . . . . . . .
Luenberger observer results . . . . . . . . .
Conclusions and further work . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Bibliography
94
Appendix
97
vii
List of Tables
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13
3.14
3.15
3.16
3.17
3.18
3.19
3.20
3.21
3.22
3.23
3.24
3.25
3.26
3.27
KM thrustmapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Back-calculated thrust mapping . . . . . . . . . . . . . . . . . . . . . .
Environmental conditions thrust agility, no yawing moment . . . . . . . .
Environmental conditions thrust agility, with yawing moment . . . . . . .
DP 4 corner - environmental conditions calm sea . . . . . . . . . . . . .
DP 4 corner - maneuver, calm sea . . . . . . . . . . . . . . . . . . . . .
DP 4 corner - environmental conditions at sea, test 1 . . . . . . . . . . .
DP 4 corner - maneuver, at sea - test 1 . . . . . . . . . . . . . . . . . . .
DP 4 corner - environmental conditions, at sea - test 2 . . . . . . . . . . .
DP 4 corner - maneuver, at sea - test 2 . . . . . . . . . . . . . . . . . . .
Pure DP, calm sea - test 1 . . . . . . . . . . . . . . . . . . . . . . . . . .
Pure DP, at sea - test 1 . . . . . . . . . . . . . . . . . . . . . . . . . . .
Surge step tests, forward speed - thrust and shaft speed time constants . .
Surge step tests, forward speed - speed time constants . . . . . . . . . . .
Surge step tests, backward speed - thrust and shaft speed time constants .
Surge step tests, backward speed - speed time constants . . . . . . . . . .
Sway starboard step tests - thrust and shaft speed time constants for sb
thruster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.28 Sway step tests, main propellers, starboard - speed time constants . . . . .
3.29 Sway port step tests - thrust and shaft speed time constants for port thruster
viii
18
18
20
22
24
24
25
26
26
27
28
28
29
29
30
30
31
31
32
32
33
33
36
36
37
37
39
39
41
3.30
3.31
3.32
3.33
3.34
3.35
3.36
Sway step tests, main propellers, port - time constant

Yaw step tests - time constant 0-50 . . . . . . . . . .
Yaw step tests - time constant 100-0 . . . . . . . . .
Yaw step tests - rudder time constant, starboard 3 kn .
Yaw step tests - rudder time constant, port 3 kn . . .
Yaw step tests - rudder time constant , starboard 6 kn
Yaw step tests - rudder time constant, port 6 kn . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
41
44
44
44
44
44
45
4.1
4.2
4.3
4.4
4.5
4.6
GPS coordinates . . . . . . . . . . . . . . . . . . . .
Results, lever arm - observer, simulated GPS data . .
Results, lever arm - adaptive, simulated GPS data . .
Results, lever arm - observer, experimental GPS data
Results, lever arm - adaptive, experimental GPS data
MRU coordinates . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
62
64
66
68
69
72
ix
List of Figures
1.1
Thruster retrofit, illustration . . . . . . . . . . . . . . . . . . . . . . . . .
2.1
Thruster configuration, Gunnerus - before retrofit . . . . . . . . . . . . .
12
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13
3.14
3.15
3.16
3.17
3.18
Agility plot, thrust time constant, no yawing moment - half thrust .

Agility plot, thrust time constant, no yawing moment - full thrust .
Agility plot - thrust time constant, with yawing moment . . . . . .
DP 4 corner, illustration . . . . . . . . . . . . . . . . . . . . . . .
NE-plot of DP 4 corner test 1, at sea . . . . . . . . . . . . . . . .
NE-plot of DP 4 corner test 2, at sea . . . . . . . . . . . . . . . .
Step test in surge - forward velocity, part 1 . . . . . . . . . . . . .
Step test in surge - forward velocity, part 2 . . . . . . . . . . . . .
Step test in surge - backward velocity, part 1 . . . . . . . . . . . .
Step test in surge - backward velocity, part 2 . . . . . . . . . . . .
Step test in sway - starboard (sway, surge, total speed) . . . . . . .
Step test in sway - starboard (total speed, shaft speed fb, power fb)
Step test in sway - port (sway, surge, total speed) . . . . . . . . .
Step test in sway - port (total speed, shaft speed fb, power fb) . . .
Step test in yaw - starboard at 3kn . . . . . . . . . . . . . . . . .
Step test in yaw - port at 3kn . . . . . . . . . . . . . . . . . . . .
Step test in yaw - starboard at 6kn . . . . . . . . . . . . . . . . .
Step test in yaw - port at 6kn . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
21
21
22
23
25
26
35
35
36
37
38
39
40
40
42
42
43
43
4.1
4.2
4.3
4.4
4.5
4.6
4.7
Plot of p, q, and r for GPS case study . . . . . . . . . .

Lever arm coordinates, l1 , observer, simulated GPS data
Lever arm coordinates, l2 , observer, simulated GPS data
NE plot of P0 , observer . . . . . . . . . . . . . . . . . .
Lever arm coordinates, l1 , adaptive, simulated GPS data
Lever arm coordinates, l2 , adaptive, simulated GPS data
NE plot of P0 , adaptive observer . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
62
64
64
65
66
66
67
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
4.8
4.9
4.10
4.11
4.12
4.13
Lever arm coordinates, l1 , observer, experimental GPS data . . . . .

Lever arm coordinates, l1 , adaptive, experimental GPS data . . . . .
Comparison between recorded and simulated GPS data, no added arm
Comparison between recorded and simulated GPS data, added arm . .
Plot of p, q, and r for MRU case study . . . . . . . . . . . . . . . . .
Lever arm coordinates MRU case study . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
68
69
70
71
72
73
5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
5.10
5.11
5.12
Position and integral state error first order system, no sampling .

Lyapunov function and integral state error zoom, no sampling .
Position and integral state error first order system, with sampling
Lyapunov function and integral state error zoom, with sampling
Position error DP, no sampling . . . . . . . . . . . . . . . . . .
Velocity error DP, no sampling . . . . . . . . . . . . . . . . . .
Integral error DP, no sampling . . . . . . . . . . . . . . . . . .
Lyapunov function and zoom on integral error, no sampling . . .
Position error DP, with sampling . . . . . . . . . . . . . . . . .
Velocity error DP, with sampling . . . . . . . . . . . . . . . . .
Integral error DP, with sampling . . . . . . . . . . . . . . . . .
Lyapunov function and zoom on integral error, with sampling . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
87
87
88
88
90
90
91
91
92
92
93
93
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
xi
Definitions
BODY frame
a reference frame fixed to the craft, where the x axis points forward
(from aft to fore), the y axis to starboard, and the z axis points down (Fossen, 2011)
Dynamically positioned vessel
a free-floating vessel which maintains its position
(fixed location or predetermined track) exclusively by means of thrusters (Fossen, 2011)
Maneuvering
2011)
the study of ship motion in the absence of wave excitation (Fossen,
NED frame
a reference frame with x axis to pointing to the true North, y axis pointing
East, and z axis pointing down, normal to the Earth. Considered interial for a marine vessel operating in a local area (Fossen, 2011)
Seakeeping
the study of motion of a marine craft on constant course and speed when
there is wave excitation (Fossen, 2011)
List of acronyms
DOF
Degree of freedom
DP
fb
Dynamic positioning
feedback
FOPDT
GNSS
First order plus dead time

Global navigation satellite system
GPS
Global positioning system
INS
Inertial navigation system
LTV
Linear time varying
xii
MRU
N
Motion reference unit

North
NE
Northeast
NED
North, east, down
NW
Northwest
PE
Persistence of excitation
PID
Proportional, integral, derivative
rpm
revolutions per minute
rps
S
revolutions per second

South
sb
starboard
SE
Southeast
SW
Southwest
UGAS
UGpAS
Uniformly globally asymptotically stable

Uniformly globally pre-asymptotically stable
xiii
Nomenclature
orientation vector in NED-coordinates
body fixed velocity vector
body fixed angular velocity vector
vector of Euler angles
state vector
output (response)
time delay
time constant
Notation
Lower case bold letters, such as x indicate a vector, and uppercase bold letters such as A
represent matrices.
Time derivatives are indicated by the dot notation, such that x is the time derivative of
x. For discrete updates the notation x+ is used to indicate the value of x after a discrete
update.
xiv
Chapter
Introduction
1.1
Background
Dynamic positioning (DP) of marine vessels are a well-established commercial application

which enables automatic position and heading control in various operational scenarios (see
for instance (Kongsberg Maritime, 2006) and (Rolls-Royce Marine, 2009)). It is also an
established and active field of research (Fossen, 2011). Therefore, developments in control
theory and sensor technology that may improve performance and/or reduce the economical
effort of designing, verifying, and operating such systems are of interest.
This thesis considers three topics which may have such implications:
Identification of low-speed vessel models and analysis of actuation system performance,
estimation of GNSS and MRU lever arms based on vessel maneuvering data,
hybrid integral action for discontinuous updates of the control systems environmental disturbance rejection.
Although each of these topics will be given a more in depth introduction in their respective
chapters, the two first topics apply an extensive data set gathered during two sea trial
campaigns using the NTNU owned research vessel R/V Gunnerus.
The reason for these sea trials is that R/V Gunnerus has a planned retrofit of the thruster
system. Figure 1.1 below show a sketch of the old and new thruster configuration. The
illustration show the current thruster configuration with two propellers and rudders in the
stern. These will be replaced with two prototype permanent magnet azimuth thrusters.
There is also a tunnel thruster in the bow, that will not be replaced.
1
Chapter 1. Introduction
Figure 1.1: Thruster retrofit, illustration - thrust region before (left), and after (right) the retrofit
1.1.1
R/V Gunnerus sea trials
The first sea trial was performed in calm sea in the Trondheimsfjord in August 2013, and
the second one in open sea outside Kristiansund in November 2013. The experimental test
campaigns ware performed by NTNU and MARINTEK. The tests relevant for this master
thesis were performed in order to compare the response and performance of the current
thruster configuration of R/V Gunnerus to the retrofit. This include the performance of
the thruster system itself, and the performance in DP mode. When the retrofit is finished,
there will be more sea trials.
The tests performed relevant for the thesis include several step tests in surge, sway, and
yaw, where the thruster response, and the combined thruster and vessel hull response were
of interest. Thrust agility tests were performed. This results in a plot showing the time
constants of the of the thruster responses in different directions. A step in thrust force is
applied for heading of 0 , 45 , 90 , 135 , 180 , 135 , 90 , and 45 .
There were performed several stationkeeping tests, referred to as pure DP in this thesis.
The reported values include the average thrust force, the power consumption, and the
2
1.2 Main contributions
position accuracy. Another maneuver, called the DP 4 corner, which is more thoroughly
explained in Section 3.5.2 was performed. This maneuver tests the performance of the DP
system in low speed. The accumulated thrust force, and the accumulated power is logged,
along with the duration of the maneuver(s).
All the above tests are meant to give a good platform for comparison when evaluating the
performance of the retrofit compared to the old configuration. The tests are all explained
more thoroughly in Section 3.5, before the results are presented.
Also, experimental data from some of the performed maneuvers from these sea trials will
be used in the lever arm estimation problem.
In the thesis the first sea trial in August is referenced as either the sea trial in calm sea
(or calm water), the first sea trial, or the sea trial performed in August. Similarly is
the sea trial from November either references as the sea trial at sea (or open sea), the
second sea trial, or the November sea trial.
During the two sea trials MARINTEK has logged data, and for the first sea trial the log
from the Dynamic Positioning (DP) system was collected. This log was collected by a
Kongsberg Maritime representative, and sent to the MARINTEK log. This log is only
available from the first sea trial. The data collected exclusively by MARINTEK is referenced as MARINTEK data, or data collected by MARINTEK, and the data from the
DP system is referenced as the DP data.
1.2
Main contributions
For the system identification part of Chapter 3 simple models for identification are proposed, and the response of thruster system, and the combined response of the thruster
system and the vessel are reported. Data from several DP maneuvers are reported, such
that the data can be compared to R/V Gunnerus data after retrofit.
In Chapter 4 the main contribution are algorithms that estimate several GNSS (or MRU)
sensor lever arms, by using a simple setup without including the INS, and for the GNSS
case the algorithms use the vessel velocity, and the raw data from GNSS anteannas only.
In Chapter 5 the main contribution are the rigorous stability proofs of the hybrid system
for both the first order linear system, and the DP system. In this thesis limits for how large
these discrete updates can be has been found, and it is proven to be stable.
1.3
Thesis structure
In Chapter 2 modeling relevant for the other chapters is covered. This includes modeling
of vessel kinetics and kinematics, and thruster dynamics. High fidelity models are covered,
and also the decoupled, and simplified models needed for system identification in Chapter
3.
3
Chapter 1. Introduction
In Chapter 3 system identification by step tests are discussed, before suggested models
used for the identification are proposed. In Section 3.5 the step tests are analyzed, and
other test results are presented, such as the thrust agility tests, DP 4 corner, and pure DP
tests (stationkeeping).
Chapter 4 discuss and analyze the lever arm estimation problem, and results on experimental data are presented. In Chapter 5 the hybrid integral action is presented, and analyzed,
before simulation results are shown on a first order linear system, and a DP system.
Chapter 3 - 5 each discuss a separate topic, and each chapter has a conclusion including
suggestions for further work within the respective chapters. This is because each chapter
presents results and they are separate in their topic.
Chapter
Modeling
2.1
Vessel models
The modeling of following section is to a high degree based on material from Fossen
(2011).
Vessel kinetics is normally expressed in BODY-frame, and the kinematics in North-EastDown (NED)-frame. There exists numerous models with varying complexity, dependent
on operation. The main difference is for high- and low-speed applications, where different
forms of hydrodynamic damping dominate for different speed regimes, and the Coriolis
and centripetal terms become negligible for low speed.
2.1.1
High fidelity model
Kinematics
Vessel kinematics is given as (Fossen, 2011)
= R(),
(2.1)
where

= N
>
R61 ,
(2.2)
contain the North, East, and down positions, and the angular orientation (Euler angles) in
roll (), pitch (), and yaw (). The velocity vector is given as

= u
>
R61 ,
(2.3)
5
Chapter 2. Modeling
where u, v, and w are the velocities in surge (x-direction), sway (y-direction), and heave
(down-direction), respectively. The other half of is given by the angular velocities in roll
(p), pitch (q), and yaw (r). The vector is given in NED-coordinates, and in BODYcoordinates. The rotation matrix R() transforms the coordinate frame between BODY
to NED and is given as
cc sc + css ss + ccs
R() = sc cc + sss cs + ssc R33 ,
(2.4)
s
cs
cc
where c() = cos(), s() = sin().
Also, the time derivative of the rotation matrix is given by
R()
= R()S()
where S() is a skew-symmetric matrix given by
0 r
0
S() = r
q p
q
p .
0
(2.5)
(2.6)
Kinetics
The vessel kinetics is generally given as (Fossen, 2011)
MRB + CRB () + MA (r )r + CA (r )r + D(r )r + g = + wind + wave .
(2.7)
The terms in the equations are
relative velocity r = c , where c is the velocity of the current,
interia terms: MRB (rigid body) and MA (added mass)
Coriolis and centripetal terms: CRB (rigid body), and CA (due to added mass),
damping forces: D(r ) = (Dp + Dv ), where Dp is the linear potential damping,
and Dv contains the viscous damping (vortex shedding, skin friction),
restoring forces: g (hydrostatics),
environmental forces:
wind : wind forces,
wave : first order oscillatory waves forces and second order waves forces
(mean drift, slowly varying drift forces, sum-frequency forces) (Faltinsen,
1990).
current: the forces from the current are included in D(r ) due to the relative
velocity r ,
6
2.1 Vessel models
propulsion forces, : the forces generated by the thrusters of the vessel

The Coriolis terms are present in the equation because the kinetics is expressed in BODY
frame, which is a rotating moving reference frame with respect to an inertial reference
frame (Srensen, 2012).
2.1.2
Uncoupled surge dynamics
From Fossen (2011) the uncoupled surge dynamics for a longitudinally symmetric ship,
where m is the mass, Xu is the added mass in surge, Xu and X|u|u | are linear and nonlinear
damping, respectively, can be written
(m Xu )u Xu ur X|u|u |ur |ur = ,
(2.8)
where comprise of the external forces, and the control input. The nonlinear damping
will dominate for higher vessel speeds, and the linear damping dominates for low speeds.
2.1.3
Sway-yaw subsystem
For a constant surge velocity u u0 , a linear sway-yaw subsystem, known as the second
order Nomoto Model, can be written (Fossen, 2011)
M + N (u0 )vr = b
(2.9)
with

m Yv
M=
mxg Yr

mxg Yr
R22 ,
Iz Nr
(2.10)
and N R22 contain the speed dependent terms from C(v) and the linear damping
matrix, D from Eq. (2.7) (Fossen, 2011). The rudder angles are collected in the vector .
Let the in Eq. (2.9) be a scalar. To relate this to the tests with R/V Gunnerus, that has
two rudders, assume that the rudders have an equal angle, and thus modelled as a single
rudder. Looking at the transfer function from the rudder angle to the yaw-rate, it can be
written as (Fossen, 2011)
Ky (1 + T3 s)
r
(s) =
.
(1 + T1 s)(1 + T2 s)
(2.11)
Also, for the sway motion a similar relationship is found as

v
Kv (1 + Tv s)
(s) =
,
(1 + T1 s)(1 + T2 s)
(2.12)
where r is the yaw rate, and v is the sway velocity, Ky , Kv are the steady state gains, and
T1 , T2 , and T3 are the time constants.
7
Chapter 2. Modeling
First order approximation

Fossen (2011) defines a equivalent time constant so that from (2.11) the equivalent time
constant Tr is defined as
Tr := T1 + T2 T3 ,
(2.13)
and from (2.12) the equivalent time constant Tv is defines as

Tv := T1 + T2 Tv ,
(2.14)
and hence, the models given by (2.11) and (2.12) can be approximated by the first order
models as
Ky
r
(s) =
,
1 + Tr s
(2.15)
Kv
v
(s) =
.
1 + Tv s
(2.16)
and
2.1.4
Low speed models
For low speed applications, like dynamic positioning, the linear damping will dominate the
nonlinear part (Fossen, 2011), so D(r ) D, where D is the linear damping matrix.
The Coriolis and centripetal terms can be neglected, also due to low speed.
DP (3 DOF) - Linearized model
For dynamic positioning of a vessel it is common to restrict the workspace to 3 degrees of
freedom (DOF). Only the horizontal plane (surge, sway and yaw) is considered, and this
gives and as

>
= N E

>
= u v r ,
(2.17)
(2.18)
and the rotation matrix reduces to
cos() sin() 0
R() = R() = sin() cos() 0 .
0
0
1
(2.19)
In the linearized DP model the slowly varying drift forces, mean drift forces, and the current are all collected into a bias b. Since current is captured in the bias, the relative velocity
vector is not included in the model anymore (it becomes superfluous) (Fossen, 2011).
2.2 Thruster dynamics
The kinematics and kinetics for the linearized 3 DOF DP model is written as (Fossen,
2011)
= R()
T
M + D = R ()b + + wind + wave1

b = T 1 b + wb ( or b = wb ),
(2.20)
(2.21)
(2.22)
where wb is white Gaussian noise (Srensen, 2012), and wave1 comprise only of the first
order wave forces. The bias force is modeled by a Gauss-Markov process above, and could
also be modeled as a white noise process (Srensen, 2012), which is given in parenthesis
in Eq. (2.22). Here b is a slowly varying disturbance or bias force, and the linear damping
matrix satisfies D > 0.
2.2
Thruster dynamics
This section is based on Srensen (2012), and Smogeli (2006). Please see these references
for more in-depth analysis of thruster control, especially other control strategies such as
torque control, and combined torque and power control. Only speed control will be covered in the following, both because this is the most common strategy, and because it is the
algorithm used for thruster control at Gunnerus (Schultz, 2013).
Thruster control is an important aspect, and good control algorithms can save fuel, and
wear and tear of the equipment. Based on the force demand calculated from the high level
DP controller, or from manual control, the thrust allocation distributes the force demand
to each thruster. Based on the required force given by the thrust allocation, the low level
thruster controller calculates the desired shaft speed of the propeller. Other possible control
strategies include torque control, and power control.
The low level thruster control scheme consists of five main components (Srensen, 2012).
There is a reference generator to ensure physically feasible input, a core controller (feedback controller), a friction feedforward term, a intertia feedforward term, and a torque
saturation block to limit the commanded torque within allowed (and physical) limits. This
commaned torque, Qc is fed into the electric motor, and the motor dynamics can be modelled as a first order model such as
1
Q m =
(Qm Qc ),
Tm
(2.23)
where Qm , Qc , and Tm are the actual motor torque, the commanded motor torque, and the
first order time constant for the motor dynamics, respectively. The shaft dynamics is given
as
Is = Qm Qa K ,
(2.24)
where Is , Qa , K , and are the shaft inertia, the propeller load torque, the friction coefficient, and the shaft speed respectively. Finally, hydrodynamic loss effects will affect the
9
Chapter 2. Modeling
delivered thrust. The thrust force actually delivered after losses have been accounted for
is called the actual thrust Ta . For a DP vessel typical loss effects are (Srensen, 2012)
in-line velocity fluctuations,
cross-coupling drag,
Coanda effect,
ventilation,
in-and-out-of-water-effects,
thruster-thruster interaction.
For R/V Gunnerus ventilation is only a problem in rough weather when the propellers
comes close to the surface, and similarly for in-and-out-of-water-effects (really bad weather).
The tunnel thruster is only effective for low-speed, and will become inefficient at higher
speed due to cross-coupling drag.
The relationship between actual thrust force Ta [N ] and actual torque Qa [N m], with the
shaft speed of the thruster are given by (Srensen, 2012)
Ta = sign(n)KT D4 n2 ,
5 2
Qa = sign(n)KQ D n ,
(2.25)
(2.26)
kg
where n[ 1s ] is the shaft speed, [ m
3 ] is the water density, D[m] the propeller diameter, and
KT [] > 0, and KQ [] > 0 are the thrust and torque coefficients, respectively. Values for
KT and KQ are found from open water tests. Both coefficients are a function of a number
of parameters, the advance ratio, the expanded blade area (ratio), the pitch ratio, and the
number of blades. For Qa the Reynold number and the maximum thickness of the blade
with respect to the chord length is also of importance. The advance ratio Ja is given as
Ja =
Va
,
nD
(2.27)
where Va is the inflow velocity on the propeller.

The power consumption can be found from (Srensen, 2012)
Pa = 2nQa = sign(n)2KQ D5 n3 .
2.2.1
(2.28)
Simplified model
According to Strand et al. (2001), full scale experiments have shown that the thruster
dynamics can be modelled as a first order process with quite satisfactory results. This
gives the thruster dynamics as
Tm = A1
thr (Tm Tc )
10
(2.29)
2.2 Thruster dynamics
where Tm and Tc are obtained and commanded thrust, respectively. Athr R33 is
a diagonal matrix containing the time constants for the thrust in surge, sway, and yaw
direction, respectively. This model is applicable for the step test analysis, and will be used
in the system identification process.
2.2.2
Thrust allocation
The DP controller (or manual control by levers) calculates desired force and moment to
be applied in surge, sway, and yaw. Thrust allocation is the process of distributing this
force demand on the thrusters of the vessel, and the low level thruster controller map
this to a desired shaft speed of the individual thrusters (Srensen, 2012). In the case of
an overactuated vessel, the thrust allocation problem becomes an optimization problem.
Most DP vessels are overactuated. Generally, the thrust force is written as (Fossen, 2011)
= T ()Ku,
(2.30)
where is a vector of the azimuth angles. Let n be the number of DOFs, and r be
the number of thruster, then the thrust configuration matrix T () Rnr geometrically
describes the thruster position and their orientation with respect to the centre of rotation.
The matrix K Rrr is a diagonal matrix containing the thrust coefficients. For speed
controlled thrusters, K and u Rr are given as
K = diag{k1 ,
=
k2 ,
diag{KT 1 D14 ,
kr }
KT 2 D24 ,
|n1 |n1
|n2 |n2
u=
... .
|nr |nr
KT r Dr4 },
(2.31)
(2.32)
For R/V Gunnerus (before retrofit), there are two main fixed pitch propellers (Srensen,
2012) with rudder in the stern, and a tunnel thruster in the bow (no azimuth thrusters).
Given the following definitions
u1 : port main propeller
(2.33)
u2 : starboard main propeller
(2.34)
u3 : tunnel thruster, bow,
(2.35)
and let l1 , l1 and l3 be the moment arms in yaw for u1 , u2 and u3 , respectively. For F1 =
k1 u1 , F2 = k2 u2 , and F3 = k3 u3 , the thrust configuration is shown in Figure 2.1.
11
Chapter 2. Modeling
Figure 2.1: Thruster configuration, Gunnerus - before retrofit (bow to the right)
Be defining the rudder angles as

= 1
T
R2 ,
(2.36)
where 1 and 2 are the rudder angles of the port and starboard rudder, respectively, the
thrust configuration matrix is written as
b11 ()
b12 ()
0
b22 ()
1 .
(2.37)
T = b21 ()
b31 ()l1 b32 ()(l1 ) l3
If the surge speed is small, and the tunnel thruster is active, the yaw motion will not depend
linearly on the rudder angle as in the sway-yaw subsystem of Section 2.1.3.
12
Chapter
Model Identification
3.1
Introduction
Simple system identification tests as step tests, and the logging of other maneuvers such
as different DP maneuvers, are useful to get a rough overview over the performance of the
DP system, and of the thruster configuration. R/V Gunnerus has a planned thruster retrofit,
and results from simple tests can be used to compare the responsiveness, and effectiveness
of the two thruster systems.
Simple identification of performance of the DP system, and time constants found in step
tests could also function as a first benchmark for simulators, and the data found could give
an indication of the simulator performance.
In this chapter the step tests are analyzed using the plot function in MATLAB, and then
the results of the tests are read off manually, using MATLAB plots.
This chapter is structured such that in Section 3.2 system identification by step tests are
discussed, before suggested models used for the identification are proposed in Section 3.3.
Some data that are needed for the tests are unavailable, and this is discussed in Section
3.4. In Section 3.5 the step tests are analyzed, and other test results are presented, such as
the thrust agility tests, DP 4 corner, and pure DP tests (stationkeeping).
3.2
3.2.1
System identification by step tests

First order system
A step test is a system subject to a sudden change in the input reference, a jump/step from
one constant value to another. This section takes its results from Palm (2009) and Smith
13
Chapter 3. Model Identification
and Corripio (2006). For a general 1st order system of the form
ay + by = u,
u
T y + y =
b
(3.1)
(3.2)
where T is the time constant, ( T = ab ). If u is a step input, the response (output) y is given
as
t
y(t) = y(0)e T +
t
U
(1 e T ),
b
(3.3)
where U is the magnitude of the input (the step). This implies two useful results that can
be applied to extract information from a step response of a 1st order system;
as steady state is reached (t ) the steady state gain, Kp can be found from
Kp =
y
,
u
(3.4)
where y and u are the change in output and input, respectively, due to the step
input. For Eq. (3.3) Kp is
Kp =
U
b
y(0)
yss y(0)
=
,
U
U
(3.5)
for yss and y(0) representing the steady state output value and output initial value,
respectively, the time constant can be found from the response curve as the time
difference from the step was applied to the response reaches y(0) + (1 e1 )[yss
y(0)]. This corresponds to the time it takes for the response to reach about 0.632 of
its steady state value.
2nd order system
A general second order system can be written as
Kp
,
(T1 s + 1)(T2 s + 1)
(3.6)
where Kp is the steady state gain, and T1 and T2 are the time constants. It is hard to get
the time constants from a second order response curve. However, there are methods for
finding both time constants for a second order system, such as Smiths method from
Seborg et al. (2004) and Smith (1972). Smiths method utilizes a normal step test, and
looks at the time at 20% and 60% of the response. Another method is that of Rangaiah
and Krishnaswamy (1994) where an overdamped second order process is estimated from a
step test by the use of 3 points, 14%, 55%, 91% (Seborg et al., 2004). However, according
to Smith and Corripio (2006), methods like these have low precision, since the amount of
information contained in a step test is too low. These methods will not be further analyzed.
The first order plus dead time (FOPDT) will instead be of main concern for the step tests.
14
3.3 Proposed identification models for the step tests
First order plus dead time approximation (FOPDT)

The first order plus dead time approximation method has a transfer function G(s) of the
form
G(s) =
Kp s
e ,
Ts + 1
(3.7)
where and T are the dead time and the time constant, respectively. The steady state
gain Kp is found from (3.4). The recommended method of Smith and Corripio (2006) for
approximating a higher order curve as a FOPDT was initially proposed by Smith (1972).
This method is based on the idea that the response curve and the approximated model
should coincide at the two points where the rate of change is high. The two chosen points
are ( + T3 ) and ( + T ). Using the notation of Smith and Corripio (2006) let t1 = + T3 ,
1
and t2 = + T . In other words, t1 is at about (1 e 3 ) 0.283 of the response, and t2 is
1
at 1 e 0.632 of the response. This renders the following solution for the dead time
and the time constant T
3
(t2 t1 )
2
= t2 T
T =
(3.8)
(3.9)
Approximating higher order systems with first order time delay functions
When the system is of higher order, but one time constant is by far the largest, then the
system can be approximated as a first order system with a time delay (Seborg et al., 2004).
According to Smith and Corripio (2006) the following rule of thumb can be used to approximate a higher order model into a first order model: preserve the time constant of the
highest value, and add the other time constants to the time delay. This is the same as
the first order Taylor series approximation given by Seborg et al. (2004) for as one of the
lower time constants, as
1
es .
s + 1
(3.10)
The thruster dynamics is significantly faster than the vessel motion, and this motivates a
model structure of first order plus dead time systems (FOPDT) for the step tests that are
performed. The delay can be used to estimate the time constant of the thruster dynamics.
3.3
Proposed identification models for the step tests
As mentioned in Section 3.2, the information that can be obtained in a step test is limited.
The most advanced result that to some certainty can be extracted is a FOPDT approximation. This could allow for a second order model, where one of the time constants could be
included in the delay. This is the approximation that will be used for all the step tests in
the analysis, and the proposed models are discussed below.
15
The steady state gains are mentioned in the models, but will not be collected in the identification results. The responsiveness of the system is of interest, so the time constants
and delay time are the interesting parameters to compare responsiveness of this thruster
configuration with the retrofit.
In the Laplace-domain (Palm, 2009), let T(s), u(s), n(s), v(s), and r(s) represent thrust,
surge speed, shaft speed of the thrusters, sway speed, and yaw rate [ /s], respectively.
Let T1 , T2 , and so on, represent time constants. All the suggested models below will be
identified identified by the FOPDT approximation of Section 3.2.1
3.3.1
Thrust and shaft speed time constants
For the thrust identification, two time constants are of interest. The time constant for the
shaft speed, and the time constant for the obtained thrust force. The obtained thrust force
is calculated from Eq. (2.25) from the measured shaft speed, since the thrust force is not
measured. These two quantities will reach steady state at the same time, but have different
time constants. This is due to the quadratic mapping between shaft speed and thrust force
from Eq. (2.25).
By the discussion in Section 2.2 the dynamics from commanded thrust force to obtained
shaft speed is of second order, that consists of a first order motor model, and first order
shaft dynamics in cascade. The dynamics from commanded to obtained thrust force is
also of second order since the mapping from shaft speed to thrust force is static. It is also
proposed a first order model by Section 2.2.1, since real life experiments show that this
approximation is satisfactory. In the following, a first order plus dead time approximation
of Section 3.2.1 is applied, and that will capture both time constants of the second order
model. If the delay and the time constant are added together, then this could approximate
the time constant of the first order thruster model. The transfer function from commanded
thrust, Tc (s) to obtained shaft speed n(s) is
Kp
n(s)
=
,
Tc (s)
(1 + T1 s)(1 + T2 s)
(3.11)
and the first order plus dead time approximation is given as

Kp
n(s)
=
es .
Tc (s)
1 + T3 s
(3.12)
The models for thrust force are the same. The first order plus dead time approach will be
used to find the time constant and the dead time of the thruster system.
3.3.2
Speed time constants
In the following let Kp , Kv , and Kr be steady state gains, a time delay, and T a time
constant.
16
3.4 Uavailable data
Proposed model for surge step tests

Using the decoupled surge dynamics of Eq. (2.8), and a first order thruster dynamics of
Eq. (2.29), letting the thruster dynamics be the delay in a FOPDT approximation, the
dynamics from commanded shaft speed n(s) to obtained surge speed u(s) can be written
as
u(s) =
Kp s
e n(s)
1 + Ts
(3.13)
where the model takes n(s) as input and the surge speed u(s) as output. The thruster
dynamics is included as a delay. The delay can further be modelled using the first order
Taylor series approximation of Eq. (3.10) as
es
1
,
(s + 1)
(3.14)
where is an estimate of the time constant of the thruster dynamics.
Proposed model for the sway step responses

The proposed model for the sway step tests is a model from commanded thrust (or commanded shaft speed), T (s) to obtained sway speed v(s). As discussed in Section 2.1.3 a
first order model can be approximated. A delay time is added to model the response of the
first order thruster dynamics (Section 2.2.1),
v(s) =
Kv s
e T (s).
1 + Ts
(3.15)
Proposed model for the yaw step direction test responses

The proposed model for the yaw step tests are a model from commanded moment, W(s)
[kNm] to obtained yaw rate, r [ /s], and as above, by the discussion in Section 2.1.3 a first
order model can be approximated, a delay time is added to model the first order thruster
dynamics, such that the model becomes
r(s) =
3.4
Kr s
e W (s).
1 + Ts
(3.16)
Uavailable data
Some data are unavailable from the sea trials, and it makes the data gathered in Section
3.5 incomplete. This section summarize what data that is lacking for thrust and power, and
how this is handled.
17
3.4.1
Thrust data
Thruster mapping for the main thrusters

For the sea trials done in August, the calculated thrust force is gathered from the DP
system. However, for the November sea trials, no thrust force estimates are available.
Kongsberg Maritime have supplied their mapping (Schultz, 2013) for the thrusters at Gunnerus, and the mapping they use for the two main thruster is similar to Eq. (2.25), where
c = KT D4 , and it is written as
Ta = cn2 ,
(3.17)
where Ta is in [kN ], and n in rotations per second [ 1s ], and c is a constant coefficient, that
vary with the sign of the shaft speed. The value of c for bollard pull (Smogeli, 2006) and
free run condition is shown in Table 3.1 (Schultz, 2013)
Table 3.1: KM thrustmapping
Bollard condition
c
n>0
6.6
n<0
-4.1
n>0
9.4
n<0
-6.4
Free run condition

c
Also, in order to have a good estimate of the thrust mapping in stationkeeping and DP
maneuvers, the results of the pure DP logs from August has been back-calculated to find
the mapping used. The measurement series of thrust force and shaft speed are used to find
the constant coeffcient c for both positive and negative shaft speed, and the values found
are included in Table 3.2. Note that the values are in between the bollard pull, and free run
condition values.
Table 3.2: Back-calculated thrust mapping from
n>0
8.13
n<0
-5.8
For the zero speed DP maneuvres completed in November, the mapping of Table 3.2 could
have been applied to find an estimate of the thrust force for the two main propellers.
Unfortunately, MARINTEK only measure the absolute value of the shaft speed for November data, so the rotational direction of the shaft is unknown. From the pure DP manuevers
performed in calm sea (August) in Section 3.5.3, often one propeller has positive rotation
while the other has negative (see Table 3.11 and 3.12 where one thruster generate less
thrust than the other. This is simply due to rotational direction, since they have about the
18
3.5 Identification results
same rotational speed (magnitude)). Therefore, in order to have an estimate of the thrust
force, a mean of the two c values will be used. So in the mapping from n2 [rps2 ] to thrust
force [kN ] from Eq. (3.17) a value of
c=
8.13 + 5.8
= 6.97,
2
(3.18)
will be used. This gives a rough estimate of the thrust force. However, since the mapping
does not take thrust losses into account, it is already quite rough. This mapping is used
to find an estimate of the thrust force for the November sea trials for the DP 4 corner in
Section 3.5.2, and for pure DP in Section 3.5.3.
Tunnel thruster
The thrust force for the tunnel thruster is available for the sea trials at calm sea (August),
but not for the trials at sea in November. The shaft speed of the tunnel thruster is however
logged in November, but Kongsberg did not have the mapping for this thruster since it
is delivered by another vendor (Schultz, 2013), and the mapping for this thruster has not
been obtained, unfortunately. The thrust mapping could have been back-calculated from
August thrust data, but unfortunately, the shaft speed was not logged at that time.
3.4.2
Power data
For power consumption MARINTEK measured the shaft torque for the two main propellers, and found the power from Eq. (2.28). This is for both the sea trials. For the tunnel
thruster the power consumption is not measured, and is therefore not available.
3.5
Identification results
For both the pure DP tests of Section 3.5.3, and the DP 4 corner of Section 3.5.2, the
environmental conditions are logged for each maneuver. That is, wind speed and direction
is found from a measurement series, and the size, and direction of the waves are manually
read off just before the DP maneuver, and the values were collected from a weather buoy
that MARINTEK logged. The DP current is the value estimated by the DP system, and
is read off the DP operator screen.
3.5.1
Thrust agility tests
As mentioned in Section 1.1.1, for the thrust agility tests a step in thrust force is applied
for heading of 0 , 45 , 90 , 135 , 180 , 135 , 90 , and 45 . This gives a plot
of responsiveness of the thruster system. Due to lack of setpoint measurements, power
feedback from MARINTEK is used to determine when the steps were applied.
19
For the agility plots the total response is considered. That is, the time constant found
include the dead time, such that the total system response is plotted in the agility plots.
That is, the time constant found is the time of 0.632 of the response (see Section 3.2).
Time constants from the tunnel thruster are highly inaccurate, because of low sampling (1
Hz). Both the start of the step, and the time where the response is at 0.632 of the response
is quite uncertain. Still, the time constants found give an indication on the responsiveness
of the tunnel thruster, so they are included.
Agility tests with no yawing moment applied
For no yawing moment, two thrust agility tests were performed. One where half of maximum available thrust was applied, and one where maximum thrust was applied. The
environmental conditions for both these two tests were about equal. The wind speed and
direction has been found from the DP data, and the current and waves were found before
the tests were performed. The current is the DP current read from the DP system on board
the vessel (the operator station).
The environmental conditions are found in Table 3.3, and the thrust time constants are
found in Figure 3.1, and 3.2 for half and full thrust, respectively. On the radial axis time
in seconds is plotted, and the time constant is shown for the two main thrusters (port and
startboard), and the tunnel thruster. For direction 0 , and 180 the tunnel thruster is not
applicable.
Table 3.3: Environmental conditions thrust agility, no yawing moment
Environment
Wind, speed
Average
Std.
3.3
0.5
Wind, direction [ ]
Average
Std.
22.9
9.0
Waves
Hs [m]
Direction [ ]
0.1
South/SW
Current
Speed [m/s] Direction [ ]
0.2
318
20
Figure 3.1: Agility plot, thrust time constant, no yawing moment - half thrust. Wind of 3.3 m
s
(purple)
(green) and current of 0.2 m
s
(green)
Figure 3.2: Agility plot, thrust time constant, no yawing moment - full thrust. Wind of 3.3 m
s
and current of 0.2 m
(purple)
s
21
Agility test with applied yawing moment

A thrust agility tests was performed with an applied yawing moment of 100%. The thrust
applied was 100% of available thrust force. The environmental conditions has been found
similarly to Section 3.5.1, and is found in Table 3.4, and the thrust time constants are found
in Figure 3.3.
Table 3.4: Environmental conditions thrust agility, with yawing moment
Environment
Wind, speed
Average
Std.
4.3
0.5
Wind, direction [ ]
Average
Std.
26.8
6.9
Waves
Hs [m]
Direction [ ]
0.1
S/SW
Current
0.2
318
(green) and
Figure 3.3: Agility plot - thrust time constant, with yawing moment. Wind of 4.3 m
s
current of 0.2 m
(purple)
s
22
3.5.2
DP 4 corner tests
In order to compare the existing thruster configuration with the retrofit, a maneuvre called
DP 4 corner was performed. This is to compare responsiveness of the thrusters, position
accuracy, thrust usage, and power consumption. This maneuvre is illustrated in Figure 3.4,
and it consists of four parts:
1. Ahead, heading same: From DP steady state with heading north, change desired
position 50 m ahead, keeping same heading,
2. Crab port, heading same: At the instant the vessel reach the operating circle of
the new setpoint, position setpoint is changed to 50 m port (west) of the vessel,
3. Astern, heading change: At the instant the vessel reach the operating circle of
the new setpoint, the new desired setpoint is 50m eastern (south) and heading set to
90 (west),
4. Astern, heading same: AS the instant the vessel reach the operating circle of the
new setpoint, new setpoint is set to be 50 m astern (east, back to original position),
with heading still 90 (west).
Figure 3.4: DP 4 corner, illustration
DP 4 corner - calm sea

During the sea trials in calm sea (August) one DP 4 corner was performed. The log from
the DP system is available, but MARINTEKs log contains only half of the maneuver. The
DP data is used to find the duration. The setpoints in position error are available, so the
duration estimates are quite accurate. However, a lot of measurements are unavailable for
the first half of the maneuver. Even though the DP system logged shaft speed, and thrust
force for the sea trials in August, these values were collected by MARINTEK as the trial
23
was ongoing. Because of this, thrust force for the first half of this maneuver is unavailable,
but power consumption is available for the entire maneuver, since this was collected from
the DP system log. See Table 3.5 for the environmental conditions, and Table 3.6 for the
thrust and power consumption.
Table 3.5: DP 4 corner - environmental conditions,calm sea - 3035
Environment
Wind, speed [m/s]
Average
Std.
5.3
0.3
Wind, direction [ ]
Average
Std.
42.2
2.2
Waves
Hs [m]
Direction [ ]
0.2
N
Current
0.3
147
Table 3.6: DP 4 corner maneuver - calm sea

Maneuvre
Duration [s]
134
214
142
134
Accumulated thrust force [kNs]

SB
Port Tunnel
Total
770 1230
1080
3080
600 1100
960
2660
Accumulated power [kJ]

SB
Port
Tunnel Total
690
600
1290
4750
4380
9130
3230
4520
7750
3620
3790
7410
Ahead, heading same

Crab port, heading same
Astern, heading change
Astern, heading same
Total
624
12290
13290
25580
DP 4 corner - at sea
Two DP 4 corner tests were performed in the November sea trials, and only MARINTEK
logged (no logs from the DP system itself). MARINTEK performed extensive logs containing quite good North and East position measurements, shaft speed measurements of
all thrusters, and a calculated power consumption of the two main thrusters. Setpoints are
therefore unavailable, but in order to find the duration of each maneuver, the north-east
plot is used in combination with the plots of north position, east position and heading. The
duration is not very accurate, but should be correct to within some seconds of margin.
The N-E plot of the first test is shown in Figure 3.5. Note that in the last part of the
maneuver there is an overshoot in east-position before the vessel returns to desired east
position. For the logs in Table 3.7, and 3.8, the duration is stopped when the correct east
24
value is reached the first time (neglecting the overshoot), to get a better estimate of the
power consumption as if there were no overshoot. It should also be noted that the vessel
does not perfectly track the square as it was supposed to do.
2020
2010
North position [m]
2000
1990
1980
1970
1960
1950
2950
2960
2970
2980
2990
3000
East position [m]
3010
3020
3030
Figure 3.5: NE-plot of DP 4 corner test 1, at sea
Table 3.7: DP 4 corner - environmental conditions, at sea - test 1
Environment
Wind, speed [m/s]
Average
Std.
2.9
0.5
Wind, direction [ ]
Average
Std.
178.0
7.5
Waves
Hs [m]
Direction [ ]
1.5
345
Current
0.2
N/NW
25
Table 3.8: DP 4 corner - maneuver, at sea - test 1

Maneuvre
Duration [s]
117
180
163
85

SB
Port Tunnel
Total
310
290
600
1430 1910
3340
530
570
1100
1850 2270
4120

SB
Port
Tunnel Total
630
410
1040
4960
6720
11680
1510
1190
2700
8770 12630
21400
Ahead, heading same

Total
545
4120
15870
5040
9160
20950
2020
2010
North position [m]
2000
1990
1980
1970
1960
1950
1940
2940
2950
2960
2970
2980
2990
East position [m]
3000
3010
3020
Figure 3.6: NE-plot of DP 4 corner test 2, at sea
Table 3.9: DP 4 corner - environmental conditions, at sea - test 2
Environment
Wind, speed [m/s
Average
Std.
2.8
0.4
Wind, direction [ ]
Average
Std.
174.3
6.8
Waves
Hs [m]
Direction [ ]
1.5
345
Current
0.4
N/NW
26
36820
Table 3.10: DP 4 corner - maneuver, at sea - test 2

Maneuvre
Duration [s]
Ahead, heading same

Total
3.5.3
101
207
120
118

SB
Port Tunnel
Total
310
370
680
1050 1300
2350
640
460
1100
2890 3240
6130

SB
Port
Tunnel Total
660
850
1510
3500
3670
7170
1800
900
2700
13820 18600
32420
546
4890
19780
5370
10260
24020
43800
Pure DP tests
For the pure DP tests, the vessel is simply trying to maintain position for a five minute
interval. There were three tests performed in calm water (in August), and nine tests at
sea, with different heading. In August these tests were performed right before other maneuvers, such as the DP 4 corner maneuver of Section 3.5.2. For the pure DP maneuvers
performed in at sea, the table names indicate a heading of the vessel. This heading name
is the observed (by the operator) heading relative to the mean environmental forces, and
therefore it might deviate from the direction of the environmental forces reported in the
table.
The tests are performed at slightly varying environmental conditions. The intention of the
logging of the thrust, and power usage is to perform the same maneuver after the retrofit,
and compare.
The thrust force in the tables below is Ta , the delivered thrust (see Section 2.2). This
is why two average thrust values can be different when power consumption is about the
same; the propellers rotate in opposite directions, and forward is optimal.
Results - calm sea
For the test in August, the North and East measurements are sampled too low, and give a
poor estimate of the standard deviation in the position error. The heading measurements
on the other hand, are useful, and give a indication of DP performance.
27
Table 3.11: Pure DP, calm sea - test 1
Environment
Wind, speed [m/s]
Average
Std.
2.6
0.4
Wind, direction [ ]
Average
Std.
18 .0
4.6
Waves
Hs [m]
Direction [ ]
negligible
Current
Speed [m/s]
Direction [ ]
0.3
325 (during)-334(after)
DP system performance
Position acc. (Std.) [m]
North
East
Heading
Average [ ]
Std [ ]
62.1
0.4
Average thruster force [kN]
SB thruster Port thruster Tunnel thruster
5.5
3.9
2.7
Total average: 12.1 kN
2860
3500
Total acc. power:: 6360 kJ (ex. tunnel thr.)
Environment
Wind, speed [m/s]
Average
Std.
2.2
0.4
Wind, direction [ ]
Average
Std.
3.3
9.4
Waves
Hs [m]
Direction [ ]
small ripples
Current
Speed [m/s]
Direction [ ]
0.2
223 (before)-241 (after)
28
North
East
Heading
Average [ ]
Std [ ]
130.2
0.4
4.0
5.5
2.9
3100
3280
Environment
Wind, speed [m/s]
Average
Std.
3.9
0.3
Wind, direction [ ]
Average
Std.
43.8
2.6
Waves
Hs [m]
Direction [ ]
0.2
North
Current
0.1
235
North
East
Heading
Average [ ]
Std [ ]
358.8
0.79
0.67
0.65
0.53
129
545
Results - at sea
For the November sea trials the position measurements are better, and they are therefore
included in the table. As mentioned in Section 3.4.1, the thrust force is not available for
these sea trials, and the thrust mapping of Section 3.4.1 is used to find the thrust force
estimate.
Table 3.14: Pure DP, at sea - test 1, relative heading 0 degrees
Environment
Wind, speed [m/s]
Average
Std.
5.2
1.1
Wind, direction [ ]
Average
Std.
79.1
7.7
Waves
Hs [m]
Direction [ ]
2.1
N-NW
Current (meas. by Kongsberg)
Speed [m/s]
Direction [ ]
0.4
North
East
0.8
0.6
Heading
Average [ ]
Std [ ]
336.1
0.8
1.2
1.0

860
413
29
Environment
Wind, speed [m/s]
Average
Std.
3.3
0.5
Wind, direction [ ]
Average
Std.
72.5
13.9
Waves
Hs [m]
Direction [ ]
2.3
340
Speed [m/s]
Direction [ ]
0.3
98
North
East
0.6
0.9
Heading
Average [ ]
Std [ ]
65.5
0.8
4.3
3.9
Total average: 8.2 kN (ex. tunnel thr.)

3660
3590
Environment
Wind, speed [m/s]
Average
Std.
4.2
0.6
Wind, direction [ ]
Average
Std.
190.8
7.5
Waves
Hs [m]
Direction [ ]
1.7
14
Speed [m/s]
Direction [ ]
0.36
336
30
North
East
0.7
0.3
Heading
Average [ ]
Std [ ]
11.0
0.78
5.0
5.6

3660
3590
Environment
Wind, speed [m/s]
Average
Std.
4.3
1.0
Wind, direction [ ]
Average
Std.
204.8
5.6
Waves
Hs [m]
Direction [ ]
1.7
14
Current
0.2
247
North
East
0.5
0.5
Heading
Average [ ]
Std [ ]
101.7
0.9
1.4
2.0

1060
1240
Environment
Wind, speed [m/s]
Average
Std.
4.2
0.5
Wind, direction [ ]
Average
Std.
215.3
14.0
Waves
Hs [m]
Direction [ ]
1.7
14
Current
0.25
300
North
East
0.5
0.4
Heading
Average [ ]
Std [ ]
55.3
0.7
2.8
2.4

Accumulated power [kJ] over 5 min
2020
1090
31
Environment
Wind, speed [m/s]
Average
Std.
4.7
0.7
Wind, direction [ ]
Average
Std.
210.3
7.1
Waves
Hs [m]
Direction [ ]
1.7
14
Current
0.3
191
North
East
0.5
0.4
Heading
Average [ ]
Std [ ]
146.0
0.8
5.5
5.6

5480
5560
Environment
Wind, speed [m/s]
Average
Std.
4.1
0.5
Wind, direction [ ]
Average
Std.
212.3
9.0
Waves
Hs [m]
Direction [ ]
1.7
14
Current
0.21
231
32
North
East
0.6
0.3
Heading
Average [ ]
Std [ ]
190.1
0.7
2.0
1.6

1220
810
Environment
Wind, speed [m/s]
Average
Std.
6.9
0.6
Wind, direction [ ]
Average
Std.
166.0
8.1
Waves
Hs [m]
Direction [ ]
1.9
West
Current
0.5
240
North
East
0.7
0.8
Heading
Average [ ]
Std [ ]
241.7
0.7
1.1
1.5

790
590
Table 3.22: Pure DP, at sea - test 9, relative heading -90 degrees
Environment
Wind, speed [m/s]
Average
Std.
6.9
0.7
Wind, direction [ ]
Average
Std.
167.6
7.3
Waves
Hs [m]
Direction [ ]
1.85
West
Current
0.6
200
3.5.4
North
East
0.8
0.9
Heading
Average [ ]
Std [ ]
150.4
2.8
23.3
21.2

29620
28140
Step tests
There were several step tests performed, all at the calm water sea trials in August. One
intention is to find the responsiveness of the thrusters. That is, the time constant for the
shaft speed, and for the thrust force. Also, the goal is to find the time constants for the vessel response in surge, sway, and yaw. The shaft speed and thrust time constants are found
for surge and sway. The tunnel thruster response is not included in this analysis since
the values for the tunnel thruster are of low sampling rate, and it is difficult to determine
33
when steps are applied. In the agility plots of Section (3.5.1) the tunnel thruster response
(including dead time) is plotted, and it can be verified that it is faster than the main propellers. Therefore, the response of the main propellers has the highest time constant, and
the response of the tunnel thruster is not of great interest.
For the surge and sway step tests the proposed model of Eq. (3.12) is applied to find the
shaft speed time constant and dead time, and also the thrust force time constant and dead
time.
For surge, the proposed model of Eq. (3.13) from commanded thrust (or shaft speed) to
surge speed is applied for the speed time constants and delay time. The same approach is
used for the sway speed time constant, by using the model of Eq. (3.15). For the yaw step
tests, time constants and delays are found for the yaw rate according to Eq. (3.16). Also,
the rudder time constants are found according to the FOPDT approach of Section 3.2.1.
For all the tables below, the time constant is represented by T , and Tn is the time constant
for the shaft speed, and Tthrust is the time constant for the thrust. Similarly is a time
delay (or dead time), and n is the delay for shaft speed, and thrust is the delay for thrust.
For the rudder time constants, Tsb and Tport represent the time constants for starboard and
port rudder, respectively. Similarly for the delays sb and port .
About how tests are read

For all the step tests no setpoints are available. To figure out when steps are applied, the
curves for power feedback is applied, since power responds quickly. In the cases where
the power curve is unclear the shaft speed feedback curve is used.
Surge step tests

Several step tests in surge were performed for both forward and backward speed. For
forward speed the step tests were run for a step in thrust magnitude [%] of 05%, 510%,
10 15%, 15 25%, 25 40%, 40 0%, 0 80% and 80 0%. Figure 3.7 and 3.8 show
the surge velocity, shaft speed feedback, and power feedback for the forward speed step
tests. Table 3.23 present the time constant and delay time for the shaft speed and thrust
feedback, whereas Table 3.24 show the speed time constant and delay. Some of the step
tests results are omitted, due to uncertainty, or for the case of 40 0%, and 80 0% they
are not included since they did not reach steady state before the tests were finished.
For backward speed the step tests were run for a step in thrust magnitude [%] of 0 5%,
5 10%, 10 15%, 15 20%, 20 0%, 0 40%, and 40 0%. Figure 3.9 and 3.10 show
the surge velocity, shaft speed feedback, and power feedback for the backward speed step
tests. Table 3.25 present the time constant and delay time for the shaft speed and thrust
feedback, whereas Table 3.26 show the speed time constant and delay. As with the forward
speed test results, the tests with poor data is not included in the tables.
34
Surge velocity
m/s
2
1
0
200
400
600
800
time [s]
1000
1200
1400
1600
600
800
time [s]
1000
1200
1400
1600
600
800
time [s]
1000
1200
1400
1600
80
Shaft speed feedback (sb)
rpm
60
40
20
0
200
400
50
Power feedback (sb)
kW
40
30
20
10
0
0
200
400
Figure 3.7: Step test in surge - forward velocity. Surge velocity (top, blue), shaft speed feedback for
starboard thruster (red), and power feedback for starboard thruster (black). Power and shaft speed
feedback is used to indicate when the input is applied. The steps shown in the figure are test number
0-5,5-10,10-15,15-25, 25-40, and 40-0 (ref. Table (ref. Table 3.23, and 3.24).
Surge velocity
m/s
4
2
0
2
50
100
150
200
time [s]
250
300
350
400
200
rpm
150
100
50
0
50
100
150
200
time [s]
250
300
350
400
300
Power feedback (sb)
kW
200
100
0
100
50
100
150
200
time [s]
250
300
350
400
Figure 3.8: Step test in surge - forward velocity. Surge velocity (top, blue), shaft speed feedback for
starboard thruster (red), and power feedback for starboard thruster (black). Power and shaft speed
feedback is used to indicate when the input is applied. The step shown in the figure is for 0-80 (ref.
Table 3.23, and 3.24).
35
Table 3.23: Surge step tests, forward speed - thrust and shaft speed time constants
Step mag. [%]

0-80
25-40
15-25
10-15
Tn [s]
5.5
1.1
1.0
1.4
n [s]
2.9
2.0
1.9
11.6
Tthrust [s]
4.7
2.2
1.4
0.3
thrust [s]
5.6
1.5
1.8
12.7
Table 3.24: Surge step tests, forward speed - speed time constants
Step mag. [%]

0-80
25-40
15-25
10-15
5-10
0-5
T [s]
16.4
34.2
59.2
135.0
63.3
65.5
0.2
[s]
8.5
1.4
1.8
7.2
5.0
6.0
Surge velocity
m/s
0
0.2
0.4
0.6
200
400
600
time [s]
800
1000
1200
600
time [s]
800
1000
1200
60

rpm
40
20
0
200
400
15
Power feedback (sb)
kW
10
5
0
5
200
400
600
time [s]
800
1000
1200
Figure 3.9: Step test in surge - backward velocity. Surge velocity (top, blue), shaft speed feedback
for starboard thruster (red), and power feedback for starboard thruster (black). Power and shaft speed
feedback is used to indicate when the input is applied. The steps shown in the figure are test number
0-5,5-10,10-15,15-20, and 20-0 (ref. Table Table 3.25, and 3.26).
36
m/s
0
1
Surge velocity
2
50
100
150
200
250
300
350
400
450
time [s]
80
rpm
60
40
20
0
50
100
150
200
250
300
350
400
450
time [s]
60
Power feedback (sb)
kW
40
20
0
20
50
100
150
200
250
300
350
400
450
time [s]
Figure 3.10: Step test in surge - backward velocity. Surge velocity (top, blue), shaft speed feedback
for starboard thruster (red), and power feedback for starboard thruster (black). Power and shaft speed
feedback is used to indicate when the input is applied. The step shown in the figure is for 0-40 (ref.
Table 3.25, and 3.26).
Table 3.25: Surge step tests, backward speed - thrust and shaft speed time constants
Step mag. [%]

0-40
15-20
10-15
5-10
20-0
Tn [s]
1.6
4.4
2.6
0.5
2.7
n [s]
2.7
1.2
4.3
5.4
2.3
Tthrust [s]
1.6
4.4
0.6
0.8
1.2
thrust [s]
3.4
1.3
6.4
5.3
2.4
Table 3.26: Surge step tests, backward speed - speed time constants
Step mag. [%]

0-40
15-20
10-15
5-10
20-0
T [s]
30.7
36.6
76.6
64.6
51.4
[s]
3.1
5.5
0.6
1.4
9.0
Sway step tests

For the sway step tests, step tests were performed for starboard and port side. For the
starboard tests steps of thrust magnitude of [%] of 0 25%, 0 50%, 0 75%, 0
37
100% and 100 0% were run. For port side, all the same tests except 0 75% were
performed.
As seen from Figure 3.11 and 3.13 the surge velocity is considerable for the step tests.
This is probably because of the thruster configuration. The vessel has one tunnel thruster,
and two propellers in the stern, and the vessel needs a certain surge velocity to have full
flexibility to move in sway. If only the tunnel thruster is used it would induce a moment
that have to be counteracted by the stern propellers.
Due to the high surge velocity the total speed is used. That is,
v 2 + u2 ,
(3.19)
where u and v are the surge and sway velocity, respectively. The plots for total speed,
shown in Figure 3.11 and 3.13 indicate a clearer step response shape, but for some of the
step tests, it does not seem like steady state is reached. The values found are therefore of
lower quality. The values reported for starboard are found in Table 3.27 for the thrust, and
shaft speed time constants, and in Table 3.28 for the speed time constants. For the port
step response the values are found in Table 3.29 for the thrust time constants, and in Table
3.30 for the speed time constants.
Sway velocity
m/s
0.5
0
0.5
m/s
0.5
200
400
600
800
time [s]
1000
1200
1400
1600
600
800
time [s]
1000
1200
1400
1600
600
800
time [s]
1000
1200
1400
1600
Surge velocity
0
0.5
1
0
200
400
1.5
Total speed
m/s
1
0.5
0
200
400
Figure 3.11: Step test in sway - starboard. Sway velocity (top, blue), surge velocity (red), and total
speed (black). The steps shown in the figure are test number 0-25,0-50,0-75,0-100, and 100-0 of
(ref. Table 3.27, and 3.28).
38
1.5
Total speed
m/s
1
0.5
0
200
400
600
800
time [s]
1000
1200
1400
1600
600
800
time [s]
1000
1200
1400
1600
600
800
time [s]
1000
1200
1400
1600
200
rpm
150
100
50
0
200
400
400
Power feedback (sb)

kW
200
0
200
200
400
Figure 3.12: Step test in sway - starboard. Total speed (top, blue), shaft speed feedback for starboard
thruster (red), and power feedback for starboard thruster (black). Power and shaft speed feedback
is used to indicate when the input is applied. The steps shown in the figure are test number 0-25,050,0-75,0-100 (ref. Table 3.27, and 3.28).
Table 3.27: Sway step tests, starboard speed - thrust and shaft speed time constants for sb thruster
Step mag. [%]

0-100
0-75
0-50
0-25
100-0
Tn [s]
4.7
3.9
2.5
2.1
4.2
n [s]
3.3
2.7
3.0
1.2
1.6
Tthrust [s]
4.4
3.5
2.5
1.2
2.8
thrust [s]
4.9
4.1
3.8
2.7
1.1
Table 3.28: Sway step tests, main propellers, starboard - speed time constants
Step mag. [%]

0-100
0-75
0-50
0-25
T [s]
10.7
12.5
15.6
14.4
[s]
2.7
2.6
1.1
4.0
39
0.5
Sway velocity
m/s
0
0.5
1
100
200
300
400
500
time [s]
600
700
800
900
1000
0.6
Surge velocity
m/s
0.4
0.2
0
0.2
100
200
300
400
500
time [s]
600
700
800
900
1000
300
400
500
time [s]
600
700
800
900
1000
0.8
Total speed
m/s
0.6
0.4
0.2
0
100
200
Figure 3.13: Step test in sway - port. Sway velocity (top, blue), surge velocity (red), and total speed
(black). The steps shown in the figure are test number 0-25,0-50,0-100 (ref. Table 3.29, and 3.30).
0.8
Total speed
m/s
0.6
0.4
0.2
0
rpm
150
100
200
300
400
500
time [s]
600
700
800
900
1000
400
500
time [s]
600
700
800
900
1000
400
500
time [s]
600
700
800
900
1000
100
50
0
100
200
300
300
Power feedback (sb)
kW
200
100
0
100
100
200
300
Figure 3.14: Step test in sway - port. Total speed (top, blue), shaft speed feedback for starboard
thruster (red), and power feedback for starboard thruster (black). Power and shaft speed feedback
is used to indicate when the input is applied. The steps shown in the figure are test number 0-25,050,0-100 (ref. Table 3.29, and 3.30).
40
Table 3.29: Sway port step tests - thrust and shaft speed time constants for port thruster
Step mag. [%]

0-100
0-50
0-25
100-0
Tn [s]
4.5
2.4
1.0
4.2
n [s]
2.9
2.9
3.4
1.5
Tthrust [s]
4.0
2.7
1.1
2.7
thrust [s]
5.0
3.2
3.4
1.0
Table 3.30: Sway step tests, main propellers, port - time constant
Step mag. [%]

0-100
0-50
0-25
T [s]
9.0
12.8
18.8
[s]
2.5
1.3
1.0
Yaw step tests
Four yaw step tests were performed. Two where the vessel turned starboard (at speeds of
3 knots, and 6 knots), and two where the vessel turned port (at speeds of 3 knots, and 6
knots). For all the yaw step tests both the rudder angles are first stepped to 20 degrees,
and then the rudders are stepped to their maximum angle (about 47 degrees). Notice for
instance from Figure 3.15 that the step in rudder angle from 20 to 47 degrees gives an
initial pulse in yaw rate, but does not significantly change the steady state yaw rate value.
The yaw rate time constant found in Table 3.31 is for the step response from 0 to 20 degrees
rudder angle, and for Table 3.32 for 47 to 0 degrees rudder angle. For the step from 20
to 47 degrees a time constant is not found since the response is so unclear and does not
resemble a normal step response. The rudder angle change is used to indicate when a step
is applied. All the steady state yaw rates oscillate after the new value is reached, and the
steady state value is chosen as the observed mean value of this oscillation.
41
2
1.5
deg/s
1
0.5
0
Yaw rate
0.5
0
200
400
600
time [s]
800
1000
1200
400
600
time [s]
800
1000
1200
400
600
time [s]
800
1000
1200
20
deg
0
20
40
60
Rudder angle starboard

0
200
deg
20
40
Rudder angle port

60
200
Figure 3.15: Step test in yaw - starboard at 3kn. Yaw rate (top, blue), rudder angle starboard (red),
and rudder angle port (black). The steps shown in the figure are for rudder angle steps of 0-50 (ref.
Table ), 50-100, and 100-0 (ref. Table 3.31, 3.32, and 3.33).
deg/s
0
1
2
Yaw rate
0
50
100
150
200
250
300
350
400
450
250
300
350
400
450
250
300
350
400
450
time [s]
60
deg
40
20
0
20
0
50
100
150
200
time [s]
60
deg
40
20
0
Rudder angle port
20
0
50
100
150
200
time [s]
Figure 3.16: Step test in yaw - port at 3kn. Yaw rate (top, blue), rudder angle starboard (red), and
rudder angle port (black). The steps shown in the figure are for rudder angle steps of 0-50 (ref. Table
), 50-100, and 100-0 (ref. Table 3.31, 3.32, and 3.34).
42
deg/s
Yaw rate
2
1
0
100
200
300
time [s]
400
100
200
300
time [s]
400
100
200
300
time [s]
400
500
600
20
deg
0
20
40
60

500
600
20
deg
0
20
40
60
Rudder angle port

500
600
Figure 3.17: Step test in yaw - startboard at 6kn. Yaw rate (top, blue), rudder angle starboard (red),
and rudder angle port (black). The steps shown in the figure are for rudder angle steps of 0-50,
50-100, and 100-0 (ref. Table 3.31, 3.32, and 3.35).
deg/s
0
1
2
3
Yaw rate
0
50
100
150
200
250
300
350
400
450
250
300
350
400
450
250
300
350
400
450
time [s]
60
deg
40
20
0
20
0
50
100
150
200
time [s]
60
deg
40
20
0
Rudder angle port
20
0
50
100
150
200
time [s]
Figure 3.18: Step test in yaw - port at 6kn. Yaw rate (top, blue), rudder angle starboard (red), and
rudder angle port (black). The steps shown in the figure are for rudder angle steps of 0-50, 50-100,
and 100-0 (ref. Table 3.31, 3.32, and 3.36).
Yaw rate time constant The yaw rate time constants for 0 50[%] and 100 0[%] are
shown in Table 3.31 and 3.32, respectively.
43
Table 3.31: Yaw step tests - time constant 0-50
Step 0-50 [%]

Starboard, 3kn
Port, 3kn
Starboard, 6kn
Port, 6kn
T [s]
4.7
5.6
2.9
3.0
[s]
2.1
2.0
2.0
1.9
Table 3.32: Yaw step tests - time constant 100-0
Step 100-0 [%]

Starboard, 3kn
Port, 3kn
Starboard, 6kn
Port, 6kn
T [s]
10.5
7.4
4.7
5.1
[s]
5.2
5.6
6.1
5.6
Rudder time constant Based on the plot of rudder angle, and assuming that the steps
were applied when the rudder starts to change angle, the time constant is found for the
rudder response for the different yaw maneuvers in the tables below.
Table 3.33: Yaw step tests - rudder time constant, starboard 3 kn
Step mag. [%]

0-50
50-100
100-0
Tsb [s]
1.8
2.4
4.2
sb
0.8
1.0
1.4
Tport [s]
1.7
2.3
4.2
port [s]
1.3
1.0
1.5
Table 3.34: Yaw step tests - rudder time constant, port 3 kn
Step mag. [%]

0-50
50-100
100-0
Tsb [s]
1.8
2.4
4.2
sb
0.9
1.0
1.4
Tport [s]
2.0
2.4
4.4
port [s]
0.7
1.2
1.5
Table 3.35: Yaw step tests - rudder time constant, starboard 6 kn
Step mag. [%]

0-50
50-100
100-0
44
Tsb [s]
1.8
2.1
4.2
sb
1.0
1.2
1.5
Tport [s]
1.8
2.3
4.1
port [s]
1.0
1.1
1.5
3.6 Conclusions and further work
Table 3.36: Yaw step tests - rudder time constant, port 6 kn
Step mag. [%]

0-50
50-100
100-0
3.6
Tsb [s]
1.8
2.3
4.1
sb
0.9
1.1
1.4
Tport [s]
1.7
2.3
4.1
port [s]
1.0
1.0
1.4
Conclusions and further work
The data logging of this chapter is simple, and the results give a rough indication of performance of the system. The different maneuvers give an overall indication of DP performance, the thruster performance, and how the thruster system in combination with the
vessel hull behave.
The agility tests give a simple and graphical measure of thruster response time. A DP 4
corner for the new thruster configuration could give good indications on power consumption, and total (accumulated) thruster force, and the same for the pure DP maneuvers. It
is a drawback that much data for the tunnel thruster is lacking. Because of this, the DP 4
corner would only give a rough performance comparison. Also, the thrust mapping used
is somewhat simple, and does not take loss effects into account. This means that the thrust
force calculated is more accurate under normal (no significant thrust losses), and zero
speed conditions. For the November sea trials, the rotational direction of the shaft speed is
not known, so the thrust data is not very reliable for these maneuvers, but it gives a rough
indication.
For further work the same tests should be run with the new thruster configuration, and
the new tests should be performed under similar environmental conditions. For the agility
plots, the DP 4 corner, and the pure DP maneuvers the same data as the data reported in
this chapter should be reported for the new system to compare. Similarly for the step tests
it could be interesting to see how the new thurster system responds, and how the vessel
responds. For the step tests, especially for the sway step tests, longer tests should be run,
such that it is certain that steady state is reached.
45
Chapter
Lever arm estimation for GNSS and

MRU sensors
This Chapter considers the topic of estimating the body frame position of the GNSS and
MRU sensors. For the GNSS position an observer design and an adaptive scheme are
proposed and analyzed. The estimation designs are tested using numerical simulations
and experimental data from the Gunnerus sea trials. A similar observer is proposed for the
MRU positions, and experimental data from the sea trials are used to test the observer.
4.1
Introduction
In automatic motion control of marine vessels we usually consider set-point stabilization.

This means that the vessel is actuated such that a defined point P0 in the vessel body
frame converges to and tracks a predefined position or path in the NED frame. As it is
often impractical or impossible to place motion sensing equipment in P0 , it is important to
know accurately where the sensors are placed. This makes it possible to compensate for
the motion difference between the sensor location and P0 . If this is inaccurate it may lead
to degraded position tracking performance and increased fuel consumption. Today, these
lever arms are measured at dock using a laser surveying technique. Although accurate, it is
time consuming, expensive, and potentially subject to human errors as the laser equipment
must be moved several times for each lever arm. Therefore, this Chapter proposes to
estimate the lever arms of both GNSS antennas and MRUs numerically using the sensor
measurements from a maneuver with sufficient rotation and movement of the vessel. This
is attractive as it removes the human influence and does not requiring expensive time at
dock.
For GPS and INS integration extensive research has been performed in Tang et al. (2009),
46
4.1 Introduction
Hong et al. (2005), Hong et al. (2006), Hong et al. (2004), a single antenna GPS with
accurate measurements is used in combination with a low-grade inertial measurement unit
(IMU). In Hong et al. (2005) the observability of the error states in the INS/GPS integrated
system is studied, and in Hong et al. (2006) experimental studies verify the lever arm estimation of the GPS antenna using a setup on a car. In Batista et al. (2011) observability of
linear motion quantities in navigation systems is considered, and sufficient and necessary
conditions for observability are derived.
4.1.1
Scope
The scope of work is to analyze and test the two following methods for simultaneously
estimating multiple lever arms for GNSS sensors:
A Luenberger observer design
An adaptive observer design
These methods will be analytically derived and their stability properties will be analyzed.
As a part of this, observability and persistence of excitation investigations will be performed to investigate what level of motion and perturbation is required for the lever arm
estimates to converge to the true values.
The observability of the system is considered in Section 4.2.1, and requirement for the
problem to be observable is found. A design for a (Luenberger) observer is proposed, and
stability is proven (Section 4.2.1). In Section 4.2.2 an adaptive observer is proposed, and
stability and convergence of the adaptive observer is shown in Section 4.2.2. Further, the
persistence of excitation criteria is investigated in Section 4.2.2.
Also, in Section 4.3 the problem setup for MRU lever arms are formulated, and in Section
4.3.1 the observability of the dynamics is considered. In Section 4.3.2 an observer is
proposed.
For the GNSS lever arms two case studies with data from the sea trials with R/V Gunnerus
are performed in Section 4.4. In the first case study (Section 4.4.1) data from a sea trial
is used, but the GPS data is simulated, and in the second 4.4.1 experimental GPS data is
used. A case study for the MRU lever arms is performed in Section 4.4.2. Data from a sea
trial with R/V Gunnerus is used.
In the following, the rotation matrix will be viewed as a signal, such that the system equations, instead of being nonlinear, can be treated as linear time varying (LTV). The GNSS
and MRU problems are treated separately, and it is assumed that all needed measurements
are available.
47
Chapter 4. Lever arm estimation for GNSS and MRU sensors
4.2
GNSS
As discussed in Section 4.1 P0 is normally a pre-defined point on the vessel. From this
point the lever arms to the different GNSS antennas are measured by surveying with laser
equipment. Then the measurements from the GNSS antennas are used to find the NED
position of P0 as the vessel is in operation.
When the GNSS measurements are used in a estimation algorithm to find the lever arms,
the P0 position that is found is the NED-position of the rotation point of the vessel. Depending on the operation, or load condition this rotation point would move.
4.2.1
Luenberger observer design
The GNSS positions in the North-East-Down (NED) frame is given as

>
PGN SSi = PGN SSiN PGN SSiE PGN SSiD R3 ,
(4.1)
where i denotes a GNSS sensor. For m GNSS sensors this can be written
PGN SS1 (t) = P0 (t) + R()l1
PGN SS2 (t) = P0 (t) + R()l2
..
.
(4.2)
PGN SSm (t) = P0 (t) + R()lm

where

lxi
li = lyi ,
lzi
i = 1 . . . m,
(4.3)
is a vector that contains the body fixed coordinates of a lever arm, and P0 (t) is the NEDposition of the vessel rotation point as discussed in Section 4.1. The rotation matrix R()
between the BODY and the NED frame given by Eq. (2.4) of Section 2.1.1.
Eq. (4.2) can be written as

P0
l1

x = l2 R3n , n = m + 1,
..
.
lm
PGN SS1
PGN SS2
y=
R3m ,
..
.
PGN SSm
48
(4.4)
(4.5)
4.2 GNSS
where x consists of the states to be estimated, and y is the output.

As mentioned in Section 4.1.1 the Euler angles (, , and ) Fossen (2011) are viewed
as an input signals to the system, so it is a linear time varying (LTV) system, instead of
a nonlinear system, setting R() := R(t). Also, p, q, r, and p,
q,
and r are considered
input signals.
The dynamics of P0 can be represented as (Fossen, 2011)
P 0 = R(t)(t),
(4.6)

>
where = u v w is the linear velocity of the vessel in BODY-coordinates, assumed
measured. The lever arms are constants, and hence
li = 0, i = 1 . . . m.
(4.7)
The dynamic lever arm system can be written as

x = Ax + B(t)(t),
(4.8)
y = C(t)x.
(4.9)
The matrices of (4.8) and (4.9) are given by

A = 03n3n ,

B(t) = R(t)> 033
(4.10)
033
I33
I33
C(t) =
.
..
R(t)
033
033
..
.
R(t)
..
.
I33
033
..
.
..
.
033
033
>
3n3
033
..
.
R3m3n .
033
R(t)
(4.11)
(4.12)
Observability assessment
In this section an observability criterion for the system given by Eq. (4.8) - (4.9) will be
investigated. Because the A matrix of (4.8) is zero, the state transition matrix is identity.
This gives the following observability gramian for the system as (Chen, 2009),
Z
t1
W0 (t0 , t1 ) =
C( )> C( )d,
(4.13)
t0
and if W0 (t0 , t1 ) is nonsingular, the system is observable (Chen, 2009).

49
C(t)> for (4.8) - (4.9) is
I33
R(t)>
C(t)> =
033
.
..
033
and C(t)> C(t) is given as
(n 1)I33
R(t)>
..
.
C(t)> C(t) =
R(t)>
R(t)>
..
.
I33
033
R(t)>
..
.
..
I33
033
..
.
033
033
R(t)>
R(t)
I33
033
..
.
R(t)
..
.
033
..
.
..
033
I33
033
R3n3m ,
R(t)
033
..
3n3n
.
.
R
033
I33
(4.14)
(4.15)
From (4.12) it can be found that if R(t) has full rank, then C(t) and C(t)> also have full
rank, i.e. 3m. From Chen (2009) it is found that the rank of C(t)> C(t) R3n3n will at
most be 3m, as
rank(AB) min(rank(A), rank(B))
(4.16)
where both A and B are matrices. In other words the matrix is not full rank in the first
place, and need to add rank as time increases (Anderson et al., 1986). To find the observability condition(s), Theorem 6.O12 from Chen (2009) is used.
Theorem 1 (Thm 6.O12 (Chen, 2009)). Let A(t) and C(t) be continously differentiable,
then the n-dimensional pair (A(t), C(t)) is observable at t0 if there exists a finite t1 > t0
such that
N0
N1
rank . = n
..
Nn1
where Nm+1 = Nm (t)A(t) +
with N0 = C(t).
d
dt Nm (t)
m = 0, 1 . . . n 1
For the system of (4.8) - (4.9) both C(t) and A(t) are continuously differentiable, and
C(t) has rank 3m given that R(t) has full rank. Hence, given C(t) full row rank, a rank
of 3 is lacking to fulfil rank n of theorem 6.O12 (C(t) has three more columns than rows).
Thus, N1 , ...Nn1 must contribute with 3 more independent rows for the observability
condition to be satisfied.
N1 :
N1 = N0 (t)A(t) +
= 033 +
50
d
N0 (t)
dt
d
C(t),
dt
(4.17)
(4.18)
4.2 GNSS
033
0
C = 33
..
.
033
R(t)S(t)
033
033
..
.
033
R(t)S(t)
..
.
..
.
..
.
033
033
..
.
,
033
R(t)S(t)
(4.19)
where it has been used that R(t)

= R(t)S(t) from Eq. 2.5, and S(t) is given by Eq.
(2.6).
N2 :
N2 = N1 (t)A(t) +
d
N1 (t)
dt
(4.20)
= C(t),
033
= 033
C
.
..
033
(4.21)
R(t)[S(t)2 + S(t)]
033
033
..
.
R(t)[S(t)2 + S(t)]
..
.
..
.
..
.
033
033
033
..
.
033
R(t)[S(t)2 + S(t)]
(4.22)
giving
I33
I33
..
.
I33
0
33
" #
0
N0
33
N1 = .
..
N2
033
033
033
..
.
033
R(t)
033
033
..
.
033
R(t)S(t)
R(t)
..
.
033
033
..
.
033
R(t)[S(t)2 + S(t)]
R(t)S(t)
..
.
033
033
..
.
R(t)[S(t)2 + S(t)]
..
.
033
..
.
..
.
033
..
.
..
.
033
..
.
..
.
033
033
..
.
033
R(t)
033
..
033
R(t)S(t)
033
..
033
R(t)[S(t)2 + S(t)]
(4.23)
where
0
S = r
q
r
0
p
q
p .
0
(4.24)
For the system to be observable, the requirement is that there exists a time t1 such that

N0
rank(N1 ) = 3n.
N2
51
Since S(t) is a skew-symmetric matrix, it is always singular. The observability require

ment is therefore that there has to exists a time t1 where rank (R(t)[S(t)2 + S(t)])
= 3.
Since det R(t) = 1, R(t) is always nonsingular (Chen, 2009), and the observability re
quirement reduce to that [S(t)2 + S(t)]
need to have full rank. This is summarized in the
proposition below.
Proposition 1. For the system of (4.8) and (4.9) to be observable at time t0 , there have to
exist a time t1 > t0 where
rank[S(t)2 + S(t)]
= 3.
(4.25)
Example 1. For a constant yaw rate r = 1, p, q, r,

q = 0, and p = 0.1 at t1 , S(t)2 + S(t)
becomes
1 0
0
S(t)2 + S(t)
= 0 1 0.1
0 0.1
0
which has full rank.
By trying to maneuver the vessel with a constant yaw rate r, observability will be assured
if there also exists some acceleration in roll or pitch. This will in practice always be the
case.
Luenberger observer design and stability analysis
For the system dynamics of equations (4.8), and (4.9) the following observer is proposed
x
= B(t)(t) + W C(t)> y
= B(t)v(t) + W C(t)> C(t)

x,
y
= C(t)
x,
(4.26)
(4.27)
where W = W > > 0 R3n3n , x

= xx
, y
= yy
, and with closed loop error
dynamics as
x
= W C(t)> C(t)
x.
(4.28)
In order to evaluate the stability properties of the observer design, the following lemma
from Anderson et al. (1986) is applied.
Lemma 1. Exponential stability of LTV system (Anderson et al., 1986)
For a system given by x = F (t)x, and the function F () is locally integrable. Suppose
there exits a positive definite matrix P = P > > 0 such that
P F (t) + F (t)> P N (t)> N (t)
(4.29)
for some matrix function N () and all t. Then x = F (t)x is uniformly stable in the sense
of Lyapunov.
52
4.2 GNSS
If, further the pair [F (t), N (t)] is uniformly completely observable, that is, writing (t, )
as the transition function of x = F (t)x, there exits T > 0, > 0, > 0 such that
Z t+T
(t, )> C( )> C( )(t, )d I
I
t
then x = F (t)x is exponentially stable.

From equations (4.27) and (4.28), F (t) = W C(t)> C(t) and P = W 1 gives
P F (t) + F (t)> P = 2C(t)> C(t)
= N (t)> N (t),
(4.30)
where
N (t) =
2C(t),
(4.31)
and the discussion above is summarized in Theorem 2 below.

Theorem 2. The observer design of (4.26) is exponentially stable given that the observability condition of Proposition 1 holds, according to Lemma 1.
4.2.2
Adaptive observer design
Another way to solve the lever arm estimation problem is an adaptive solution. This set
up will remove P0 as a variable to be estimated, and thus have full state measurements
available. The value of P0 can then be found once the lever arms have converged, from
Eq. (4.2) as
P0 = PGN SSi (t) R()li .
(4.32)
By taking the time derivative of (4.2), using (4.6) and R() = R(t), S() = S(t), the
lever arm problem can be formulated as
PGN SS1 (t) = R(t)(t) + R(t)S(t)l1
PGN SS2 (t) = R(t)(t) + R(t)S(t)l2
(4.33)
..
.
PGN SSn (t) = R(t)(t) + R(t)S(t)ln ,
and Eq. 4.33 can be written as
PGN SS1
PGN SS2
x(t) =
,
..
(4.34)
PGN SSn
y(t) = I3n3n x,
(4.35)
53
where x are the states to be estimated, and y is the output. This system can be written as
x = B(t)(t) + (t),
(4.36)
y = Cx,
(4.37)
where
R(t)
R(t)
B(t) = . R3n3 ,
..
R(t)
C = I3n3n ,

l1
l2
3n
=
R ,
ln
R(t)S(t)
033
..
033
.
R(t)S(t)
(t) =
..
..
.
.
033
033
(4.38)
(4.39)
(4.40)
033
..
.
R3n3n .
033
R(t)S(t)
(4.41)
Adaptive observer design and stability analysis

Let a state observer be given as
x
= B(t)(t) + (t)
+ Ly LC x
= B(t)(t) + (t)
+ LC x
,
(4.42)
=xx
. Let
= ,
and
where L R3n3n and x
x
= x x
= LC x
+ ,
=
= .
(4.43)
(4.44)
Define the following Control Lyapunov Function (CLF) (Khalil, 2002)

V =
1 >
1 > 1
x
x
+
,
2
2
(4.45)
where the constant matrix = > > 0. This gives
V = x
> [LC x
+ ]
+
> 1 ()
>
>
> 1
=
x LC x
+x

54
(4.46)
(4.47)
4.2 GNSS
For the following update law for
:= > x
(4.48)
V =
x> LC x
0.
(4.49)
V becomes
Let LC > 0, then V is negative semidefinite. Barbalats lemma (Slotine and Li, 1991) is
applied, and restated here for convenience.
Lemma 2. (Barbalats lemma) If the differentiable function f(t) has a finite limit as
t , and if f is uniformly continuous, then f (t) 0 as t .
To assess uniform continuity, Slotine and Li (1991) states that a sufficient condition for a
differentiable function is that its derivative is bounded. For the above CLF, this implies
that V should be bounded, and V is given as
V = 2
x> LC x
(4.50)
From Eq. (4.49) it is shown that V 0, and therefore V is bounded, and from (4.45) it
can be concluded that both x
and
are bounded. From (4.43), since both x
and
are
bounded, x
is bounded, and hence is V bounded. So V is uniformly continuous, and by
Barbalats lemma V 0 as t , and hence x
0 as t . This implies that
x
0 as t , and from (4.43) it follows that
0 as t .
Percistence of excitation
Consider
= 0. In order to have
= 0 as the only solution need to be persistently
excited. For a update law of the form (4.48), then persistence of excitation (PE) can be
formulated as (Slotine and Li, 1991)
Theorem 3. (Persistence of excitation) The matrix is persistently excited if there exists
, T > 0 such that t
Z t+T
( )> ( )d > I
(4.51)
t
Seeing that (t) is diagonal, condition (4.51) will only depend on R(t)S(t), and PE
criterion can be evaluated based on R(t)S(t) as the integrand (with I := I33 in (4.51)).
Looking at the integrand
(R(t)S(t))> R(t)S(t) = S(t)> R(t)> R(t)S(t)
>
= S(t) S(t).
(4.52)
(4.53)
This result is valid due to R> R = I, which is a fundamental property of the rotation
matrix (Fossen, 2011). The expression for S(t)> S(t) is given as
r(t)2 + q(t)2
p(t)q(t)
p(t)r(t)
r(t)2 + p(t)2
q(t)r(t) ,
S(t)> S(t) = p(t)q(t)
(4.54)
p(t)r(t)
q(t)r(t)
p(t)2 + q(t)2
55
and the PE criteria can be written as

Z t+T
S( )> S( )d =
t
R t+T
R t+T
[q( )2 + r( )2 ]d
t
p( )q( )d
t R
R t+T
t+T
p( )q( )d
[p( )2 + r( )2 ]d
t
t R
R t+T
t+T
t
p( )r( )d
t
q( )r( )d
R t+T
t
p( )r( )d
R t+T
(4.55)
q( )r( )d
R t+T t
[p( )2 + q( )2 ]d
t
> I33 ,
A matrix is positive definite if the leading principal minors are positive (Chen, 2009). This
gives the following three conditions for Eq. (4.55) to be satisfied
1.
Z
t+T
[q( )2 + r( )2 ]d > 0
(4.56)
2.
Z
t+T
2
t+T
[p( )2 + r( )2 ]d
[q( ) + r( ) ]d
t
"Z
#2
t+T
p( )q( )d
>0
(4.57)
3.
(Z
t+T
t+T
[p( )2 + r( )2 ]d
[q( ) + r( ) ]d
t
"Z
t+T
q( )r( )d
#2
(Z
t+T
t+T
t+T
[p( )2 + q( )2 ]d +
t+T
Z
q( )r( )d
p( )r( )d
(Z
t+T
t+T
p( )r( )d
q( )r( )d +
t
t+T
t+T
p( )q( )d
Z
p( )r( )d
p( )q( )d
[p( )2 + q( )2 ]d
p( )q( )d
Z
t+T
t+T
t+T
2
[p( ) + r( ) ]d
>0
(4.58)
Proposition 2 (PE of adaptive observer). The adaptive observer of Eq. (4.42) and (4.48)
is persistently excited if the conditions of (4.56) - (4.58) are satisfied.
56
4.3 MRU
Remark 1. If p = 0 t [t, t + T ], the minimum requirement is that r 6= 0 for some time

t in [t, t + T ], and that
"Z
#2
Z
Z
t+T
t+T
t+T
r( )2 d
q( )2 d >
q( )r( )d
(4.59)
is satisfied. If q = 0 t [t, t + T ], the minimum requirement is that r 6= 0 for some time

t in [t, t + T ], and that
"Z
#2
Z
t+T
t+T
r( )2 d >
p( )r( )d
(4.60)
is satisfied.
Example 2. For a constant yaw rate (steady turn) r, that is nonzero for a time T , assume
pitch motion is zero, and roll motion oscillates with a zero mean value. Let
p(t) = sin(t),
(4.61)
q(t) = 0,
(4.62)
r(t) = cr 6= 0 t [t, t + T ],
(4.63)
where cr is a constant. The left side of Eq. (4.60) from Remark 1 gives
Z t+T
c2r d = c2r T,
t
and the right side gives

"Z
#2
t+T
cr sin( )d
= c2r (cos(t) cos(t + T ))2
4c2r ,
such that for T > 4,
c2r T
"Z
>
#2
t+T
cr sin( )d
= c2r (cos(t) cos(t + T ))2 ,
(4.64)
and PE is satisfied.
The discussion above in concluded in Theorem 4 below.
Theorem 4. The adaptive observer design of (4.42) and (4.48) is globally asymptotically
stable if Proposition 2 is satisfied.
4.3
MRU
A MRU (Seatex, 2014) allows for the measuring of the body fixed accelerations (linear and
angular) at the placement of the MRU, and by differentiating Equation (4.2) for the GNSS
57
equation twice a similar type of equations can be obtained for the acceleration, written as
1,
A1 = A0 + R()S()2 l1 + R()S()l
(4.65)
where A1 is the linear NED-acceleration, and A0 is the linear NED-acceleration at the

reference point. The acceleration from the MRUs will have to be transformed into NEDcoordinates in order to be used in the above equation, which is chosen to be in NEDcoordinates to make it similar to the GNSS equation of (4.2). In the same way as for the
GNSS equations of Eq. (4.2), the equations for several MRUs can be written as
A1 (t) = A0 (t) + R()S()2 l1 + R()S()l
1
A2 (t) = A0 (t) + R()S()2 l2 + R()S()l
2
..
.
(4.66)
Am (t) = A0 (t) + R()S()2 lm + R()S()l

m
In the same manner as with the GNSS problem from Section 4.2 the rotation matrix R,
and the other signals, such as the angular rates, and the linear velocities are treated as input
= S(t).
signals to the system, such that R() = R(t), S( = S(t), and S()
In this
way the system is treated as a time varying linear system.
4.3.1
Luenberger observer design
Eq. (4.66) can be written on a state space form as

A0
l1

x = l2 R3n , n = m + 1,
..
.
(4.67)
lm
AM RU 1
AM RU 2
y=
R3m ,
..
(4.68)
AM RU m
where x are the states to be estimated, and y is the output.
From Eq. (4.6)
P0 = R(t) = V0 (t),
and by (time) differentiating P0 twice A0 can be found as
V 0 = R(t)S(t)(t) + R(t)(t)
= A0 (t),
2
A 0 = (R(t)S(t) + R(t)S(t))(t)
+ R(t)S(t)(t)
+ R(t)
(t).
58
(4.69)
(4.70)
4.3 MRU
As with the GNSS-setup, the system description can be written as

x = Ax + B(t)u(t),
(4.71)
y = C(t)x,
(4.72)
where
A = 03n3n ,
R(t)S(t)2 + R(t)S(t)
R(t)S(t) R(t)
033
033
033
B(t) =
..
.
033
033
033
(t)
,
u(t) = (t)
(t)
(4.73)
R3n3 ,
(4.74)
(4.75)
and
C(t) =
I33
I33
.
..
I33
033
033
..
.
..
.
..
.
..
.
033
033
033
..
.
033
2
R(t)S(t) + R(t)S(t)
R3m3n .
(4.76)
Observability assessment
Observability will be investigated in a similar way as with the GNSS setup of Section
4.2.1, and Theorem 1 will be used. As discussed in Section 4.2.1, since the A(t) matrix
C,
and so on are relevant for the analysis. In fact, it is sufficient
of (4.73) is zero, only C,
to include C, and this is shown by a simple example later (in Example 3).
The time derivative of R(t)S(t)2 + R(t)S(t)

gives
d
+ R(t)S(t)
[R(t)S(t)2 + R(t)S(t)]
= R(t)S(t)3 + R(t)S(t)S(t)
+ 2R(t)S(t)S(t)
dt
=: R(t)G(t),
(4.77)
and C can be written
033 R(t)G(t)
033
..
.
0
033
R(t)G(t)
C(t)
= 33
.
.
.
.
.
..
..
..
.
033
033
033
033
..
.
033
R(t)G(t)
R3m3n .
(4.78)
59
Then, by Theorem 1
h
N0
N1

C
= C
I33
I33
..
.
I33
033
033
.
.
.
033
033
033
..
.
..
.
033
R(t)G(t)
033
033
..
.
033
R(t)G(t)
..
.
..
.
..
.
033
..
.
..
.
033
033
..
.
033
R(t)S(t) + R(t)S(t)
.
033
..
033
R(t)G(t)
(4.79)
By the same argument as in Section 4.2.1 the observability requirement is that there has to
exists a time t1 where rank [G(t)] = 3. This is summarized in the proposition below.
Proposition 3. For the system of Eq. (4.71) and (4.72) to be observable at time t0 , there
have to exist a time t1 > t0 where
+ S(t)]
rank[S(t)3 + S(t)S(t)
+ 2S(t)S(t)
=3
Example 3. For a constant yaw rate r = 1, p = 0.1, and r = 0.1, and p, q, q,
p, q, r = 0
at t1 , and , G(t) (Eq. (4.77) becomes
0.3
1
G(t) = 1 0.3
0.1
0
0.2
0.1
0
which has full rank, so the system is observable at t0 for this excitation.
4.3.2
Observer equations and stability
Given the observer dynamics below

x
= B(t)u(t) + W C(t)> y
>
= B(t)u(t) + W C(t) C(t)

x,
(4.80)
(4.81)
the error dynamics is given as

x
= W C(t)> C(t)
x,
(4.82)
which is similar to the error dynamics of the observer of Section 4.2.1. Therefore, following the same approach as Section 4.2.1, if if Proposition 3 is satisfied, exponential stability
is concluded.
60
4.4 Case studies
4.4
Case studies
Although the theoretical foundation derived above guarantees convergence of the lever
arms for the various designs under the stated assumptions, it does not evaluate performance
and robustness. To shed light on this, case studies are performed for the GNSS and the
MRU lever arm estimation. Data from sea trials with R/V Gunnerus is used for the case
studies.
For the GNSS lever arms two different case studies are performed for both the Luenberger
observer, and the adaptive design. In the first case study the GPS data is simulated, using
sea trial data. That is, the linear velocity, and Euler angles (Fossen, 2011) from the sea
trial are used, and the GPS data is simulated based on that data (see Section 4.4.1). The
reason GPS data is simulated, is to avoid signal synchronization and consistency issues in
the input data. For the second case study the actual GPS data from the experiment is used.
For the MRU lever arm estimation one case study is performed. That is with the Luenberger observer of Section 4.3.2 on experimental data. All observer gains are found by
trial and error.
4.4.1
Case studies - GNSS (GPS) lever arms
Gunnerus data In the August sea trials three turning maneuvers were performed. Two
of them has a constant yaw rate, but with different rates, and the last one has a varying
yaw rate. The intention was to use these maneuvers for lever arm verification. Unfortunately, the GPS measurements from MARINTEK are poor during those maneuvers(low
resolution), and only the DP system GPS measurements are available. Since only one
GPS measurement series is available for these maneuvers, a turning circle performed by
MARINTEK is used instead. This still only gives one GPS measurement series, but all
the other signals used (angular rates and linear velocities) are collected from MARINTEK
measurements, so this reduce the problem of synchronization issues, since all the data is
collected from the same source.
For the GPS case studies, data from a turning circle (Fossen, 2011) maneuver is applied.
The vessel turns with a constant rudder angle of 20 degrees on the rudders. The maneuver was performed by MARINTEK in the November sea trials. The GPS measurements
available are the MARINTEK measurements, and the values for roll, pitch, and yaw rate
from this maneuver are shown in Figure 4.1 below. The roll and pitch rates oscillate about
a zero mean, whereas the yaw rate oscillates around a nonzero mean.
61
0.2
[rad/s]
0.1
0
0.1
0.2
50
100
150
200
Time [s]
250
300
350
400
0.1
[rad/s]
0.05
0
0.05
0.1
50
100
150
200
Time [s]
250
300
350
400
[rad/s]
0.06
0.04
r
0.02
0
50
100
150
200
Time [s]
250
300
350
400
(a)
Figure 4.1: Plot of p, q, and r for GPS case study
Installed GPS antennas With the standard coordinate system used for ship navigation
(Fossen, 2011) (x (pos forward), y (positive sb., z (pos down), the GPS antenna from the
DP system, and the GPS antennas used by the seapath have the following coordinates,
measured by surveillance (Parker, 2013).
Table 4.1: GPS coordinates
Antenna
GPS DP
Seapath fwd
Seapath aft
x [m](pos fwd)
0.202
1.571
-0.924
Coordinates
y [m](pos stb)
-0.611
0.691
0.574
z [m](pos down)
-13.891
-13.520
-13.567
R/V Gunnerus had 3 GPS antennas installed during the sea trials, but only two measurement series are available from the data. One from the DP system, and one from MARINTEK (Seapath).
62
4.4 Case studies
Case study 1 - Simulated GPS data
In the first case study the GPS data for two GPS antennas are simulated. The simulations
are done with experimental data for and , and the GPS data is generated from Eq. (4.2)
PGP Si = P0 + R(t)li ,
(4.83)
where P0 is found from integration of (Eq. (4.6))

P0 = R(t).
(4.84)
Lever arms with coordinates as seapath fwd, and seapath aft from Table 4.1 are added
in Eq. (4.83). The initial condition of the estimated arm coordinates are

>

>
1.0m 0.3m 13.0m , and 0.5m 0.9m 14.0m . Simulation results for
the Luenberger, and the adaptive observer are shown below.
Luenberger observer results The lever arm convergence results are shown below in
Figure 4.2 and 4.3, and Table 4.2 summerize the results. The estimation of P0 is shown
in Figure 4.4 where a North-East plot is shown. The estimated P0 -values are initialized at
the correct value of P0 , so the estimated and measured P0 values are very similar during
the entire simulation.
The gains used in the simulation are
W1
W = 033
033
033
W2
033
033
033 ,
W2
(4.85)
where
W1 = I33 ,
1.8 0
W2 = 0 1.8
0
0
(4.86)
0
0 .
1.5
(4.87)
63
2
1.573
Pos [m]
1.572
1.57
20
40
60
80
100
Time [s]
120
140
160
180
0.693
1 xd
l1 y
l1 y d
l1 z
l1 z d
Pos [m]
Pos [m]
l1 x
l
l 1 x
l1 xd
1.571
0.692
l 1 y
l1 y d
0.691
8
0.69
10
20
40
60
80
100
Time [s]
120
140
160
180
Pos [m]
13.519
12
13.5195
l 1 z
l1 z d
13.52
13.5205
14
0.5
1.5
Time [s]
2.5
13.521
20
40
60
(a)
80
100
Time [s]
120
140
160
180
(b)
Figure 4.2: Lever arm coordinates, l1 , observer, simulated GPS data
2
0.922
Pos [m]
0.923
0.925
2 xd
20
40
60
80
100
Time [s]
120
140
160
180
0.576
Pos [m]
Pos [m]
l2 x
l2 xd
l2 y
l2 y d
l2 z
l2 z d
l 2 x
l
0.924
0.575
l 2 y
l2 y d
0.574
8
0.573
20
40
60
80
100
Time [s]
120
140
160
180
Pos [m]
10
12
13.5665
l 2 z
l
13.567
2zd
13.5675
14
0.5
1.5
Time [s]
(a)
2.5
13.568
20
40
60
80
100
Time [s]
120
140
160
180
(b)
Figure 4.3: Lever arm coordinates, l2 , observer, simulated GPS data
Table 4.2: Results, lever arm - observer, simulated GPS data

Lever arm
1
2
64
x [m]
Avg
Std. [104 ]
1.5712
5.6065
-0.9238
5.6065
Lever arm coordinates

y [m]
z [m]
Avg
Std. [104 ]
Avg
Std. [105 ]
0.6912
5.8248
-13.5194
1.5667
0.5742
5.8248
-13.5664
1.5667
Collected from interval
20 - 180 [s]
20 - 180 [s]
4.4 Case studies
40
20
North [m]
20
P 0, m e a s
P 0, e s t
40
60
80
100
40
30
20
10
10
20
30
East [m]
(a)
Figure 4.4: NE plot of P0 , observer
Adaptive observer results The gains used by the adaptive observer are
= diag{45, 35, 80, 45, 35, 80} 1000,
(4.88)
L = 1000I66 .
(4.89)
The initialization of the lever arms are similar for the Luenberger, and adaptive observer,
but the adaptive observer is tuned quite high such that it has a large deviation at the beginning. As seen from Figure 4.5 and 4.6, and Table 4.3 lever arms converge quite well, but
with higher standard deviations than the Luenberger observer.
65
20
2.2
Pos [m]
15
2
l 1 x
l1 xd
1.8
1.6
10
1.4
40
60
80
100
120
140
160
180
Time [s]
l1 x
l1 xd
l1 y
l1 y d
l1 z
l1 z d
2
Pos [m]
Pos [m]
1.5
l 1 y
l1 y d
1
0.5
40
60
80
100
120
140
160
180
Time [s]
Pos [m]
13.45
10
13.5
l 1 z
l1 z d
13.55
13.6
15
20
40
60
80
100
Time [s]
120
140
160
13.65
40
180
60
80
100
120
140
160
180
Time [s]
(a)
(b)
Figure 4.5: Lever arm coordinates, l1 , adaptive, simulated GPS data
25
0.2
Pos [m]
20
15
0.4
0.8
1
40
10
60
80
100
120
140
160
180
Time [s]
2
l2 x
l2 xd
l2 y
l2 y d
l2 z
l2 z d
Pos [m]
Pos [m]
l 2 x
l2 xd
0.6
1.5
0.5
40
l 2 y
l2 y d
60
80
100
120
140
160
180
Time [s]
13.4
Pos [m]
10
15
20
20
40
60
80
100
Time [s]
120
140
160
180
13.6
l 2 z
l2 z d
13.8
14
40
60
80
100
120
140
160
180
Time [s]
(a)
(b)
Figure 4.6: Lever arm coordinates, l2 , adaptive, simulated GPS data
Table 4.3: Results, lever arm - adaptive, simulated GPS data

Lever arm

y [m]
z [m]
Std.
Avg
Std.
Avg
Std.
0.0137 0.6842 0.0108 -13.5169 0.0205
0.0143 0.559 0.0116 -13.5636 0.0206
x [m]
1
2
66
Avg
1.5787
-0.9132
150 - 180 [s]

150 - 180 [s]
4.4 Case studies
40
20
North [m]
20
P 0, m e a s
P 0, e s t
40
60
80
100
40
30
20
10
0
East [m]
10
20
30
40
(a)
Figure 4.7: NE plot of P0 , adaptive observer
Case study 2 - Experimental GPS data

For this case study experimental GPS data is used. As mentioned in the beginning of
Section 4.4.1 only one GPS measurement series is available, so only one lever arm will
be estimated. The GPS data is expected to be translated to a chosen point in the vessel
(a chosen P0 ), such that at least the lever arm coordinate in y and z-direction should
converge to zero, but not necessarily in x-direction, if the chosen P0 does not coincide
with the rotation point of the vessel.
Luenberger observer results The initial conditions of the lever arm estimate are set as

>
3.0m 2.0m 3.0m .

W =
W1
033

033
,
W2
(4.90)
where
W1 = 0.003I33 ,
(4.91)
W2 = diag{0.012, 0.0084, 0.012}.
(4.92)
67
The lever arm convergence is shown in Figure 4.8, and values for x and z-coordinate
averages are given in Table 4.4.
3
Pos [m]
1 xd
50
100
150
200
Time [s]
250
300
350
400
2
Pos [m]
Pos [m]
l1 x
l1 xd
l1 y
l1 y d
l1 z
l1 z d
l 1 x
l
l 1 y
l1 y d
50
100
150
200
Time [s]
250
300
350
400
Pos [m]
50
100
150
200
Time [s]
250
300
350
400
l 1 z
l1 z d
2
3
50
100
(a)
150
200
Time [s]
250
300
350
400
(b)
Figure 4.8: Lever arm coordinates, l1 , observer, experimental GPS data
Table 4.4: Results, lever arm - observer, experimental GPS data

Lever arm

x [m]
y [m]
z [m]
Avg Std. Avg Std.
Avg
Std.
1.17 0.26
-0.05 0.05
130-360 [s] (x), 200-380 [s] (z)
The y-coordinate values are omitted in Table 4.4. The y-coordinate oscillates about zero,
but it is not very clear. Higher tuning of the observer give more noise, but it is reasonable
to assume that the y-value it should converge to is close to zero. Also, for the x-coordinate,
it is difficult to precisely know the value it should converge to, but it appears to be close to
1.
Adaptive observer results The initial conditions of the lever arm estimate are set as

>
3.0m 2.0m 1.0m .
The gains used by the adaptive observer are
= diag{1, 1.2, 15},
(4.93)
L = diag{0.1, 0.1, 0.01}.
(4.94)
The lever arm convergence is shown in Figure 4.9, and values for x and y-coordinate
averages are given in Table 4.5.
68
4.4 Case studies
3.5
3
Pos [m]
2.5
0
50
l1 x
l1 xd
100
150
200
250
300
Time [s]
l 1 x
l1 xd
l 1 y
l1 y d
l 1 z
l1 z d
0.5
1.5
Pos [m]
Pos [m]
1.5
1
l 1 y
l1 y d
0.5
0
0.5
50
100
150
200
250
300
Time [s]
0.5
Pos [m]
0.5
l1 z
l1 z d
0.5
1
1.5
50
100
150
Time [s]
200
250
300
1.5
50
100
(a)
150
200
250
300
Time [s]
(b)
Figure 4.9: Lever arm coordinates, l1 , adaptive, experimental GPS data
Table 4.5: Results, lever arm - adaptive, experimental GPS data
Lever arm

x [m]
y [m]
z [m]
Avg Std.
Avg
Std. Avg Std.
1.17 0.05 -0.003 0.02
170-300 [s] (x), 280-300 [s] (y)
The matrix need high tuning in y and z to have convergence within the time of the maneuver performed. Convergence is the z-coordinate is not good. As with the Luenberger
Observer result the x-coodinate seemingly converges to a value around 1.
The y-coordinate seems to converge, from Figure 4.9, but after 300s, the value for y goes
away from zero, possibly due to low excitation in the interval after 300 seconds.
Simulated vs recorded GPS data

The plot in Figure 4.10 and show the north and east data for simulated, and recorded GPS
data. The data is simulated without any added lever arm. That is, the simulated GPS
data is the simulated P0 . The plots are quite similar, but at some parts they deviate a lot,
especially around 300-350s for north.
The x-coordinate of the lever arm in case study 2 on recorded GPS data (Section 4.4.1)
indicate that the x-coordinate could be close to 1. Because of this, a second plot is shown
in Figure 4.11, where the simulated GPS data now have a lever arm of 1m in x-direction,
and zero in y- and z-direction. Still, the differences of the data are present.
69
40
20
Pos, N [m]
0
20
r ec
sim
40
60
80
100
50
100
150
200
Time [s]
250
300
350
400
150
Pos, E [m]
100
r ec
sim
50
50
50
100
150
200
Time [s]
250
300
350
400
Figure 4.10: Comparison between recorded and simulated GPS data, no added arm
70
4.4 Case studies
40
20
Pos, N [m]
0
20
r ec
sim
40
60
80
100
50
100
150
200
Time [s]
250
300
350
400
150
Pos, E [m]
100
r ec
sim
50
50
50
100
150
200
Time [s]
250
300
350
400
Figure 4.11: Comparison between recorded and simulated GPS data, added arm
4.4.2
Case study: MRU lever arms - Gunnerus data
The data used for the MRU case study is different from the one used for the GPS case
studies. The reason for this is the update rate of the MARINTEK data. In August MARINTEK logs with an update rate of 100Hz, and the MRUs also operate at 100Hz, whereas
in November MARINTEK logs with 50Hz. The GPS has an update rate of 1Hz, so it does
not matter whether data from November or August is used, but for the MRUs a higher update rate is desirable. The data series that is used is one of the maneuvers intended for
lever arm estimation. However, this maneuver is shorter than the one used for the GPS
case studies.
The values for p, q, and r are shown in Figure 4.12. The yaw rate r appears to be somewhat
constant, but it is higher than for the GPS case study of Section 4.4.1.
71
x 10
[rad/s]
2
0
2
4
50
100
150
Time [s]
3
[rad/s]
x 10
0
q
2
4
50
100
150
Time [s]
0.15
[rad/s]
0.1
0.05
0
0.05
50
100
150
Time [s]
(a)
Figure 4.12: Plot of p, q, and r for MRU case study
Experimental measurements from MRUs in Gunnerus are used, and the acceleration measurements are transformed to NED-frame. There are no measurements of ,
and available. The values for are found from Eq. (4.69) as
= R> A0 S,
(4.95)
and the value of

will be neglected from the observer equations in the case study.
Installed MRUs
There were 4 MRUs installed in R/V Gunnerus at the time of the sea trials, and in the case
study data from two of them will be used. The lever arms of those two MRUs are given
in Table (4.6) below. measured by surveillance (Parker, 2013).
Table 4.6: MRU coordinates
Antenna
MRU 1
MRU 2
72
x [m](pos fwd)
0.358
14.978
Coordinates
y [m](pos stb)
0.804
0.039
z [m](pos down)
4.321
0.568
Luenberger observer results

MRU accelerations from MRUs with lever arms of coordinates

>
0.358m 0.804m 4.321m ,
and

14.978m
>
0.568m ,
0.039m

are used. The initial condition of the estimated arm coordinates are 30m

>
and 30m 30m 30m .
30m
30.0m

033
,
W2
W1
033
W =
(4.96)
where
W1 = I33
(4.97)
W2 = diag{1, 0.8, 2} 40000.
(4.98)
The results of the observer are shown in Figure 4.13.
30
20
l 1 x
l1 xd
10
Pos [m]
Pos [m]
30
0
10
25
l 2 x
l2 xd
20
15
50
100
10
150
50
Time [s]
l 1 y
l1 y d
Pos [m]
Pos [m]
20
20
l 2 y
l2 y d
10
0
50
100
10
150
50
Time [s]
150
30
10
l 1 z
l1 z d
20
Pos [m]
Pos [m]
100
Time [s]
30
150
30
10
30
100
Time [s]
10
20
l 2 z
l2 z d
10
0
50
100
150
10
Time [s]
(a)
50
100
150
Time [s]
(b)
Figure 4.13: Lever arm coordinates l1 (a), and l2 (b) MRU case study
4.5
For the GPS case studies where GPS data is simulated it can be seen from Section 4.4.1
that the lever arms and P0 are clearly observable (and the adaptive observer is PE). Both
73
>
observers have convergence in the lever arms, and the maneuver is sufficient for observability. The z-coordinate in the Luenberger observer have a (small) steady state offset, but
this coordinate seems to converge for the adaptive observer, although it is a bit noisy.
For the GPS case studies where the recorded GPS data is used both for the Luenberger and
the adaptive observer the lever arms converge close to the expected value. This becomes
clear when the lever arm estimates are initialized far from the expected value. However,
where the estimates of the lever arms were very accurate for the simulated GPS data, they
are quite inaccurate for the experimental GPS data.
This inaccuracy could be explained by the difference in the simulated, and the experimental GPS data, as shown in Figure 4.10 and 4.11. There are small deviations between the
GPS values, and this could be enough for the algorithms to become less accurate.
From Figure 4.13 is observed that the MRU estimation problem is observable. The lever
arms converge quite well, and is better than the GPS lever arms when using experimental
data. This is probably due to the fact that the MRUs have an update rate of 100Hz,
whereas the GPS has an update rate of 1Hz.
For both the GPS and MRU case studies it would be interesting with a longer maneuver,
such that both observer algorithms could have time to converge with lower tuning. It could
also be interesting to look into the differences between the simulated GPS data, and the
actual.
For the GPS case study on experimental GPS data the x-coordinate of the lever arm seemingly converges to a value close to 1m for both the Luenberger, and the adaptive observer.
It is difficult to say whether this is the case, because these results change slightly depending on the tuning of the observers. This again motivates a longer maneuver such that low
tuning can be used on both the observer designs, and more certain results can be obtained.
74
Chapter
Hybrid integral action for DP

5.1
5.1.1
Introduction
Motivation
The motivation behind the hybrid integral action proposed in the following chapter is a
Dynamically Positioned (DP) marine vessel experiencing large unknown disturbances. DP
vessels normally experience wave loads, wind loads, and current. The loads that the integral action part of the controller normally compensate for are slowly varying forces, almost
constant for given periods of time. Because of this, the integral action is normally tuned
very low, such that it does not induce unnecessary oscillation.
However, if the vessel experiences sudden load changes such as ice loads, or a mooring
line that snaps, the integral action spends a long time (about 20 minutes) reaching the new
steady state value. This chapter proposes a method that improves the convergence for the
integral action when it is subject to large sudden disturbances that are constant some time
after impact. The proposed method augments the standard PID controller with a hybrid
integral action law.
5.1.2
Literature review
In El Rifai and El Rifai (2009), hybrid resetting of integral action for a PID controller is
discussed. Based on the sign of the position state and the integral state, the integral action
value is reset if they are of opposite sign, thereby reducing transient behavior of the closed
loop plant. For Prieur et al. (2012) a hybrid high-gain observer is constructed to reduce the
peaking behavior of the observer on a second order planar nonlinear system. Trajectories
with peaking are projected into areas without peaking behavior.
75
Chapter 5. Hybrid integral action for DP
A framework with several continuous controller and observer-pairs are proposed for hybrid
control of DP vessels in Nguyen et al. (2007). The operational window of a DP vessel is
extended by switching to different observer-controller pairs depending on the sea state.
This approach could have been used to augment the plant with a controller with a high
gain for the continuous integral action, and then switching back to a controller with lower
integral action gains when the integral error is small.
The method of El Rifai and El Rifai (2009) is similar to that of the following chapter in
that it uses the sign of the integral value and the states to determine when jumps can occur.
However, the goal of the following chapter is to use information about the states to update
the integral value, not only reset the integral value when the signs changes. This will be
especially useful when large constant disturbances should quickly be compensated for by
integral action.
5.1.3
Scope
The objective of the chapter is to improve performance of the PID controller when a system is subject to large disturbance changes that remain constant for some time (step disturbance). A PID-controller is augmented with a hybrid (Goebel et al., 2012) integral action
law that changes the integral action value at discrete instances (jumps). When the absolute
value of the error in states are small, jumps are no longer allowed. This discrete change
in integral action value allows higher convergene of the integral action, with no, or small
overshoot. This will be developed for both a first order linear system, and a DP system.
Section 5.2 considers the mathematical modeling. Section 5.2.1 the preliminaries for hybrid control theory and the Lyapunov stability theory needed is summarized. The stability
conditions is then derived for both the first order system (Section 5.2.2), and the DP system
(Section 5.2.3). In these sections a theorem concludes the stability conditions for the hybrid system. In Section 5.3 there are case studies for both a linear system (Section 5.3.1),
and a DP system (Section 5.3.2).
5.2
5.2.1
Mathematical modelling
Preliminaries
The theory for hybrid control theory is based on Goebel et al. (2012). The benefit of the
theory is that continuous and discrete dynamics can be combined, and stability can be
proven. The aspects relevant for the approach of the chapter is mentioned below, but for a
more in-depth analysis of hybrid control theory, it is referred to Goebel et al. (2012).
Continous dynamics, here called flow, given generally by a differential inclusion F (x)
is allowed on the flow set C. The discrete dynamics, here called jumps, given by the
difference inclusion G(x) is allowed on the jump set D. In the following only a differen76
5.2 Mathematical modelling
tial equation f (x), and difference equation g(x) will be used instead of F (x), and G(x)
respectively.
In order to prove stability of the systems in this paper, Theorem 3.18 of Goebel et al.
(2012), using a Lyapunov function is applied. This theorem proves global uniform (pre)asymptotic stability of the system, and requires the Lyapunov function to decrease in
value for both flow and for a jump. Since the jump augmentation of the PID controller
intends to improve performance, it makes sense to demand that V (x) also decreases in
jumps. To prove stability for a set A Theorem 3.18 in its most basic form is applied,
and restated below for convenience. Please note that pre-asymptotically stable includes
solutions not being complete (Goebel et al., 2012), and since this is of no concern for the
systems considered, pre-asymptotically stable is the same as asymptotically stable for
this chapter.
Theorem 5 (Goebel et al. (2012) Theorem 3.18). (Sufficient Lyapunov conditions)
Let H = (C, F, D, G) be a hybrid system and let A Rn be closed. If V is a Lyapunov
function candidate for H and there exists 1 , 2 K , and a continuous positive definite
function such that
1 (|x|A ) V (x) 2 (|x|A ) x C D G(D)
(5.1)
hV (x), f i (|x|A ) x C, f F (x)
(5.2)
V (g) V (x) (|x|A ) x D, g G(D)
(5.3)
then A is uniformly globally pre-asymptotically stable for H.

In the above theorem |x|A is the distance to the set A. In this paper the V notation will be
used instead of hV (x), f i, and for Eq. (??) V (g) and V (x) represents the value of V (x)
after and before a jump respectively. All stability proofs later in the chapter will use this
theorem to prove stability.
About the approach
The key for the proposed hybrid integral action in Section 5.2.2, and 5.2.3 is the location of
the jump set. It will be shown that when within the jump set, the proposed jump rule will
guarantee decrease of the Lyapunov function. The clue is therefore to restrict the jump set
to a set that will depend on the sign and size of the error in the integral value, and the sign
and error in other states. By construction of the problem, the integral error (the difference
between the unknown disturbance and the integral value) is unknown. Therefore, the
dynamic equations will be used to find an expression for the integral error, and a new jump
set will be formulated based on known states, and an estimate of a state derivative (that will
be found by sampling). Flow can occur in the entire state space. The jump rule proposed
is a linear jump rule, and jumps are proportional to the error in other state variables.
5.2.2
First order linear systems
In this approach a general first order system subject to an unknown constant disturbance
is presented. A control law with proportional control, and integral action is proposed, and
77
the closed loop system dynamics is derived. Then the flow set, flow map, and the jump
map are defined, before Theorem 6 defines a jump set such that the jumps always decrease
the Lyapunov function. The theorem also gives conditions for stability in flow, such that
the combined system is stable.
Consider the first order system with an unknown constant disturbance d as input to the
system. The system is written as
z = az + d + u
d = 0,
(5.4)
(5.5)
where a > 0. Let zd be the desired z-value, and by selecting the control input u as
u = zd + azd kp (z zd ) d,
(5.6)
with kp > 0, the closed loop error dynamics becomes

z = (a + kp )
z + d = a0 z + d
d = d d = d = ki z,
(5.7)
(5.8)
a0 = a + kp > 0.
where d = d d,
Below the flow set, flow map, and the jump map are defined. The jump set is defined in
Theorem 6. The states are defined as

z
x
x= = 1 .
(5.9)
x2
d
.
Flow set
The flow set is the entire state space, so C is given as

C = x R2 .
(5.10)
Flow map
or the flow map f (x) is given
From (5.9), (5.7), and (5.8) the time derivative of the state x,
as
" #
0

z
a0 x1 + x2
a 1 x1
x = =
=
= Ax = f (x).
(5.11)
ki x1
ki 0 x2
d
78
Jump map
In the proposed jump map x1 remain the same, and x2 is updated based on the x1 -value,
such that
+

x1
x1
1 0 x1
+
x = + =
=
= x = g(x)
(5.12)
x2 x1
1 x2
x2
Theorem 6. Given a linear system with x R2 of (5.9), the flow set C given by (5.10),
and the closed loop flow map f (x) given by (5.11), and jump map g(x) given by (5.12),
where the constants a0 , ki , > 0, let 1 , 2 K , and a Lyapunov function be given as
1 (|x|A ) V (x) = x> P x 2 (|x|A ),
(5.13)
where
P =

p1
p2

p2
= P > > 0,
p3
(5.14)
is set such that the Lyapunov function decreases in flow. That is,
V (x) (|x|A ) x C.
(5.15)
Let
:= 2p3 > 0,
(5.16)
(5.17)
> 0,
(5.18)
> 0,
(5.19)
:= p3 2p2 > 0,
where and are constants. For the jump set given by

( + )x21 + x22 )
2
,
D = x R : |x1 | , x1 x2
(5.20)
then
V (g(x)) V (x) |x|2 ,
(5.21)
and by Theorem 5 the set A = {0, 0} is uniformly globally pre-asymptotically stable for
H = (C, f, D, g).
Proof. Consider the Lyapunov function given by (5.13). The time derivative of the Lyapunov function gives

V (x) = x> P x + x > P x
(5.22)
>

>
= x (P A + A P )x ,
(5.23)
79
so for
P A + A> P < 0,
V (x) < 0
(5.24)
x C,
(5.25)
so if the matrix P is set such that P A + A> P < 0, the Lyapunov function decrease in
flow.
The value of the Lyapunov function after a jump, V (g) is
V (g) = x> > P x,
(5.26)
such that
V (g) V (x) = x> [> P P ]x,
2
p3 2p2
= x>
p3

p3
x.
0
(5.27)
Expanding (5.27) gives

V (g) V (x) = (2 p3 2p2 )x21 2p3 x1 x2 ,
= x21 x1 x2 ,
(5.28)
and by inserting for (5.20) the inequality of (5.28) becomes

V (g) V (x) |x|2
x D,
(5.29)
so by Theorem (5.3), uniform global pre-asymptotic stability is guaranteed.

For Theorem 6 knowledge of the integral error is used to find the jump set. The integral error is not known, since the integral action is used to compensate for an unknown
disturbance. Therefore, the results of Theorem (6) is not applicable for a practical implementation. In Remark 1, knowledge of the system dynamics is used to estimate the
integral error, and to find a jump set based on x1 and x 1 .
Remark 1 (Practical implementation). For a practical implementation of Theorem 6, the
jump set D can not depend on x2 . From the flow map of (5.11) x1 x2 can be written as
x1 x2 = a0 x21 + x1 x 1 ,
such that the jump set of (5.20) can be rewritten as

( + a0 )x21 + x22 )
2
D = x R : |x1 | , x1 x 1
.
(5.30)
(5.31)
The value of can be set arbitrarily small, such that Eq. (5.31) does in practice not
depend on the value of x2 .
80
5.2.3
DP system
Similar to the first order system, a closed loop system with integral action in the controller
will be derived. Then the flow set, flow map, and jump map are defined. Thereafter
Theorem 7 will define stability of the system by specifying a jump set that depends on the
states of the system.
Consider the linearized DP system of Section 2.1.4, that has kinematics and kinetics given
as
= R(),
(5.32)
>
M = D + R() b + ,
b = 0,
(5.33)
(5.34)
where b is considered a constant disturbance or bias force (from Eq. (5.34 )), and the linear
damping matrix satisfies D > 0, and the mass matrix is assumed to have the following
= 0. The system now contain a rotation matrix, makproperties M = M > > 0, and M
ing it nonlinear. However, the jump map used will be linear, and similar to the approach
used in Section 5.2.2.
By using a backstepping approach with integral action (Skjetne and Fossen, 2004), the
x-variables can be defined as

R>
x1
x = (, t) = x2 ,
(5.35)
x3
b
where (, t) is a virtual control law to be defined later. An integral state
b is augmented
to the plant, and its dynamics are given as
b = Ki R()x2 ,
(5.36)
with Ki = Ki> > 0. The other variables are defined as

= d ,
b = b
b, where
d is the desired position, and the desired velocity is zero (d = 0). By setting the virtual
control law (, t) as
(, t) = Kp x1 + R()> d ,
(5.37)
with Kp = Kp> > 0, and the actual control input as

= x1 Kd x2 R()>
b + D + M ,
(5.38)
where Kd = Kd> > 0. This results in the closed loop continuous dynamics
x 1 = rSx1 Kp x1 + x2
(5.39)
>
M x 2 = x1 Kd x2 + R() x3
x 3 = Ki R()x2 ,
(5.40)
(5.41)
81
and this is shown in Appendix A.

From the states defined in (5.35), and the continuous closed loop dynamics of (5.39 - 5.41),
the flow set, the flow map, and the jump map are defined below. The jump map is defined
similarly to that of Section 5.2.2, and the jump set is defined in Theorem 7.
Flow set
Flow should be allowed in the entire state space, and the flow map C is given as

C = x R9 .
(5.42)
Flow map
Rewriting the equations (5.39 - 5.41), the flow map f (x) is given as
x 1 = rSx1 Kp x1 + x2
(5.43)
>
M x 2 = x1 Kd x2 + R() x3
x 3 = Ki R()x2 .
(5.44)
(5.45)
Jump map
For the jump map chosen the jump in integral action is proportional to x2 . That is, the
same state variable used in the continuous integral action.
For a constant > 0,
= diag{, , } > 0,
the jump map is given as
+

x1
x1
I
x
0
x+ = x+
=
=
2
2
x
x
0
x+
3
2
3

0
0
x1
I
0 x2 = e = g(x)
I
x3
(5.46)
(5.47)
Theorem 7. Given the closed loop DP system with x R9 of (5.35), the flow set given
by (5.42), and the closed loop flow map given by (5.43 - 5.45 ), and jump map given by
= 0, Ki = K > , 1 , 2 K , and a Lyapunov function be
(5.47), let M = M > , M
i
given as
1 (|x|A ) V (x) =
1
1
1 >
x x1 + x>
M x 2 + x>
K 1 x3 2 (|x|A ),
2 1
2 2
2 3 i
such that V (x) < 0, x C can be shown.

82
(5.48)
Let
:= 2P3 > 0,
(5.49)
(5.50)
:= P3 2P2 > 0,
1 > 0,
(5.51)
2 > 0,
(5.52)
> 0,
(5.53)
where 1 , 2 , and are constants. For the jump set given by

D = x R9 : |x1 | 1 , |x2 | 2 ,
1
x>
2 Ki x 3
(5.54)

> 1
>
>
x K x2 +
[x x1 + x>
2 x2 + x3 x3 ] ,
2 2 i
2 1
(5.55)
then
V (g(x)) V (x) |x|2 ,
(5.56)
and the set A = {0, 0, 0} is uniformly globally pre-asymptotically stable for H = (C, f, D, g).
Proof. Stability in flow:
The time derivative of the Lyapunov function given by (5.48) gives
1
V (x) = x>
1 + x>
2 + x>
3
1x
2 Mx
3 Ki x
>
>
= x>
1 [rSx1 Kp x1 + x2 ] + x2 [x1 Kd x2 + R() x3 ]
1
+ x>
3 Ki [Ki R()x2 ]
>
= x>
1 Kp x1 x2 Kd x2 0,
(5.57)
and the continuous dynamics is uniformly globally stable (UGS) (Khalil, 2002), and to
proove UGAS Barbalats Lemma (Lemma (2)) is applied, since the system is time varying.
The double time derivative of V (x), V (x) is
V (x) = 2x>
1 2x>
2.
1 Kp x
2 Kd x
(5.58)
From Eq. (5.48) it is known that x1 , x2 , and x3 are bounded. From Eq. (5.43) is can be
seen that x 1 is bounded, and from Eq. (5.44) x 2 is bounded, and hence, V (x) is bounded,
and
>
lim V = lim (x>
1 Kp x1 x2 Kd x2 ) = 0
Since x2 goes to zero as time goes to infinity, then

lim x 2 = 0,
83
and by Eq. (5.44)

lim x3 = 0,
and UGAS can be concluded for flow x C.

Stability in jumps:
The value of the Lyapunov function after a jump V (g) is given as
1 >
1
1
x x1 + x>
M x2 + (x3 x2 )> Ki1 (x3 x2 )
2 1
2 2
2
1 >
1 >
= x1 x1 + x2 M x 2
2
2
1 > 1
1
1
1
> >
> >
+ [x3 Ki x3 x>
3 Ki x2 x2 Ki x3 + x2 Ki x2 ],
2
V (g) =
(5.59)
such that
1
1
1
> >
> >
V (g) V (x) = x>
3 Ki x2 x2 Ki x3 + x2 Ki x2 ,
(5.60)
and since = I,
1
1
2 >
V (g) V (x) = 2x>
2 Ki x3 + x2 Ki x2 .
(5.61)
For
1
x>
2 Ki x3
>
> 1
>
x2 Ki x2 +
[x x1 + x>
2 x2 + x3 x3 ],
2
2 1
(5.62)
V (g) V (x) |x|2 < 0.
(5.63)
As with the approach of the first order system, can be arbitrary small, so it does not
matter that x3 is unknown.
Since UGAS can be proved for both flow and jumps, UGAS (UGpAS) for the hybrid
system is concluded.
Remark 2 (Practical implementation). For a practical implementation of Theorem 7, the
1
jump set D can not depend on x3 . From the flow map equation of (5.44) x>
2 Ki x3 can
be written as
1
1
>
x>
2 + x1 + Kd x2 ]
2 K i x3 = x2 Ki R()[M x
such that the jump set of (5.55) can be rewritten as

D = x R9 : |x1 | 1 , |x2 | 2 ,
>
1
>
x>
2 x>
[x1 x1 + x>
2 Ki R()M x
2 KI x2 +
2 x2 + x3 x3 ]
2
2

1
x>
2 Ki R()[x1 + Kd x2 ]
(5.64)
(5.65)
(5.66)
(5.67)
As in Remark 1, the value of can be set arbitrarily small, such that Eq. (5.67) does in
practice not depend on the value of the integral state error x3 .
84
5.3 Case studies
5.3
5.3.1
Case studies
First order linear system
For this example both the cases when knowledge of x2 is assumed known, and when x2 is
estimated from sampling of x 1 are simulated. Consider the first order system
z = z + d + u
d = 0,
(5.68)
u = zd + zd 9(z zd ) d,
(5.70)
(5.69)
and control input u as
so that the closed loop continuous dynamics becomes

z = 10
z + d
d = ki x1 .
(5.71)
(5.72)
Let the state vector be given as

z
x
x= = 1 .
x2
d
(5.73)
Typically the integral action would have a low tuning, so set ki = 1, and a0 is given as 10
by (5.71), so x is written as

10 1 x1
x =
= Ax,
(5.74)
1 0 x2
and for

P =

1
0.1
0.1
1
(5.75)
the condition given by (5.24) is satisfied, and the Lyapunov function decrease in flow.
For this example, the values chosen for , and are
= 3,
(5.76)
= 0.3,
(5.77)
such that Eq. (5.16), and (5.17) gives = 6, and = 9.6, and the jump set given by
(5.20) becomes

3 2
(x1 + x22 ) .
(5.78)
D = x R2 : |x1 | 0.3, x1 x2 1.6x21 +
85
Note that 3 can be arbitrary small, so it is not considered in the simulation.
The jump set for when x2 is not assumed known, given by (5.31) becomes

D=
x R : |x1 | 0.3, x1 x 1
8.4x21

3 2
2
+
(x + x2 ) .
1
(5.79)
For both simulations a step of magnitude 100 at time t = 5s, and a step of 50 at time
t = 30 is applied.
Even though it is not necessary in terms of stability, a dwell-time (Goebel et al., 2012)
(minimum amount of time between jumps) of T = 0.1s is demanded between each jump
in the simulations. For the simulation run with x2 assumed known, it would work fine
without a dwell time, but if is set small, then all jumps necessary would be performed
with a magnitude 2 , which could require a lot of jumps. A method without dwell-time
could be computationally unfeasible. Also, a method without a dwell-time constraint could
be less robust, since it allows for an infinite number of jumps in no time, which means that
is could respond too quickly to noise. However, if the dwell-time is slower than the noise,
this would work as a noise filter. For the sampled version a dwell-time is needed in order
to sample values of x1 .
For the second simulation run where x 1 is found by sampling, a sampling period of Ts =
0.03s is applied. Let x01 be the value of x1 at previous sampling instant. Then x 1 is found
from
x1 =
x1 x01
.
Ts
(5.80)
White noise of power 0.1 is added to make the sampling process more realistic. It is not
crucial that the value for x 1 is correct, but it should be close in magnitude to the true value,
at least when close to the boundary of the jump set. This is because x 1 is just used as a
jump condition (defines the jump set), and the value itself is not used in any feedback.
86
5.3 Case studies
Results first order system - with knowledge of x2
Position error. With jumping (blue) and without (green)

15
Position error
10
10
15
20
25
Time [s]
30
35
40
45
50
Integral state error. With jumping (blue) and without (green) red lines indicate jumps
Integral state error
100
50
50
100
10
15
20
25
Time [s]
30
35
40
45
50
Figure 5.1: Position and integral state error first order system, no sampling
Integral state error zoom. With jumping red lines indicate jumps
100
80
60
40
20
0
20
5.2
5.4
5.6
Time [s]
5.8
6.2
6.4
Lyapunov function value for the hybrid solution red lines indicate jumps
Lyapunov function
10000
8000
6000
4000
2000
0
5.2
5.4
5.6
Time [s]
5.8
6.2
6.4
Figure 5.2: Lyapunov function and integral state error zoom at time t = 5s, no sampling
87
Results first order system - without knowledge of x2
Position error. With jumping (blue) and without (green)

15
Position error
10
10
15
20
25
Time [s]
30
35
40
45
50
45
50
Integral state error. With jumping (blue) and without (green)
100
50
50
100
10
15
20
25
Time [s]
30
35
40
Figure 5.3: Position and integral state error first order system, with sampling, and white noise.
Integral state error zoom. With jumping red lines indicate jumps
100
80
60
40
20
0
20
Time [s]
Lyapunov function
10000
8000
6000
4000
2000
0
Time [s]
Figure 5.4: Lyapunov function and integral state error zoom at time t = 5s, with sampling, and
white noise.
88
5.3 Case studies
5.3.2
Example - DP
In this example, knowledge of x3 will be assumed known for the first simulation run. For
the second simulation, x2 will be sampled to find an estimate of x 2 . Some noise will be
added to this simulation, to make the sampling process more realistic.
Consider the closed loop DP system of Eq. (5.43 - 5.45) as
x 1 = rSx1 Kp x1 + x2
(5.81)
>
M x 2 = x1 Kd x2 + R() x3
x 3 = Ki R()x2
(5.82)
(5.83)
with
M = diag{450, 450, 100}
(5.84)
Kp = diag{100, 100, 50}
(5.85)
Ki = diag{10, 10, 5}
(5.86)
Kd = diag{500, 500, 200},
(5.87)
and , 1 , and 2 chosen as

= 25,
(5.88)
1 = 0.0001,
(5.89)
2 = 0.01.
(5.90)
For both simulations a step of magnitude
3000
b = 2000 ,
1000
(5.91)
is applied at time t = 10s.

By the same reasoning as for the example for the first order system of Section 5.3.1 a
dwell-time is added between jumps here as well. The dwell chosen is T = 0.1s.
Similar to the example for the first order system of Section 5.3.1 a dwell-time is needed
in order to sample values of x2 . For the second simulation run where x 2 is found by
sampling, a sampling period of Ts = 0.03s is applied. Let x02 be the value of x2 at
previous sampling instant then x 2 is found from
x 2 =
x2 x02
.
Ts
(5.92)
In the second simulation noise of power 20 has been added in all DOF.
89
Results DP - with knowledge of x3

For the results below is says (N ), (E), (Heading), on the plots, but this is not entirely
the case. The x1 state is
transformed to body coordinates, and x2 is not equal to , but
in the plots x1 , and x2 are called position and velocity.
Position error [m]
Position error [m]
Position error, p1 (N). With jumping (blue) and without (green)

0.1
0.1
10
20
30
40
50
60
Time [s]
Position error, p2 (E). With jumping (blue) and without (green)
70
80
10
20
30
40
50
60
Time [s]
Position error, p (Heading). With jumping (blue) and without (green)
70
80
70
80
20
30
40
50
60
Time [s]
Velocity error, v (E). With jumping (blue) and without (green)
70
80
20
30
40
50
60
Time [s]
Velocity error, v3 (Heading). With jumping (blue) and without (green)
70
80
70
80
0.04
0.02
0
0.02
Position error [rad]
0.1
0.1
10
20
30
40
Time [s]
50
60
Figure 5.5: Position error DP, no sampling
Velocity, v (N). With jumping (blue) and without (green)

Velocity error [m/s]
10
5
0
5
10
Position error [rad/s]
Position error [m]
4
2
0
2
10
10
20
30
40
Time [s]
50
60
Figure 5.6: Velocity error DP, no sampling
90
5.3 Case studies
2000
1000
0
1000
10
10
10
20
30
40
50
60
Time [s]
Integral state error, (E). With jumping (blue) and without (green)
70
80
20
30
40
50
60
70
Time [s]
Integral state error (Heading). With jumping (blue) and without (green)
80
2000
1000
0
1000
Integral state error, (N). With jumping (blue) and without (green)
3000
1000
500
0
500
20
30
40
Time [s]
50
60
70
80
15
16
15
16
Figure 5.7: Integral error DP, no sampling
Integral state error, (N) Zoom. With jumping red lines indicate jumps
3000
2500
2000
1500
1000
500
0
500
Lyapunov function
10
11
12
Time [s]
13
14
x 10
1.5
0.5
10
11
12
Time [s]
13
14
Figure 5.8: Lyapunov function and zoom on integral error, no sampling
91
Results DP - without knowledge of x3 - sampling of x2
Position error [m]
Position error [m]
Position error, p1 (N). With jumping (blue) and without (green)

0.1
0.1
10
20
30
40
50
60
Time [s]
Position error, p2 (E). With jumping (blue) and without (green)
70
80
10
20
30
40
50
60
Time [s]
Position error, p (Heading). With jumping (blue) and without (green)
70
80
70
80
20
30
40
50
60
Time [s]
Velocity error, v (E). With jumping (blue) and without (green)
70
80
20
30
40
50
60
Time [s]
Velocity error, v3 (Heading). With jumping (blue) and without (green)
70
80
70
80
0.04
0.02
0
0.02
Position error [rad]
0.1
0.1
10
20
30
40
Time [s]
50
60
Figure 5.9: Position error DP, with sampling
Velocity, v (N). With jumping (blue) and without (green)

1
Velocity error
10
5
0
5
10
Velocity error
4
2
0
2
10
10
Velocity error
20
30
40
Time [s]
50
60
Figure 5.10: Velocity error DP, with sampling
92
2000
0
2000
10
10
10
20
30
40
50
60
Time [s]
Integral state error, (E). With jumping (blue) and without (green)
70
80
20
30
40
50
60
70
Time [s]
Integral state error (Heading). With jumping (blue) and without (green)
80
3000
2000
1000
0
1000
Integral state error, (N). With jumping (blue) and without (green)
4000
1500
1000
500
0
500
20
30
40
Time [s]
50
60
70
80
Figure 5.11: Integral error DP, with sampling
Integral state error, (N) Zoom. With jumping red lines indicate jumps
4000
3000
2000
1000
0
1000
10
12
14
16
18
16
18
Time [s]
6
Lyapunov function
x 10
1.5
0.5
10
12
14
Time [s]
Figure 5.12: Lyapunov function and zoom on integral error, with sampling
5.4
For both the first order linear system, and the DP system the hybrid integral action improves performance for large constant disturbances. The state estimates, and the integral
action error converge much faster than in the case without jumps. After a step the integral
action converge to the new value quickly, and (almost) without overshoot. When sampling
93
is used to find the jump set, the performance decrease is negligible compared to the ideal
case, and the noise does not induce problems.
For further work the algorithm should be tested on more realistic data, either a comprehensive numerical ice simulation, or a model test for a DP vessel subject to large step
changes in disturbances. It would be interesting to see how this approach could improve
performance (offset from desired position), thrust, and power consumption. One other
interesting aspect for a model test is how performance is affected by sensor noise. For
the DP system an alternative to sampling the velocity is to use accelerometers to find the
accelerations.
94
Bibliography
Anderson, B. D. O., Bitmead, R. R., Johnson, Jr., C. R., Kokotovic, P. V., Kosut, R. L.,
Mareels, I. M., Praly, L., Riedle, B. D., 1986. Stability of Adaptive Systems: Passivity
and Averaging Analysis. MIT Press, Cambridge, MA, USA.
Batista, P., Silvestre, C., Oliveira, P., 2011. On the observability of linear motion quantities
in navigation systems. Systems & Control Letters 60 (2), 101110.
Chen, C.-T., 2009. Linear System Theory and Design - International Edition, 3rd Edition.
Oxford University Press.
El Rifai, K., El Rifai, O., 2009. Design of hybrid resetting pid and lag controllers with
application to motion control. In: Advanced Intelligent Mechatronics, 2009. AIM 2009.
IEEE/ASME International Conference on. IEEE, pp. 685692.
Faltinsen, O., 1990. Sea loads on ships and offshore structures. Cambridge University
Press.
Fossen, T. I., 2011. Handbook of marine craft hydrodynamics and motion control. John
Wiley & Sons.
Goebel, R., Sanfelice, R. G., Teel, A. R., 2012. Hybrid Dynamical Systems: modeling,
stability, and robustness. Princeton University Press.
Hong, S., Lee, M. H., Chun, H.-H., Kwon, S.-H., Speyer, J. L., 2005. Observability of
error states in gps/ins integration. Vehicular Technology, IEEE Transactions on 54 (2),
731743.
Hong, S., Lee, M. H., Chun, H.-H., Kwon, S.-H., Speyer, J. L., 2006. Experimental study
on the estimation of lever arm in gps/ins. Vehicular Technology, IEEE Transactions on
55 (2), 431448.
Hong, S., Lee, M. H., Kwon, S. H., Chun, H. H., 2004. A car test for the estimation
of gps/ins alignment errors. Intelligent Transportation Systems, IEEE Transactions on
5 (3), 208218.
95
Khalil, H. K., 2002. Nonlinear systems, 3rd Edition. Prentice hall Upper Saddle River.
Kongsberg Maritime, 2006. Kongsberg K-Pos Dynamic Positioning - Optimizing complex
vessel operations. Brochure.
Nguyen, T. D., Srensen, A. J., Tong Quek, S., 2007. Design of hybrid controller for
dynamic positioning from calm to extreme sea conditions. Automatica 43 (5), 768785.
Palm, W. J., 2009. System dynamics, 2nd Edition. McGraw-Hill Higher Education.
Parker, M. A., 08 2013. R/v gunnerus summary report offset survey and gyro-calibration
additional survey of em3002s survey of dps232, seapath, mru5 and mru5+. Tech. rep.,
NTNU.
Prieur, C., Tarbouriech, S., Zaccarian, L., et al., 2012. Hybrid high-gain observers without
peaking for planar nonlinear systems. Proceedings of CDC 2012.
Rangaiah, G. P., Krishnaswamy, P. R., 1994. Estimating second-order plus dead time
model parameters. Industrial & engineering chemistry research 33 (7), 18671871.
Rolls-Royce Marine, 2009. Icon DP - A positioning product from Rolls-Royce. Brochure.
Schultz, T., 2013. personal communication.
Seatex, K., 03 2014. Mru 5+. Tech. rep., KM.
Seborg, E. D., Edgar, F. E., Mellichamp, D. A., 2004. Process Dynamics and Control,
2nd Edition. John Wiley & Sons.
Skjetne, R., Fossen, T. I., 2004. On integral control in backstepping: Analysis of different
techniques. In: American Control Conference, 2004. Proceedings of the 2004. Vol. 2.
IEEE, pp. 18991904.
Slotine, J.-J. E., Li, W., 1991. Applied Nonlinear Control. Prentice hall, Englewood Cliffs.
Smith, C. A., Corripio, A. B., 2006. Principles and practice of automatic process control,
3rd Edition. Springer, London.
Smith, C. L., 1972. Digital computer process control. Intext Educational Publishers Scranton.
Smogeli, . N., 2006. Control of marine propellers. Ph.D. thesis, Dept. of Marine Technology, Norwegian Univ. of Sci. and Tech.
Srensen, A. J., 2012. Marine Control Systems - Lecture notes. Department of Marine
Technology, Norwegian Univ. of Sci. and Tech.
Strand, J. P., Srensen, A. J., Ronss, M., Fossen, T. I., 2001. Position control systems
for offshore vessels. In: The Ocean Engineering Handbook. Ed. Ferial El-Hawary Boca
Raton CRC Press LLC. 2001.
Tang, Y., Wu, Y., Wu, M., Wu, W., Hu, X., Shen, L., 2009. Ins/gps integration: global
observability analysis. Vehicular Technology, IEEE Transactions on 58 (3), 11291142.
96
Appendix
By differentiating Eq. (5.35), using the dynamics from Eq. (5.32), (5.33), and (5.36),
and inserting for the virtual control (, t) given by (5.37), and the controller given by
(5.38), the resulting closed loop dynamics for x becomes
>
x 1 = R()
( d (t)) + R()> ( d (t))
= rSR()> ( d (t)) + R()> d (t))

= rSx1 + x2 + (, t) R()> d (t))
= rSx1 Kp x1 + x2 ,
M x 2 = M M
= D + R()> b(t) + M
= x1 Kd x2 R()> (
b b)
= x1 Kd x2 R()> x3 ,
such that the closed loop dynamics is given as
x 1 = rSx1 Kp x1 + x2
M x 2 = x1 Kd x2 + R()> x3
x 3 = Ki R()x2
97

Topic Dynamic Positoning

Uploaded by

Copyright:

Available Formats

Topic Dynamic Positoning

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Topic Dynamic Positoning

Uploaded by

Copyright:

Available Formats

Topics in dynamic positioning

System identification, GNSS and MRU lever

Svenn Are Tutturen

Norwegian University of Science and Technology

PROJECT DESCRIPTION SHEET

Svenn Are Tutturen

Marine control engineering

Thesis title (Norwegian):

Tema innen dynamisk posisjonering: Systemidentifikasjon, estimering

Thesis title (English):

Topics in dynamic positioning: System identification, GNSS and MRU

Faculty of Engineering Science and Technology

Faculty of Engineering Science and Technology

As specified by the administration.

Lever arm estimation for GNSS and MRU sensors

Hybrid integral action for DP

Sway step tests, main propellers, port - time constant

Thruster retrofit, illustration . . . . . . . . . . . . . . . . . . . . . . . . .

Thruster configuration, Gunnerus - before retrofit . . . . . . . . . . . . .

Agility plot, thrust time constant, no yawing moment - half thrust .

Plot of p, q, and r for GPS case study . . . . . . . . . .

Lever arm coordinates, l1 , observer, experimental GPS data . . . . .

Position and integral state error first order system, no sampling .

the study of ship motion in the absence of wave excitation (Fossen,

First order plus dead time

Global positioning system

Inertial navigation system

Linear time varying

Motion reference unit

North, east, down

Proportional, integral, derivative

revolutions per minute

revolutions per second

Uniformly globally asymptotically stable

orientation vector in NED-coordinates

body fixed velocity vector

body fixed angular velocity vector

vector of Euler angles

Dynamic positioning (DP) of marine vessels are a well-established commercial application

R/V Gunnerus sea trials

1.2 Main contributions

High fidelity model

2.1 Vessel models

propulsion forces, : the forces generated by the thrusters of the vessel

Uncoupled surge dynamics

Also, for the sway motion a similar relationship is found as

First order approximation

and from (2.12) the equivalent time constant Tv is defines as

Low speed models

and the rotation matrix reduces to

2.2 Thruster dynamics

M + D = R ()b + + wind + wave1

where Va is the inflow velocity on the propeller.

2.2 Thruster dynamics

u2 : starboard main propeller

u3 : tunnel thruster, bow,

Be defining the rudder angles as

System identification by step tests

Chapter 3. Model Identification