The identification of linear ship steering dynamics using maximum likelihood parameter estimation

FRAN

STATENS SKEPPSPROVNINGSANSTALT

(PUBLICATIONS OF THE SWEDISH STATE ShIPBUILDING_{EXPERIMENTAL TANK)}

Nr75

_{GUTEBORG 1975}

THE IDENTIFICATION OF

LINEAR

SHIP STEERING DYNAMICS USING

MAXIMUM LIKELIHOOD PARAMETER

ESTIMATION

BY

KARL JOHAN AsTROM AND CLAES G KLLSTROM

Lund Institute of Technology

AND

NILS H NORRBIN AND LENNART _BYSTROM

ISBN-9 1-38-03032-2 ISSN 03'73-4714

PREFACE

Since long the establishment and the analysis of input-output relations has been a way to define the function of dynamic systems, as an alternative to that of the synthesis of ele-ment characteristics obtained by direct measureele-ments. In ship and ship model testing routines definite manoeuvres, such as a zig-zag test, furnish such _{helm-to-heading} rela-tions, but the obstacles of force non-linearities, of outer disturbances and of measure-ment noise have seriously hampered an analysis beyond the first-order approximations to the characteristics of the dynamics.

In recent years and in other fields of application "system _{identification" has emerged}

as a tool not only for the evaluation of the coefficients of the equations of

_motion

-i e the parameter estimation _{- but as a technique, which includes the design of an} experiment, the selection of the structure of the proper model to fit this experiment, and the numerical validation of the results. The Department of Automatic Control at Lund Institute of Technology has played an important role in the advancement of the proba-bilistic identification theory.

The present volume in the series of Publications from the_{Swedish State Shipbuilding} Experimental Tank forms the first and welcome document on the results of a

co-opera-tion between SSPA and the Department of Automatic Control in Lund in a field of growing importance to ship model basins, ship designers, and ship operators.

SYNOPSIS

The project design for steering qualities of a new ship as well as the prediction and

control of its behaviour when at sea may all be based on records of the motions ofscale

models or previous ships, or even on the observation of the particular ship in the

immediate past. In recent years new techniques of stochastic system identification have appeared as alternatives to the deterministic methods hitherto applied to the analysis of free-sailing experiments.

The present report covers the initial results of co-operation on ship steering problems

between the Swedish State Shipbuilding Experimental Tank (SSPA) and the Lund

Institute of Technology, Department of Automatic Control (LTH), where the theory of

the maximum likelihood parameter estimation has earlier been developed for a variety of

dynamic processes.

The role of system identification within ship steering analysis is discussed, the

kinematics of ships are high-lighted and the mathematical models for the ship dynamics are reformulated. The system identification complex is introduced with special emphasis

on parameter identifiability and estimation, and on the requirements to be placed on

input-output models.

A general program package developed for the computation of the maximum likelihood

estimate of the parameters of either a continuous or a discrete model is described. _A

special subroutine containing three different types of ship models is included.

A program for calculating the sway velocity from heading and Decca coordinates is put

forth to use the special information from manoeuvring tests.

Finally several examples of estimations of the hydrodynamic coefficients and the components of the second order transfer function, relating the yaw rate to the rudder

angle, are given.

SHIP KINEMATICS AND SHIP DYNAMICS

This chapter will discuss the steering and manoeuvring characteristics of ships and the way "models" are used as an aid to understanding.

2.1 Mathematical Models, Scale Models, and Ships

Ship model basins routinely test "scale models" prior to the building of a ship of new design, mostly for propulsive performance only; less frequently these models are used to

support the prediction of manoeuvring qualities. Whenever available the results from

scale model manoeuvring tests should therefore be analyzed _{not only with a view to}

predict the behaviour of the particular ship prototype but to add to the collected and

much-needed knowledge of manoeuvring characteristics related to hull geometry.

Shipyards run a number of delivery or acceptance trials with every new ship; at least

one ship of each class will be subject to some manoeuvring trials. There_{are several}

Manoeuvra-bility Committee of the International Towing Tank Conference presented a new ITTC code, which will pay proper attention to the needs of ship officers who will handle the ship as well as to those of the analysts, who willcorrelate the results with pre-trial

predic-tions and, again, extract as much as possibleof general information from the data.

Ship trials are too expensive and time is often too short for adequate information to be

obtained. Steering tests that can be performed in the course of a routine passage are

strongly desired. Such tests may also be used to acquire supplementary information on

the characteristics in shallow water areas, etc.

The prediction, analysis and correlation of manoeuvring performance, and the transformation and storage of information on manoeuvring characteristics are possible only by reference to the equations of motion for( or to the "mathematical model" of) the ship-and-screw-and-rudder system. Together with telemotor system and steering engine

this first-mentioned system forms the forward leg of the steered ship loop.

The ultimate information asked for by the naval architect or control engineer is the

numeric value of each of the individual coefficientsappearing in the mathematical model

chosen. The method to find these coefficients from the analysis of experiments is

"parameter identification". As will be seen in the next chapter the probabilistic identi-fication theory has the potential of derivingthe model structure as well, which motivates

the term "system identification".

As far as possible the same model structure should be used in experiment analysis and

prediction synthesis.

This structure must be soundly based on the

dynamics and

hydrodynamics of the ship. Depending on the type of experiment available and on the completeness and accuracy of the data there will stillbe a need for mathematical models of varying complexity, with and without non-linearelements. In first-order or approximate

models the coefficients to be identified are generally of a composite form, built up from combinations of the original coefficients in the "complete" equations of motion. (See below.) Such models have been widely used in the deterministic identification of ship characteristics by "equation error" type methods Nomoto (1960), Norrbin (1965)

-and by the phase-plane analysis technique - cf Bech -and Smitt (1969).

The diagram in Fig 2.1

has been compiled to indicate the role of parameter

identification in ship manoeuvring applications. A special note on scale-model-to-ship

predictions is pertinent to the subject.

2.2 A Note on Scale Model Experiments and Data Extrapolation

Scale model experiments may be divided into two principally different groups, one

including force measurements on constrained or "captive" models, the other including tests with free-sailing self-propelled models that are manned or remotely controlled, or

auto-piloted.

The principle of the Planar Motion Mechanism (PMM) was introduced by Gertler and Goodman (cf Goodman (1960)) for measurements on submerged models, which were oscillated in well-defined motions in the vertical plane. It has since then been widely

adopted for tests with surface models; cf Strom-Tejsen &Chislett (1966), etc.

FScaLe

modeL tests 8. analysis a

Free -Sailing model tests 4. Ship trial & analysis V Ship ma dcl

corr clot ion

+

Ship service operation

Fig 2.1 Block diagram illustrating tile role of steering parameter identification in scale-model

and ship

experinient analysis, and in ship operation.

1

_.l_--__---J

I I Pre-/7 Ship

,/

Ship Parameter Ship service

design know Ledge

/'/ design -,

7///// //

trials idnt. experiment

r

a a L.. a i r----i

r---r----i

-,

I Captive p Scot. i Trial Ship i p Ship model affect r manoeuvre model model tests orr.ction LPnedctb0J 1orrslationj

L.°.!!J

I I Parameter Scale Trial Ship Ship dent. if fact manoeuvre model model correction prediction correlation c or r S lot ion

steering problems. The scale model (moving in obedience to the Froude law) may be said to be a programmed analogue computer. furnishing the model structure of the dynamics - unknown to the experimenter as well as the computer facility. In common with the ship trial as discussed below this problem-orientated free-model technique has two main

drawbacks:

The solution of a new manoeuvre problem requires a complete new experiment.

New results for a modified configuration or load condition cannot be uniquely

inter-preted as to their causes.

In addition the scale model may suffer from scaling-lawdivergences (scale effects). mainly due to screw and rudder loadings being affected by a viscous boundary layer that is "too thick" and a viscous resistance that is "toolarge". During the tests an auxiliary air

screw is sometimes used to supply a friction force allowance, but this artifice will

introduce other discrepancies. Also, the dynamics of the propulsion engine will have

characteristics which are different from those of the ship.

There is no proper way to enter corrections for scale effects into the free-model

computer program", but only to apply empirical correction factors to the over-all or

final results.

A typical correction of this kind may be, say, an empirically motivated 7.5 per cent deduction of model steady turning diameter - cf Burcher (1972); much larger corrections

could well be necessary in the transient stage, however.

ln contrast to the freesailing technique the alternative captive scale model technique -with subsequent computer predictions of full scale behaviour - does accept the proper corrections for the scale effects, applied to individual "coefficients" in the equations of motion. These corrections are based on the continuousmodel/ship correlation as well as on detailed studies of model flow phenomena, which latter add to a special "block" of

knowledge.

Obviously the access to a suitable method of parameter identification will greatly

enhance the premises of model/ship correlation, and so the free-sailing model technique will gain a new potential as an alternative to the captive model testing. (Again, more or

less the same method of identification may eventually be applied to improve the extraction of data from PMM force measurements

2.3 A Look at the Kinematics

Before discussing different mathematical models it is worth while to look for guidance from the kinematics of typical ships, manoeuvring in a horizontal plane as shown in Fig 2.2. One of these ships is a dynamically unstable full-form tanker, the other a

dynami-cally stable narrow-beam destroyer.

A ship is said to be dynamically stable on straight course or in a turn of constant

curvature if, upon a small disturbance from its steady motion, it soon resumes that same motion along a slightly shifted path, without any correcting rudder being applied.

however small the disturbance met, and it will end up moving in a circle of certain radius,

in which it will now be stable with zero rudder.

The "classical" acceptance trial manoeuvre to be performed with a new ship is the "hard-over turning circle", which is in effect a kind of step response test applied to an over-damped higher-order non-linear element. From the nautical point of view the result

is expressed in terms of minimum "reach", "advance" and "tactical diameter" at

maximum helm , which may be compared with handbook values. The old figures still

hold true, the tactical diameter being typically 3.5 ship lengths for the tanker and 7 ship

lengths for the destroyer at high speed, as illustrated in Fig 2.3. If, thus, modern

technology has not changed these handbook values very much it mainly reflects the fact that ship operators have not specified any new requirements, but also that the turning characteristics are largely inherent in the main proportions of the type of ship considered.

The differences between the two ships become even more obvious if transient be compared to transient and steady turning to steady turning. The time histories of the

turning manoeuvres are shown in Fig 2.4, also including a so called "pull-out", in which the rudder is simply returned to midship.

Based on ship lengths travelled the "reach" of the tanker is roughly twice as large as that of the destroyer, indicating a larger effective "time constant". For further reference this time constant of the equivalent first-order system will be denoted by T.

Whereas the final turning circles in Fig 2.3 are about the same, the non-dimensional constant turning path curvature L/Rc (or yaw rate ,t) is also much larger for the tanker than for the destroyer; this indicates a higher effective gain. K = which may be attributed to a smaller effective yaw-damping moment experienced by the tanker hull. In Fig 2.5 the non-dimensional yaw rate obtained from the turning manoeuvre is plotted in a diagram of steady-state characteristics 1(c). The curves shown are the loci derived from turning circles at several helm angles, and in particular from the special trial manoeuvre known as the "reversed spiral". (Smitt (1967).) The full drawn branches of the spiral

loci represent stable steady-state conditions, whereas the dotted centre branch is the unstable equilibrium condition, which is characteristic for the dynamically unstable ship.

Surely, the handling of a ship is not just turning kinematics, but it involves steering on a straight course as well as repeated changes of heading in fairways or congested areas, etc, etc; at best, it is based on an understanding of the dynamics of the ship.

A good helmsman or auto-pilot will steer a large unstable tanker on a straight course almost as well as if it were stable. The diagram in Fig 2.6 shows a time history of the ,L'() locus during a simulator test, plotted on top of the spiral characteristics. In fact it will be possible to control the ship with small helm, well inside the width of the loop, if it is observed that the yaw rate should never be allowed to develop too much before checking rudder be applied; for checking the rudder must be moved over to the opposite side of the dotted curve.

Special trial manoeuvres (such as the zig-zag test by Kraemer (1934)) have been devised

to describe the controllability and response of a ship in a relative manner, as compared with other ships or as predicted from ship model basin experiments with a free-sailing scale model. The standard or Kempf zig-zag test procedures - using 100 opposite helm

Advance Tactical diameter

-Tanker Helm execute Tanker _Destroyer

Fig 2.3 Typical turning circle tracks for a 300 m tanker and a 110 m destroyer.

1.0 _1.0

Or-35

0

300 600 900 t S _{0 300}

Fig 2.4 Schematic time histories of turning circle manoeuvres with pull-outs.

0

\

Dc st raycr t S 600 m 0 I,

L

Fig 2.5 Non-dimensional steady-state steering characteristics for the unstable tanker and the stable destroyer. The diagrams also illustrate the definition of K' (and of "effective" K) as well as one

1. 3 - 0.2 -0.3 6 50 5 a a

Fig 2.6 Time history of the straight course helm angle - yaw-rate locus ofa 200 000 tdw tanker during manual control in the SSPA simulator, plotted on top of the steady-state steering characteristics.

00 50

execute at 100 change from initial heading, or 200 helm at 20° change of heading - both result in oscillations of almost the same "period", which for different ship types may vary

within a range from 6 to 14 ship lengths. The "dynamic gain", i e the ratio of double

amplitude of yaw rate to that of helm angle, should reflect the inertia in the motion, but

the results are obscured by the influence of non-linear damping. The multitude of explicit

data collected along these lines has again failed to add very much of useful information on the dynamics of different types ofships.

2.4 The Mathematical Models

In applications to ship motion studies in general the ship may be regarded as a rigid body

with six degrees of freedom, i e three translations and three rotations. The following

relates to the analysis of the motion proper in response to the control actuations achieved at the rudder and screw.

The ship moves in the interface between two fluids and the generation of waves and

vortices are real world facts, since long well-known to be the possible sources of "memory

effects" in the dynamics. A functional approach to the handling of such effects within a linear theory has been presented by Bishop, Burcher and Price (1973). In the application of ordinary differential equations the equivalent approach would require coefficients which vary with the reduced frequency. w' w L/V. A well-documented experience of

calm-water ship manoeuvring analysis has shown that this frequency-dependence may usually be ignored in favour of a quasi-steady low-frequency theory, which more readily

accepts the unavoidable complications due to non-linear hydrodynamics.

The formal equations of motion are conveniently written down in Euler form, with reference to body axes in a right-handed system where x is positive forward and y positive

to starboard. The same convention infers that a stern rudder for lateral manoeuvringwill

produce a positive (starboard) yawing moment when switched in negative (starboard) direction. (Cf Fig 2.2.)

For the purpose of this study the problem will be greatly simplified in that the

dis-cussion is confined to motions in the horizontal plane only, and in that the influence

of heeling in a turn is ignored. Mostly this heeling is small in itself, and in other casesthe

coupling effect on sway and yaw has been shown to remain small. The turning of aship is

initiated by rudder control but the subsequent motion is

largely governed by the hydrodynamic forces due to drift (or side-slip) and yaw. This is in clear contrast tothe case of an aircraft, where side-slip is to be kept as small as possible, whereas the

horizontal turn is achieved by proper banking. In comparison with the aircraft steering problem. again, the ship problem is complicated by the dominating influence of added

inertias and large-value damping non-linearities.

Referring to Fig 2.2 the steering manoeuvre will be described by three second order differential equations for the linear translations u and v and the yaw rater Li. (The general expression for yaw rate includes coupling with heel and pitch, of course: r = 1i cos

0 cos

-

sin p.) The heading angle ,Li may indicate the attitude in relation to a datum

100

_f

200 t (sec)

!Ln.

contr.

Fig 2.8 Total and linear contribution to swayand yaw accelerations in the mathematical modelling of

4 10 (rod/s2) 0.5

&

' : Lin. contr.200

:-

I-0.5 (sec)

Fig 2.7 Total and linear contribution to sway and yaw accelerations in the mathematical modelling of a 100/100 zig-zag manoeuvre with a twin-screw container ship.

x0(t) = f Vcos(Li-j3)dt

0

y0(t) = fVsin(i-I3)dt

0

The derivation and application of the Euler equations to the ship manoeuvring problem may be found in papers by Davidson & Schiff (1946), Norrbin (1960) and Abkowitz (1964). More or less "complete" non-linear models have been put forth by Abkowitz (1964), by Eda & Crane (1965) and by Norrbin(1970).

As the parameter identification procedure presented in the following chapters will treat the linear problem only, the equations of motion pertinent to the present task are written in the form of eq (2.3) below. The mass inertias proper, displayed in the left side mem-bers, have been complemented by those parts of the hydrodynamic forces, which are usually named "added masses", and which will mean - as already pointed out - a

sub-stantial increase of virtual inertias.

The first position in each of the right side members is reserved for control forces,

which at constant speed and screw loading (or constant speed u and revs n) may be taken as proportional to helm angle ö. (If the rudder is not within the influence of the screw

race the rudder force is of course not a function of screw loading.)

The right side members in second and third positions are approximations to damping

forces in steady sway or yaw, and they include pure mass forces as well as the con-tributions due to hydrodynamics. If the forward speed u is fairly constant the four

a certain manoeuvre. From a knowledge of the initial position ship track and heading are obtained by integration over time t, or over the distance t' = I (V/L)dt in ship lengths

sailed. Thus

(2.2)

x0(t) = f (u coslli -v sinili) dt

0

y0(t) = f (u sinai + v cosil') dt (2.1)

0

1i(t) = fçlidt

The drift velocity v is related to the side-slip angle = - arc sin - and therefore the

terms may obviously be considered as linear. If speed variation is handled separately their

linearization is still valid in terms of side-slip and path curvature. Their coefficients will also appear in the analytical condition for dynamic stability on a straight course.

The mass forces in position (4) should be included in the formal linear model although

they are mostly small enough to be ignored.

All terms in the positions (5), (6) and (7) are non-linear. The non-linear "residuals" indicated in the last column positions are complicated functions of the hydrodynamic

state of motion, and they will be required in the proper description or simulationof

"tight" or radical manoeuvres. As an example, the lateral force due to side-slip no longer

increases in proportion to the drift velocity v but much faster, due to the additional contribution of cross-flow viscous drag.

(2.3)

It should be especially observed that the force balance along the x-axis is completely non-linear. The added "resistance" in a turn is dominated by the fore-and-aft component

of a "centripetal force", in which the mass force is strengthened by a force of

hydrodynamic nature. (The centripetal acceleration is usually written as V2 /R, where the

momentaneous radius is given by R = V/,1i. The acceleration thus has x and y components equal to ui,Ii and vL' respectively.)

The diagrams in Figs 2.7 and 2.8 illustrate typical records of helm manoeuvre and speed and heading response (top) as well as "linear" and total contributions to the sway and yaw accelerations in simulated 100/100 and 200/200 zig-zag tests with a container

ship. Whereas the linear analysis may possibly have some justification when applied to the

100/100 test it is quite clear that this is no longer true in the case of the 200/200 test.

If, thus, even the 100/100 zig-zag test is far from a linear manoeuvre it would be

tempting to revert to an analysis of the routine coursekeeping record. Unfortunately, now, such records often display the very-low-frequency limit cycle type steering

asso-ciated with non-linear elements as shown in Fig. 2.6; the non-linear element may be

represented by the ship hydrodynamics, the steering engine, or the helmsman's mode of control. A forced bang-bang type small-helm control executed at random intervals will

produce a nearly straight course and a response record covering a wide frequency

spectrum, however, and it will be seen to furnish an adequate input for linear analysis. (The response frequencies now considered may still be said to be "low" with respect to

the hydrodynamic frequency effects mentioned earlier.) Spec

Linear Linear Non-Linear

inertia

terms

contro'

terms damping terms

nertio

Cou pLungs inertia compLings

residuaL

terms

Position 0 2 3 1. 5 6 7

X (m-X,,)u .bx..1X.,)4+(rn+X,)v,i, +X

y (m-V) Y(u,n)6 V55uv +(Y-m)u4, +(Vr_mx)4 +V

2.5 The Linear Model

From the above it follows that the linearized model will include perturbationsin sway

and yaw only, i e the Y- and N-equations. That does not _{mean that speed loss must}

necessarily be ignored, but it can be handled separately. In general the forces involvedare

all proportional to speed squared or roughly, in normal steeringmanoeuvres, to forward speed squared, u2.

It _{is convenient to introduce non-dimensional coefficients to characterize the ship}

geometry in the dynamics. The "bis" system - Norrbin (1970) - is particularly handy in more general cases; here all forces are related to the displacement force of the ship, linear

accelerations are related to the acceleration of gravity, and the speed is_{given by the}

Froude number. In dimension-true equations the non-dimensional coefficients (or

"derivatives") appear together with multiplication factors in powers of ship length L.

On the next page the Y- and N-equations (2.4) are presented in this way. The formal forms appear at top and bottom, and the equations may be reduced as indicated in the way of the arrows. The example shown refers to a single screw tanker, where the rudder force is taken to be proportional to rudder angle 6 and to n2, i e to RPM squared. The

asymmetric lateral force, which is an effect of the single-screw action, is likewise_{taken to} be proportional to n2. When speed and RPM variations_{are moderate the rudder force will}

vary with 6 only, and the asymmetric force will be balanced by a small residual helm In the absence of outer disturbances 6

= -

10 _{may be a typical mean helm value on a}

straight course; the ship with a righthanded screw has to carry a small starboard helm.

2.6 Outer Disturbances in the Linear Model

A record of heading angles and helm positions in _{a ship experiment will exhibit the}

influence of water currents and wind loadings. Referring to Fig 2.9 it will be observed that a steady current will not be sensed by the ship hull; it must of course be corrected for, by heading upstream, when aiming for a certain position. The effects of the steady

wind, on the other hand, must be balanced by water forces through_{a combination of}

side-slip and correcting helm.

A typical diagram to show the variation of lateral force and yawing moment due to a "relative" wind of given angle YR off the bow of a modern tanker is reproduced in Fig

2.10, from van Berlekom, Tragârdh and Dellhag(l974).

During manoeuvres in a true wind of constant velocity and direction_{the magnitude as} well as the direction of the relative wind will change with the heading of the ship, and the force and moment will be complicated functions of this heading;_{especially it may be seen}

that the lateral centre of wind pressure may move rapidly as the ship turns through

certain angles. It is important to remember that the ship isno weathercock, however, and

that its tendency to luff into or bear away from the wind

_{always depends on the} combination of wind and water forces.

Within the application to the linear analysis the wind force and _{moment may be}

assumed to vary proportionally to the change of heading, so that yw_{= + YL, and}

m2 2 + I. N0/ (i_NrC _{- 0.0000074} - 0.0000037 - 0.0000037 /(J-N.) 2rw' 1N /raL22 nfl -0.40000037 3.7 10 1 1N2 an m m252 L1

f'c6

c26 329.21 0.210 c26 0.000638 c26 c26 m 1

V-i

U)i) 1 -0.525 U)j) -0.525 u4 urm) /m u4 /(rr-Y.) 2cc5 V c2/n2 a2 6 0.000382 34 2 § 0.0130 2.0 6 0.0260 6 Y6/(rn-V) 6 mS2 CV-m)/Cm-V.) U -0.314 ) - 0.314 -2.74 4) (Vr_mu)/(m_V4 ) 4) ms S -1S -1c )Nr mU)/)J- N) 4) -0.0680 -0.0077 8 8.74 4) -5.00778 U (N-mn)/)i-N.) u 4) (N,,mx&mL2 u4) -0.000778 u 4) 329.2' -0.256 U),) L Na6 Ut) a mc (N-Q/CJ-NC + + N5(i-N - 0.0012 2 -0.000145 - 0.00122 V -0.0000166 0.74 + (I) - 0.00122 - 0.0000166 U V (N.-x6)/(JNt) (J-N) U V CJ-NjJ/m L2 4) CN-x6)/mL2 P{/m L2 US 0.100 4; -0.000122 -0.00000166 U, + 90576+00424 329.21 -0040 329.22 -0.180 "V 4) 1_i N X6 L2 uV Un + m m 2 1-v-S6 Yuv US 1.67 329.2 - 0.050 329.21 -1.21 US + 1.67 -16.'.6 -0.00368 uV Cm -) I m V (y. -C/rn Yuv/m US N6 / Ci - N) 6 -0.000629 6 + - 0.000314 2.0 § -000000923 34 n2 6 N.) 6 1-N /m12 2 ccó -0.000000923 2 C + 3292 2- 0.100 6 tNCC6 2O m m 2 mc rfl2 m2 2 m 2 nfl 329.2 + 0.000 263 / m -V/ 0.000158 0.00015 8 + 0.000316 Y 1Cm-V.)0 m12 (Y.-)/(m-Y.) V /Cm-V. US V U V -0.00220 tl - 9.86 4) + - 0.00220 + - 9.86 0.0192 - 9.86 Y /(m-V C v- QI Cm - Y) V ç2 m m

vc

The homogeneous cross-current case: The ship is drifting with the mass of water with no hydrodyna-mic side-slip ( 0) or correcting helm.

The beam-wind case: Helm angle ö and side-slip _{are both adjusted to make the ship produced in force} and moment balance. The force Y' and the corresponding yaw moment vary with magnitude V and

direction _{of rclativc wind as examplified in Fig 2.10.}

Fig 2.9 _{'Drift" due to cross-current and due to beam-wind in}

cIO -5 I0 30 50 90 120 150 1.0 1 1)

..u.uRuu..Mauu.

-20 -30 30 50 90

....

uIi

j

a... a

c ¶0' ISO y()

The non-dimensional coefficients are based on product of relative wind stagnation pressure and refe-rence area L2. Symbols refer to specified deckhouse configurations. (From van Berlekom, Tragârdh

and Dellhag (1974))

13 3' '1 Wnd speed 5 3 W,nd speed

(By courtesy of the Meteorological Station of the 2nd Helicopter Division in Goteborg and of the Sw Meteorological and Hydrological Institute in Stockholni)

Fig 2.11 Examples of wind speed fluctuations during 2-hour intervals, from anemometer records 10 m above ground at Save Airport.

I

10 20 0 20

i

finite negative values. In a (relative) beam wind

0, but here the three other

coefficients have finite values. The constants corresponding to the initial average wind

load will not be distinguished from the residual rudder load; thus, say, F =

-and F2 =N'N660.

In the real world the true wind fluctuates in "strength" and direction in a way which is known to be influenced by the special meteorological conditions prevailing. Trial codes usually stipulate that ship performance tests are not to be run in wind conditions above 3 Bft, in which case the fluctuations just mentioned may be described as the effect of a homogeneous isotropic turbulence. The time histories included in Fig 2.11 will then be fairly typical of the variation of wind velocities close above a flat surface such as a wide field or the open sea. In view of the long time constants seen to characterize the lateral response of large ships it is often reasonable to assume that the variations of wind force

and moment are mainly realizations of white noise.

2.7 The Linear Model Continued

If v, r ( i) and Ill are considered as three state variables, and if the expressions of the last

section are added to eq (2.4), then the linear equations may be written in matrix form as

m-Y;

mxG-Y 0 2

mxG-N,, mk-N 0

0 0 1

YflTh

N NrmuxG N'

0 1 0 v r + Y& N, 0

6+

0 F1 F 2 (2.5)

In the next chapter this equation will appear (eq (3.1)) with the coefficients in suitable matrix notation, and formal expressions will be derived for the transfer functions, eq (3.3) and (3.5).

Whereas the control function 6(t) may itself include a heading reference the open

ship-and-screw-and-rudder loop has no inherent preference for any one heading. The characteristic equation s (2 + a1 s + a2) = 0 thus has one root equal to zero, and it may be

shown that the two other roots normally are real. By convention the root to the right on

the real axis is denoted s

= - T

. The appropriate transfer function from helm to

heading may read

K(1+T s)

G1(s) = 3 (2.6)

'

s(TTs2 +(T +T)s+ 1)

with accepted symbols following Nomoto (1957). The same time constants T1 and T2 will appear in the transfer function from helm to sway velocity:

K(1+T s)

G2(s) = Gv6 v 3v (2.7)

T1Ts2 +(T +T)s+ 1

v r ill

The two roots are negative for the stable ship, whereas s1 turns positive in theunstable

case. For the stable ship, again, K is negative according to the sign convention

_{- cf}

Fig 2.5 - and changes its sign with T1 _{. (In the Bode diagram of the transfer function}

in its frequency form the change of sign merely means an additional 1800 phase lag at

the low-frequency end.)

The formal relations between the constants in the transfer functions above and the coefficients of the original equations are listed in the table on thenext page.

In principle the tranfer function (2.6) may be identified from a set of frequency

response _{tests. Nomoto (1957) first showed that the low-frequency or first-order}

approximation (valid for w ( T ,or slightly higher)

K

G4(s) = _(2.8)

s(1+ Ts)

contributed a powerful tool to the analysis of many steering problems for stable ships. According to this approximation the effective time constant is T = T1 + T2 - T3, and

Nomoto also pointed Out that by this definition the first-order steering equation

T+=Kb

_(2.9)

predicts the correct final response of a linear system to a rudder step input. Due to the non-linearities present, however, the normal turning circle test will furnisha kind of

effective gain and time constants K and T5. (Cf Section 2.3.) In order _{to obtain true}

linear coefficients from the analysis the amplitudes of theyaw and sway motions must be

kept small. The prediction of moderate to hard manoeuvres will then require additional information. This latter is even more obvious for the unstableship.

In terms of "stability derivatives" the analytical criterion for dynamic stability on a straight course may be written as an unequality of two leverarms,

mxG Nur N

l1 =

m - _Yuv

>0

_(2.10)

provided m

-

_ur> 0, n.b. If _{ur > m, then this hydrodynamic force due to yaw must}

apply aft of the centre of gravity G in order to _{secure dynamic stability. The physical}

interpretation of this criterion is that in the motion initiated by_{a disturbance the yaw}

damping moment and the sway damping force will tend to stop the yawing. The linea-rization may in a similar way be applied to the ship in a turn of given radius and with a side-slip * 0, and it will then be seen that the stability in the _{turn is increased over that} on a straight course. It was seen earlier in this chapter that the unstable ship _{was in fact} inherently stable in a certain turn to port or starboard.

TABLE OF TRANSFER FUNCTION

CHARACTERISTICS

in

true

dimensions

in non-dim "prime and

bis" systems K Y6u K' Fn1 Y6' L6( -21uU6

trLl, cf. e2(1OJ

L I, -rN(6) /1(ç) m - 1ur Lr rn Y. ( -L 1--1 Dirnl 'UUOU K" K FL iuu*

L-('

Iux r L v v -r v uv r v Dim LI1 1 I3 m-Y mx- N ' 3 nI 3 m'-Y4' mx- N V u' - N' m-Y L- m,_eY

1Y'

u -L -t -Dim I T3 mx mk1- Nr l3uFL 13 m m'k'- N ' - V k'2- N' mx -m':, - Y -rn -1ur -1 - Yr -DimI v

L3v

12 rnx6- 1r m- V mk- Nr I'I'uF2-iT" !uI II m'x-Y mx-N rn-c m'k- N x-Y' x'-N1 1-Y k;' V mlur ( -1un mYur t,) u 2 Y m'-Y Y, m'-Y

i;- t;

Y 1-Y Y

'1r

Dim - t 12 m-Y mx6-V L mk-Nr mu6- N T1+I2 F,(I'+T') Y!-m' m'x-Y L' m'h-N m'x-N LL.' x'-YP k"2-N' r m- 1ur m- 1ur 1uv r m'-Y ni-V r 1 " 1 1ur Dim I ( 1r t) U ( I+ 12) L' L'

2.8 Approximations with Added Non-Lineanties

It was seen above that the validity of the linear relations is limited to moderate

amplitudes. In particular 'c = K6c fails to describe the steady state characteristics as

derived from spiral tests with marginally stable or unstable ships, and diagrams as in Fig

2.5 suggest that the linear relation be replaced by an expression with a cubic or abs-square term, such as

cKc'c = K6

(2.11)

As a first approximation to the non-linear problem these equilibrium conditions may

be introduced into the steering equations, which correspond to the transfer functions

(2.8) or (2.6); Norrbin (1963, 1965). In the equation

T1T2i +(T +T)Vi+

' - ,KIiIi1i=K(&+T3)

(2.12)

the presence of an explicit abs-square term may be argued for with reference to the

dominating role of viscous cross-flow resistance and to the fact that the phase difference

between yaw and sway remains fairly constant during the course of low-frequency

manoeuvres. Thus the equation may be especially useful in standard zig-zag test analysis,

and in the analysis of small deviations from a pure turning circle.

Bech & Srnitt (1969) generalized the non-linear assumption by use of an unknown function H(Li) to be derived by curve fitting to the reversed spiral results, and they also

put forth a convenient way to find the time constants from graphical analysis of the zig-zag test phase plot of ,1i (st').

In the context of the present work an equation such as (2.12) will be suggested for predictions of non-linear manoeuvres on the basis of linear characteristics, that have been

obtained through the parameter indentification procedure to be described next. It is

anticipated that similar procedures will be developed in a future for a direct evaluation of the most significant non-linear contributions to the original equations of motion, (2.3).

3. SYSTEM IDENTIFICATION

This chapter presents a brief review of system identification methods with particular

emphasis on their use in determination of ship dynamics.

3.1 System Identification Methods

System identification is a collection of techniques for obtaining a mathematical model of

a dynamical system. It includes the following steps Planning of experiments

Selection of model structure

Parameter estimation Model validation

The experiments are performed by perturbing the system inputs and recording the

system outputs. The result is a record of inputs and outputs which is the basis for the other steps in the solution of the identification problem. Experimental planning involves many factors, choice of input signals, sampling rates, and recorded outputs as well as

selection of signal generation and recording equipment. The choices of sampling rates and

the input signal wave-form are very important for the accuracy of the final result. The

selection of input signals also involves the decision if the perturbations should be

introduced as changes in the reference values of available regulators, if the regulators should be removed or if it is sufficient to use the input signals generated by the feedback itself. Many of the decisions related to experimental planning can not be done rationally without any knowledge about the oroperties of the system. Since a proper experimental planning is of paramount importance for the final result, the identification problem is often solved iteratively, by repeating the steps I through 4 many times.

A suitable model structure is usually chosen based on a priori physical knowledge. For special classes of models, e g linear systems, it is also possible to select canonical models which do not rely upon physical insight (black box models), It is a substantial problem to

decide upon a reasonable model complexity.

There are many techniques available for parameter estimation. Most of them can be

formulated as optimization problems. A criterion which defines the fit between the

model and the experimental data is postulated. The model parameters are then adjusted until the best fit is obtained. Many different criteria have been proposed like input error,

output error, prediction error etc. The result will sometimes depend strongly on the

chosen criterion as is shown by examples in the following chapters. A large number of

techniques for the optimization have also been proposed. In combination with the

possible classes of models and possible criteria this gives a rather bewildering impression on the newcomer to the field although the principles are very simple.

The purpose of the model validation is to decide if the model obtained is consistent

with the data. This is a most difficult problem. The available methods are based on common sense and statistics. The most reliable results are generally obtained if the

esti-of giving an empirical foundation for finding the effects esti-of scaling on determination esti-of

ship dynamics.

The procedure to determine ship steering dynamics is roughly the following.

Experiments are performed by perturbing the rudder. The resulting motion is observed by

measuring variables like heading angle, yaw rate and the velocity components. Under

certain circumstances the experiments can even be performed using a regular autopilot, which means that the experiments can be made during ordinarytrips with only marginal perturbations of the ship's mission. A mathematical model like the one discussed in Chapter 2 having certain unknown parameters is then fitted to the experimental data

using some parameter estimation method. The unknown parameters reflect the lacking

knowledge about the hydrodynamic forces. By applying various statistical methods it is then investigated if the assumptions made when formulating the model are consistent

with the observed data. It can also be attempted to compare models having different

complexity to see if a more complicated model gives a significantly better fit. Some further details about the procedures used are given below.

Several parameter estimation schemes have been applied to determine ship dynamics.

General aspects of partly philosophical nature are given by Schmiechen (1969). In

Kaplan, Sargent and Goodman (1972) it is proposed to use an iterative methodbased on

least squares combined with calculation of the sensitivity derivatives. The method

requires that all state variables are measured. In the same report it is also suggested to

use extended Kalman filtering when only certain components of the state are measured. The technique was applied to determine hydrofoil dynamics. The theory is based on the assumption that a continuous time record of the outputs are available. This means that

very short sampling intervals (about a second) must be used. Schmiechen (1972)

proposed an output error method. He also outlined methods for on-line updating. A

deterministic approach using the output error criterion and a gradient technique based on calculation of sensitivity derivatives is also proposed by Sandman and Kelley (1974).

They report good results on noise-free simulated data. The non-linear steering dynamics of

a small training ship was determined by van Amerongen et al (1975) using the output error criterion and a hill climbing procedure. The non-linear part, however, was estimated by a least squares regression of a polynomial to the spiral curve.

Sukselainen (1975) also used an output error method with analytic calculation of

gradients to determine a model for a radiocontrolled scale model. He experienced some

convergence difficulties but he was able to get his procedure towork using other data.

It is also proposed by Paloubis and Thaler (1972) to use a standard least squares

procedure to determine the parameters of a linear transfer function model. A continuous time approach is used and consequently it is proposed to sample 1/10 to 1/100 of the shortest system time constant when using discrete time measurements. The technique has

been applied with good success to determine parts of a ship boiler model.

The maximum likelihood method was used by Astrom and Källstrom (1972), (1973), (1976) to determine models for a freighter and a tanker. In these works a continuous time model was fitted using discrete time measurements. The theory admits long sampling intervals (10-30 seconds) which reduces the computational load. Schweppe (1973) used the maximum likelihood method to determine a model for submarine dynamics using

mation. Model validation also includes the problem of deciding between several models having different structure.

System identification played an important role in classical _{control theory. The fact} that the transfer function of a system can be determined by frequency response analysis or by transient response analysis was very important for the success of classical control theory. There has been a substantial development of system identification methods in the past 10 years. The purpose has been to obtain techniques that will give the_{state space} models, which are the starting point for applying modern control theory. This includes modelling of both the system dynamics and the disturbances. A significant _{trend has also} been to relax requirements on experiments in return for increased _{computations. This is} of course motivated by the fact that experiments are time consuming and costly but that

numerical calculations are becoming cheaper and cheaper.

There is _{a substantial literature on system identification. Good sources are the} proceedings from IFAC Symposia held in Prague 1967 and 1970 and in the Hague 1973. IEEE published a special issue entirely devoted to the identification problem in December

1974. There are also several survey papers, Aström and Eykhoff (1971), Bekey (1970) and Nieman et al (1971), as well as books, Eykhoff(1974) and Graupe (1972).

3.2 Applications to Ship Dynamics

System identification techniques have been applied to determine ship dynamicsfor a long time. Nomoto (1960) suggested a straight-forward least _{square method for identifying the}

two coefficients of eq (2.9) by using the integrated response to a succession of ramp inputs. An interesting application of frequency response methods for determination of submarine depth control dynamics was given by Garde and Persson (1960). _Submarine dynamics was determined by Booth (1975) using a pulse response method. Windett and Flower (1975) used a pseudo-random binary sequence (PRBS) input signal combined

with a correlation method to determine the roll dynamics ofa small frigate. Neither

frequency nor transient response methods have, however, been extensively applied to determine ship steering dynamics. One reason is that both methodsrequire precise inputs

which may also give large course deviations. Because of the long time _constants

measurement of the frequency response will also be time consuming.

The development of experimental techniques to determine ship dynamics has therefore followed another pattern. The spirit has been to design suitable apparatus for captive tests that permits generation of precise motions, and equipment to measure hull forces. In this way the identification problem reduces to a static parameter estimation problem. See e g Strom-Tejsen and Chislett (1966). It is of fore-most importance in the design stage of a new project, which is always the subject of a variety of scale model testing. For

prac-tical reasons the technique is limited to such scale model experiments. (Cf Chapter 2.)

The techniques for system identification developed in the field of automatic control are aimed at determining the parameters of an arbitrary dynamical system. The methods can therefore also be applied to determine models for ship dynamics. They apply both to experiments on scale models and on full size ships. Consequently they have the potential

13 I

liii

I

air1

a

r'+'b

23 21 2

OJLIJLOJ

LO

where f1 and f2 are disturbances and

ra

a

-

mxG

-i

-r

Yrmuo

Y Y F

mxG - N, mk - Nrj

[Nv Nr - muoxG N' N6 F2

[11

a a 21 (3.1) (3.2) The parameters a13 and a2 depend on the windforces; they will vanish if there are no

disturbances. The forward velocity is assumed to be constant and denoted u0.

The transfer function relating heading angle to rudder is given by

b s+b

G1(s) = 1 2

(3.3) s3 + a1s2 + a2s + a3

simulated data. Maximum likelihood estimation of a discrete time model was applied to the ship M/S Esquilino by Tiano, Piattelli and Leccisi (1973). Compare also Dagnino et al (1973). Discrete time models for rolling and steering were determined by Flower and

Towhidi (1975) using the maximum likelihood method. General aspects on the

identification of ship dynamics are given by van Leeuwen (1973) which briefly describes application to data obtained from a large simulator.

3.3 Model Structures

Different mathematical models for ship steering dynamics were given in Chapter 2. in addition to the models discussed there the state variable form of the equation (2.5) is also of interest. This model is obtained simply by solving (2.5) for the state variables v, r and

ii, i e a 12 a 22 b 13 b 23 11 21 f21₁I

r

ra

a I 11 12 dt r

[

a 21 a 22 LO 1

where the parameters are related by a

= -a

-a

1 11 22 a

= -a

a

+a a

-a

2 12 21 11 22 23

a = -a

a

+a a

3 13 21 11 23 (3.4)

b=b

1 21

b = a

b

-a b

2 21 11 11 21

The transfer function relating sway velocity to rudder is given by

c s2 +c s+c

1 2 3

G (s)

-

(3.5)

2

s3+as2+as+a

1 2 3

where a1, a2 and a3 are given by (3.4) and

c=b

1 11

c = -a

b

+a b

2 22 11 12 21 (3.6)

c

= -a

b

+a b

3 23 11 13 21

Note that if the effects of wind are neg'ected, then the transfer functions G1 and G2 are

reduced to the expressions given by equations (2.6) and (2.7) respectively.

In the experiments the signals will be sampled. The sampled model can be represented by the difference equation model

y(t)+a1y(t-1)+ ...+ay(t-n) = b1u(t-l)+ .

+X[e(t)+ce(t_1)+...+cne(t_n)]

(3.7)

where fe(t) } is a sequence of random variables which accounts for the combined effect of disturbances and measurement errors. The model (3.7) will represent the sampled version of(3.3) even if there is an arbitrary time delay in the rudder servo.

3.4 _{Parameter Identifiability}

Before it is attempted to estimate the_{parameters it will be investigated if the parameters}

of the model can at all be determined from the _{experimental data. For the particular}

models discussed here this problem is solved in Astrom and Kallström (1973, 1976). The result is that 5 _{parameter combinations can be determined from an experiment where the}

rudder and the heading angle are measured, i e the coefficients of the transfer function (3.3). If the sway velocity is also measured, 3 additional parameters corresponding to the coefficients of the numerator of G2 in (3.5) can also be obtained. As a consequence all the 8 parameters of the state variable model (3.1) can be determined if both the heading angle and the sway velocity are measured. All 14 parameters of the equations of motion,

I e m, XG, _{u0 and the 10 hydrodynamic derivatives can not be determined.}

3.5 Parameter Estimation

Having discussed the problem of identifiability, the parameter estimation problem will now be treated. Both the estimation of parameters in a state model and the estimation of

the parameters of a transfer function model

_{can be formulated as a problem of}

determining the parameters in the stochastic differential equation dx =

_{Ax dt + Bu dt + dw}

(3.8)

It

_{is assumed that the initial state is}

_{a gaussian vector with mean value m and}

covariance R0 and that {w(t), 0 t °°} is a Wiener process with incremental covariance

R1 dt which is assumed independent of the initial state. Assume that an input signal has been applied to the system and that the output has been observed at discrete times to,

t1 ...tN with a measuring device which can be characterized by

y(tk)=Cx(tk) Du(tk) + e(tk)

(3.9)

k=0,l,.. .,N

The measuxement errors te (tk)} are assumed_{to be independent and gaussian with zero}

mean and covariance R2. It is furthermore assumed that the _{measurement errors are}

independent off w(t), 0 t } and of the initial state.

The model (3.9) implies that the measuring _{instruments are such that they give an}

output signal which is the instantaneous value of a linear combination of the state

variables. The errors of measurements taken _{at different times are independent. The} equation (3.9) is a good model when the sensor dynamics is considerably faster than the system dynamics and the measurement errors _{are so small that it is not reasonable to}

average measurement signals over intervals comparable to the sampling intervals.

In the particular case of ship dynamics the model _{(3.9) is reasonable because the}

shortest time constant of interest is about 5-60 s, and the sampling interval may be

chosen smaller than that. All sensors have dynamics with time constants shorter than 1 s, and the measurement errors are about 0.10 in heading, 0.02°/s in yaw rate and 0.05 rn/s

It is thus assumed that an identifiable model (3.8), (3.9) is given and the problemis to determine the identifiable parameters from observed input-output pairs. The maximum likelihood method is used to solve the parameter estimation problem. This method is known to have several attractive properties at least for large data sets. See e g Ljung

(1974).

To obtain the likelihood function, i e the joint probability density of the observed

outputs assuming all parameters known, we introduce

y(t0) 1 Y(t1) _(3.10) y tk

y(tk_ )

y(tk)

i e a vector consisting of all outputs observed up to and including time tk.

Assume that the probability distribution of y (tk) has a density p (Ytk). It then follows from the definition of conditional probabilities that

P(Ytk) =

P(Y(tk) IYtk

l)PYtkI

(3.11)

Repeated use of this formula gives the followingformula for the likelihood function

L = p (Ytk) = p(y (tk)

I y

-

1)

p (y (tk

_-

1 IY

-2)

p(y (t1) ly (t0))p(y (t0))

The likelihood function can thus be conveniently written as a product of conditional

densities.

In the particular case of a model described by (3.8) and (3.9) all random variables are gaussian and the conditional density is also gaussian. The logarithm of the likelihood function can then be written as

N

- log L

=-k=0

log det R(tk) +

(3.13) (3.12) N 2 _k=0

T(t)R

(tk)e(tk)+const

where

c(tk) = y(tk)9(tk tk_I)

and '(tk Itkl) denotes the conditional mean of y(tk) given

Ytkl and R(tk) the

con-ditional covariance. We have

9 (t0 t-1

) = Cm + Du (t )

0 (3.15)

Also note that _c(tk), k = 0...N are the residuals or innovations of the output

process. The conditional mean '(tkftk_1) and the conditional covariance R(tk)are easily

determined recursively through the Kalman-Bucy filtering theory. Seee g Astrom (1970). We have

S(tk Itk._l) = C(tk Itk_l) + Du(tk)

(tk Itk) =

x(tkltk_l) + K(tk)e(tk)

-

(tJtk) =A(tItk)4-Bu(t)

_tk t _tk+l

K(t,K) = P(tkltk_I)CT[R2 +CP(tk

_{tk_l)C1}

(3.16)

P(ttk) = P(tk Itk_l)

_{K(tk)CP(tk tk_I)}

P(tIt) = AP(tltk) + P(tltk)AT-4-R

_tkttk+l

R(t,) = R + CP(tk Itkl)C

The computation of the likelihood function is thus easily done recursively. A description of the program LISPID, which performs the calculations, is given in the_{next chapter.}

If we are interested only in an input-output model, the parameter estimation can be

simplified significantly if the parameters of the difference equation model (3.7) are

determined directly, as was described in Astrom and Bohlin (1965). An interactive

program IDPAC, which performs this, is described in Gustavsson et al (1973). This

program minimizes the loss function

V = 2_(t) (3.17)

which is equivalent to maximizing the likelihood function. The maximum likelihood

estimate of X will be

= 2 çr (3.18)

N+l

n

where _nis the minimum value of Vn.

It should be observed that the maximum likelthood estimate can be interpreted as the estimate which minimizes the one-step prediction errors. This interpretation does not require any assumptions on the statistics of thedisturbances.

3.6 Model Validation

The purpose of the model validation is to find out if the model obtained is reasonable and at least does not violate the assumptions made when deriving the parameter estimation scheme. In the particular case this is done by analysing the autocorrelation function of the residuals and the cross correlation between the input and the residuals. To compare models of different complexity the following criterion introduced by Akaike (1972) is

also used:

AIC = - 2 log L + 2 (3.19)

where I. is the maximum of the likelihood function and i.' the number of estimated

parameters in the model. According to Akaike the quantity AIC should be minimum for the correct model structure. The maximum of the likelihood function will increase with increasing model complexity, but the term 2v in(3.19) will make AIC grow linearly with the number of parameters. Hypothesis testing analogous to the F-test used in regression

analysis is also used to decide between models having different complexity.

Analysis of the deterministic properties of the model is also an important part of the model validation. This includes investigation of pulse and step responses and calculation

of time constants, gains and frequency responses.

It is shown in the specific examples in the following chapters that the model validation procedures give valuable insight into the complexity required to explain the measured data. Indirectly this has provided an understanding of the effects of disturbances, mea-surement errors and the sensorresolutions.

Generally speaking the parameter estimation methods are very efficient in the sense that given enough parameters they will always give a model which fits a particular measurement very well. Since the methods are fairly complex because they depend on an

interplay of the model structure, the

experimental conditions and the parameter estimation method it is not easy to acquire a detailed understanding of their properties. The value of efficient model validation procedures can therefore not be overemphasized. To demonstrate this we will give an illustration of what happened when an output error method was used to estimate stability derivatives from a 200/200 zig-zag test on a

Mariner class ship, the USS Compass Island. The derivative _{was estimated to 0.00 14}

which is far away from the value 0.0 116 given in literature. The value_{could be accepted}

if the user did not have a priori knowledge. An investigation of the residuals will, however,

immediately reveal that the residuals are far from white. The data thus contradicts one of the basic assumptions required for the output error method. Similarly Akaike's criterion

will directly show that the model with_{process noise is significantly better. In many cases}

it _{is thus very important to have a procedure which allows for process noise. This may}

perhaps explain the difficulties encountered by Sukselainen _{(1975) in his use of the} output error method.

To obtain reliable model parameters it is thus important to model the process noise.

This is done throughout this report. Since output error methods are so common in

literature the calculations have also been performed for models without process noise. See Chapters 6,8,9 and 10.

4. PARAMETER ESTIMATION PROGRAM LISPID

A general program for Linear System Parameter IDentification (LISPID) has been

developed for the computer UNIVAC 1108. It has also been implemented on IBM 360. The program performs the computation of the maximum likeithood estimates described

in Chapter 3.

The program LISPID which is written in FORTRAN consists of 52 subroutines.

Including comments the program size is 9 200 statements. The program without any data storage requires a core of 64 k cells on the UNIVAC1108, if no segmentation is used, it requires a core of 25 k cells using segmentation and overlays. Additional memory space to store data is required.

4.1 General Features

The parameters, which are to be estimated, can enter the system equations in many

different ways. The dependence of these parameters must be supplied by the user in a special subroutine. The models may be given in continuous or discrete time form, and they may be time varying. It is possible to use only measurement noise in the models as well as both measurement and process noise. The former alternative is equivalent to an output error identification method.

The sampling rate may be uniform or varying, and the measurements may be instant-aneous or integrating. The maximum likelihood estimates of the parameters are ob-tained by minimizing, in a certain sense, the one-step prediction errors. It is also

possi-ble to minimize the prediction errors an arbitrary number of time steps ahead. The

program LISPID is described in some more detail in Astrom, Norrbin, Kallstrom and

Bystrom (1974), and Källström, Essebo andAstrom (1976).

4.2 Numerical Aspects

The maximum of the likelihood function is found by an optimization routine. Since it is extremely tedious to compute the gradient analytically, optimization techniques using

only the values of the loss function have been tried.

Two different algorithms have proved to

be the most suitable for this kind of

problems. In the first one the gradient is computed numerically using finite differences. Then a quasi-newton method is applied to find the optimum (see Fletcher (1972)). The other algorithm does not use numerical gradients, but gets information about the loss

function by a special search pattern (see Brent (1973)).

The actual computations are far from straightforward, and although a fast computer, a UNIVAC 1108. is used, the execution times often become rather long, especially when the observations are not equally spaced in time. The execution time also depends on the complexity of the model, the number of measurements, the number of parameters to be estimated and the accuracy of the initial parametervalues.

4.3 Ship Models

The program LISPID has been supplied with a special subroutine to be used for identi-fication of ship steering dynamics. Parameters @ of three different models of the ship

+

+ [

(t-Ti

_j

dt+dw

r(tk)

a1 0131 a1 011012

oi:j

0 0 0 0

a

₂

La

₁₂

a

-La

2

22

a2 0 o 1/a1 o 0

[ (tk

TD) L U' Model 1 (cf(2.5)) L 01

-0

L2 dv

10

V

La

VU6 09 v (t)

L0

V2 L2 0 dr

La

_VV7 L

V08

09010 r (t)

dt +

0 0 1

d'

0 1 0 _{' (t)} 0 0 0 0 1/a

+ e(t,)

(tk) (tk) (tv) +

k0,1,.

,N 0 ₀₁₅ 0 ₁₆ 0 ₀₁₅

w is a Wiener process with incremental covariance R1 dt, where

The measurement errors { e(t) } are assumed to be independent and gaussian with

Note that it is not possible to use the measurement signals v1 and v at the same time,

which explains that the order of R2 is 4 and not 5.

The initial state is given by

025 ₂

cs0

_{1 26}

1 27

-and the initial covariance matrix of state estimate errorsis

0 I 0 0 28 31 32 P(t0) = 03 ₁₀₂₉₁ 033 032 033 030 R1 = 1018 I sin 0 20 sin 0 20 0 0 0 1015 I I0 I 10 I 19 0 I0

II09 l

0

zero mean and covariance

= R2, where 0₂₁ 0 0 10 I 22 0 0 0 0 0 0 0 I 23 0 0 0 0 10241

Note that it is possible to estimate a subset of the 34 _{parameters of the model. The} other parameters are then given arbitrary fixed values. It is not necessary to have

measurements of all five output signals of the model, but attention has to be paid to the

identifiability discussion of Chapter 3.

Hydrodynamic derivatives normalized by use of either the "prime" system or the

"bis" system can be estimated in model 1. If the values of the four derivatives 01, 02, 03 and 04 are given in the "prime" system, all other derivative estimates are obtained in the "prime" system, and analogous for the "bis" system.

Model 2 V

-TO1

1 0 V2

- _02

0 1 + 'r3 dx

-3j

I

L V2 04 07 V4

a1rO6

09 0 0 dt + dw The time lag TD is computed as

TD =

T-T

= Tlsin 0341

where I is the sampling interval.

UI is an artificial unit step input signal to make it possible to estimate the bias

parameters

013 -

. The wind influence during the experiment is described by the parameters 09 and 0, (see Chapter 2). _i and 2 are conversion factors from degrees to radians and from m/s to knots, respectively. The dimensions of inputs, states and measurements are shown below.

deg v1 _knots UI

-

v2 _knots rn/s v knots F _{rad/s r deg/s} rad 1i deg dv dx 2 v (t) x2(t)

dt +

x3(t)

R1 = R = 0 2 17 P(t0) = TD =

T-r

N

10111 1012 I sm ₀₁₄ 019 020

r = TIsinO

₂₇ 018/a2 024 25 10221 026 026 10231

v(tk) =

a2 0 JO 11 0]

(tk)1

x2tk)j

+ o

o]

x3 (tk)

k=0,1,...,N

6 _(tk

_{_TD)]}

[

UI sin 04 0₁₁ 110₁₂

+ e

(tk) 0 15 0 16 0 15 016 10211 024 025

The states x2 and x3 are linear combinations of the original states . f and i. and they

have the dimensions of rn/s2 and rn/s3 resp. The transfer function relating the sway

velocity i to the rudder angle, measured in radians, is given by

05s

+

_;06

L 0452 +

G (s)

= _(4.1) 2 2 3 s3 + 2 _{+ 0}

s +

₀₃

Li

L 2 L

or, if the wind parameters 03 and 06 are equal to zero.

G2(s)

-O4S + 05

K (1 + sT)

s2 +

.1

Li

0

S+T2O2

(1+sT )(1+sT )

1 2

The two transfer functions can be compared to (3.5) and (2.7).

(4.2) Model 3 dx I di

di

+ = v3

c-1205

a12

0 0 0 1

-*;0

0 1 06 04 0 7 0

_L0

L3 L2 02 0 [o

(t-TDjl

Ul

j

dt+dw

x1(t) r(t) i (t)

dt +

r = I

kin

0 24

TD = Ir

The state x1 is a linear combination of the original states _{i and i, and it has the}

dimension of 1/s2. The transfer function relating ,1i to is given by

V2 V3

r

04 S +

G1 (s) = _(4.3)

s3 + *0s2

+2

02s ±

L

or, if the wind parameter fl3is equal to zero,

G1(s) = G6(s) = v2 v3

j

0 S + 0S K(1+sT3) V V2

s(1+sT )(1+sT )

s [s2 + -o s + - 0 J

_Li

L 2 (4.4) Cf(3.3) and (2.6).

The transfer function relating r to & is obtained by multiplying (4.3) or (4.4) by s, e g

K (1+sI)

G3(s) = (4.5)

(1+sT1)(1+sT)

Nomoto's first order model can be identified, using model 3, by giving the wind

parameter 03 and the parameters 02 and 05 the fixed value zero. The following transfer function, relating r to ö, will then be obtained (cf (2.8)):

v2 0 L

S + *01

K

I + sT

(4.6)

A summary of the three ship models, which are described in this section, is also given