July 1982 LIDS-R-1220 REAL TIME ESTIMATION AND PREDICTION OF SHIP MOTIONS USING KALMAN FILTERING TECHNIQUES by Michael S. Triantafyllou Marc Bodson Michael Athans The research was conducted at the M.I.T. Laboratory for Information and Decision Systems with support by the U.S. National Aeronautics and Space Administration Ames Research Center, under Grant NGL-22-009-124. Laboratory for Information and Decision Systems Massachusetts Institute of Technology Cambridge, MA 02139 2 the various modeling parameters on the rms error is assessed and simplifying conslusions are drawn. The models developed are used to predict the motions a few seconds ahead. An upper bound for prediction time of about five seconds is established, with the exception of roll which can be predicted up to ten seconds ahead. The effect of noise and modeling errors on the rms prediction error is investigated in detail. 3 TABLE OF CONTENTS Table of Figures .. .... ........ 4 Tables ...................... 8 Introduction .................... 9 Overview ..................... Chapter 1: 11 EQUATIONS OF MOTION . . . . . . . . . . 13 Definitions . . . . . . . . . . . . . . . . Simple Derivation . . . . . . . . . . . . . . Strip Theory ................. 13 14 18 Relation between Added Mass and Damping . . . . 19 Speed Effects .............. Frequency of Encounter . . . . . . . . . . . . Equations of Motion . . . . . . . . . . 19 20 22 Heave - Pitch Approximation . . . . . . . . . . Sway - Roll - Yaw Approximation . . . . . . . . 24 28 Chapter 2: SEA MODELING . . . . . . . . .. 51 Chapter 3: DERIVATION OF THE STATE-SPACE EQUATION . Heave Pitch Model . . . . . . . . . . . . . . Sway - Roll - Yaw Model . . Chapter 4: KALMAN FILTER AND SILUMATION . . . . . . Heave - Pitch Motions 60 60 64 69 . . . . . . . . . . . . . 69 Sway - Roll - Yaw Motions . . . . . . . . . . . 71 Chapter 5: SHIP MOTION PREDICTION . . . . . . . . Theoretical Background . . . . . . . . . . . . Reduced Number of States . . . . . . . . . . . CONCLUSIONS . . . . . . . . . . . . . . . . . . . . ACKNOWLEDGEMENTS . . . . . .. . . . . . . . . . . . 90 90 96 115 117 REFERENCES . Appendix 1: Hydrodynamic Theory . Appendix 2: Hydrodynamic Data . . . . . . . . . . . 127 Appendix 3: Simulations 129 Appendix 4: Computer Program Listing ................ . . . . . . .. . ... 121 ... . . . . 118 .. 132 4 TABLE OF FIGURES 1.1 Ship Reference Systems 1.2 Wave Profile 1.3 Hydrodynamic Problem 1.4a Radiation Problem 1.4b Diffraction Problem 1.5 Strip Subdivision of Ship Hull 1.6 Effects of Ship Speed 1.7 Angle of Incoming Waves 1.8 Heave Exciting Force on a Rectangular Barge of Length L, Beam B 1.9 Effect of Ship Speed on Heave Force and Pitch Moment at Head Seas 1.10 The Equivalent Wavelength for Heave and Pitch is X/cos 1.11 Phase Difference Between Heave and Pitch 1.12 Wave Force and Moments in Sway, Roll, and Yaw at 60 ° Angle and Zero Forward Speed 1.13 Sway Force Versus Frequency for the DD963 Destroyer and Its Finite Dimensional Approximation (dotted line) for Zero Speed and 90 ° Angle of Incidence 1.14 Roll Moment for DD963. 1.15 Yaw Moment for DD963. 1.16 Sway Force Versus Frequency for the DD963 Destroyer and Its Finite Dimensional Approximation (dotted line) Same conditions as in 1.13. Same conditions as in 1.13. for Speed U=15.5 ft/sec and 45 ° Angle of Incidence 1.17 Roll Moment for DD963. 1.18 Yaw Moment for DD963. 1.19 Sway Transfer Function for the DD963 Destroyer and Its Finite Dimensional Approximation (dotted line), Speed U = 15.5 ft/sec and 450 Angle of Incidence Same conditions as in 1.16 Same conditions as in 1.16. 1.20 Roll Transfer Function for DD963. in 1.19. 1.21 Yaw Transfer Function for DD963. in 1.19. Same conditions as Same conditions as for 5 TABLE OF FIGURES (continued) 2.1 Typical Wave Spectrum with Two Peaks 2.2 Frequency Content of to Storm Duration 2.3 Bretschneider Spectrum at Various Wind Speeds (Fully developed seas) 2.4 Bretschneider Spectrum for H 1 / 3 = 10 ft and Modal Frequency 0.72 rad/sec (sea state 3) and Its Finite Dimensional Approximation (dotted line) 4.1 Results of Heave Simulation and Its Kalman Filter Estimate (dotted line), Using Accurate Model at U=21 ft/sec and 00 Angle of Incidence and in sea state 5. 4.2 Results of Pitch Simulation and Its Kalman Filter Estimate. Same conditions as in 4.1 4.3 Results of Sway Simulation and Its Kalman Filter Estimate (dotted line), Using Accurate Model at a Wave Spectrum According U=15.5 ft/sec and 450 Angle of Incidence, and in sea state 5. 4.4 Results of Roll Simulation and Its Kalman Filter Estimate. Same conditions as in 4.3 4.5 Results of Yaw Simulation and Its Kalman Filter Estimate. Same conditions as in 4.3 4.6 Results of Heave Simulation and Its Kalman Filter The actual wave spectrum Estimate (dotted line). model frequency is 0.52 rad/sec, while the value used in ;the Kalman Filter is 0.72 rad/sec. All other parameters as in 4.1. 4.7 Results of Pitch Simulation and Its Kalman Filter Estimate. Same conditions as in 4.6. 4.8 Results of Sway Simulation and Its Kalman Filter Actual W =0.52 rad/sec, while Estimate (dotted line). in Kalman Filter Wm=0.72 rad/sec. All other parameters as in 4.3. 4.9 Results of Roll Simulation and Its Kalman Filter Estimate. Same conditions as in 4.8. 6 TABLE OF FIGURES (continued) 4.10 Results of Yaw Simulation and Its Kalman Filter Estimate. Same conditions as in 4.8. 4.11 Results ,of Sway Simulation and Its Kalman Filter Estimate (dotted line), using noisy measurements (light line) when the roll measurement is disconnected. Same other conditions as in 4.3. 4.12 Results of Sway Simulation and Its Kalman Filter Estimate (dotted line), using noisy measurements (light line) when the actual angle of incidence is 600 and the value used in the Kalman Filter is 45 ° . Same other conditions as in 4.3. 5.1 Heave Simulation Results and Its Prediction (dotted line starting at t=40 sec) for U=21 ft/sec and q=0 °, and in sea state 5. Perfect state knowledge is assumed. 5.2 Pitch Simulation Results and Its Prediction. conditions as in 5.1. 5.3 Sway Simulation Results and Its Prediction (dotted line starting at t=40 sec) and in sea state 5. 5.4 5.5 Same for U=15.5 ft/sec and ~=45 ° , Perfect state knowledge is assumed. Roll Simulation Results and Its Prediction. conditions as in 5.3. iYaw Simulation Results and Its Prediction. ditions as in 5.3. Same Same con- 5.6 RMS Prediction Error Over RMS Motion Versus Prediction Time for Heave and Pitch, U=21 ft/sec, W =0.72 rad/sec, sea state 5. 5.7 RMS Prediction Error Over RMS Motion Versus Prediction Time for Sway, Roll, Yaw. U=15.5 ft/sec, ~=45 °, sea state 5. 5.8 Heave Simulation Results, Its Kalman Filter Estimate (up to 40 sec) and Its Prediction Using the Kalman Filter Estimate (after t=40 sec). Same conditions as in 5.1. 5.9 Sway Simulation Results, Its Kalman Filter Estimate (up to 40 sec) and Its Prediction Using the Kalman Filter Estimate (after t=40 sec). Same conditions as in 5.3. 7 TABLE OF FIGURES (continued) 5.10 Roll Simulation and Prediction. as in 5.9. Same conditions 5.11 Yaw Simulation and Prediction. in 5.9. 5.12 RMS Prediction Error Versus RMS Motion Versus Prediction Time for Sway. Actual W =0.52 rad/sec, used Wm=0.72 rad/sec All other conditions as in 5.3. 5.13 RMS Error Versus Prediction Time for Roll. as in 5.12. Same condition 5.14 RMS Error Versus Prediction Time for Yaw. as in 5.12 Same conditions 5.15 RMS Prediction Error Over RMS Motion for Heave and Pitch. U=21 ft/sec, 4=0 °, sea state 5. In the prediction model the nonminimum phase zero have been omitted. Same conditions as 8 TABLES 2.1 Sea spectrum coefficients 3.1 Matrices used in the numerical applications 4.1 Poles of the heave, pitch model 4.2 Poles of the heave, pitch Kalman filter 4.3 Sensitivity of the rms error for sway, roll, yaw to changes in the model parameters 4.4 Poles of the sway, roll, yaw model 4.5 Poles of the sway, roll, yaw Kalman filter 9 INTRODUCTION The present study started as part of the effort directed toward designing an efficient scheme for landing VTOL aircraft on destroyers in rough seas. [14] A first study showed a significant effect of the ship model used on the thrust level required for safe landing. In a landing scheme therefore it would be desirable to have accurate ship models capable of providing a good real time estimation of the motions, velocities and accelerations of the landing area, resulting in safer operations and with reduced thrust requirements. The modeling is quite complex and a substantial effort is required to reduce the governing equations to a finite dimensional system of reasonable order. The study contains a first chapter on the equations of motion as derived from hydrodynamics, their form and the physical mechanisms involved and the general form of the approximation. The second chapter describes the modeling of the sea, which proved to be a crucial part of the overall problem. The third chapter describes the derivation of the state-space equations for the DD-963 destroyer. In the fourth chapter the Kalman filter studies are presented and the influence of the various parameters is assessed. 10 In the fifth chapter the feasibility of predicting the ship motions a few seconds ahead in time is studied within the present formulation. Finally the appendices provide the characteristics of the destroyer, hydrodynamic information and some computer programs used. OVERVIEW The real time estimation of the rigid body motions, velocities and accelerations of a vessel in rough seas requires accurate modeling of the wave exciting forces and the hydrodynamic coefficients of the ship. The wave forces are obtained after an integration over the ship hull of the pressure forces, so that their evalua- tion requires a seakeeping program, while their magnitude and phase represent clearly an infinitely dimensional system with non-minimum phase characteristics. The complexity of the resulting equations is due primarily to the wave formation as the vessel moves, which is a mechanism of energy dissipation and additionally it introduces memory effects. The wave spectrum contains a rather narrow band of frequencies so that an efficient approximation of the ship characteristics can be achieved within this frequency band. A DD-963 destroyer was used as the basis for the present study. First the geometric and mass properties of the vessel were analysed by the M.I.T. Ocean Engineering Department Seakeeping program and its hydrodynamic forces and coefficients were obtained. Subsequently a finite dimensional approximation was fitted in this data within the wave frequency range. Two groups of ship motions were distinguished, the heave-pitch and the roll-sway-yaw sets of motion, which up to the first order are uncoupled to each other. The parameters of the approximations are four: 12 The speed of the vessel, the direction of the waves, the significant wave height and the modal frequency of the wave spectrum. These models were used to estimate the ship motions, velocities, accelerations using noisy measurements of the motions. The Kalman filter designed for this purpose gives very good results when a relatively accurate estimate of the modal frequency of the spectrum is available. The modal frequency was found to be the most significant parameter in the overall scheme since it influences the estimation error significantly and is the most difficult to estimate. The ship speed and the wave heading are important parameters also,but can be estimated easily and accurately. The double peak spectrum, i.e. seas containing swell also, require separate treatment, because the low frequency peak is hard to estimate, while its influence is quite important. The predictability of ship motions has been investigated within the frame of the present study. First perfect state information is assumed and by propagating the prediction error covariance from zero initial value it has been established that within 25% rms error over rms motion, the prediction time is about five seconds for all motions with the exception of roll which can be predicted up to ten seconds ahead. Simulations confirmed these results. The effect of noise and modeling errors is to reduce the prediction time. Omission of the non-minimum phase zeros has a particularly pronounced effect. In summary, the approximations described in the sequel provide a good model of the quite complex ship equations of motion within the wave frequency range. The derived models can be used for a real time estimation and prediction of the ship motions and other responses using Kalman filter techniques. Computer programs have been prepared that provide the required model matrices once the parameter of the problem has been specified. 13 Chapter 1: EQUATIONS OF MOTION Definitions The rigid body motions of a ship in six degrees of freedom are shown in Figure 1.1: We define the xlzl plane to coincide with the symmetry plane of the ship, with the zl axis pointing vertically upwards when the vessel is at rest, and the Y 1 axis so as to obtain an orthogonal right-hand system,while the origin need not coincide with the center of gravity. The X oY oZo system is an inertial system with X Yo fixed on the undisturbed sea surface, o while the x y z system is moving with the steady speed of the vessel (i.e. it follows the vessel but it does not participate in its unsteady motion). Then the linear motions along the xl, Y 1 , Zl axes are surge, sway and heave respectively. In order to define the angular motions, we normally require Euler angles, in the present case, though, we consider small motions so that the tensor of angular displacements can be replaced by a vector of small angular displacements, which are roll, pitch, yaw around the x 1, Yl, Zl axes respectively. The characteristics of a ship are its slender form, i.e. L/B>>i, L/T>>l, where L is the length, B the beam and T the draft. Also, the ship is symmetric about the xz plane and near symmetric about the yz plane. For this reason I I yz xy =I =I zy yx = 0 =0 14 The value of Ixz is typically small compared with I, Iyy. The justification of using the linearity assumption is as follows: The excitation consists of wave induced forces, which include fluid inertia forces and hydrostatic forces. It is well established that the wave height to wave length ratio is small, since at a typical upper value of 1/7 the wave breakes and loses all its energy [15] (Figure 1.2). As a result, the major part of the wave force is a linear function of the wave elevation and can be obtained by a first order perturbation expansion of the nonlinear fluid equation, using the wave height to length ratio as the perturbation parameter [15]. The wave spectrum, as will be shown later, has a frequency range between typically 0.2 and 2 rad/sec. Given the large mass of the vessel, the resulting motions, within this frequency range, are of the order of a few feet, or a few degrees, so that the equations of motion can be linearized. The only motion that requires attention is roll, because due to the slender form of the ship, the rolling motion may become large, in which case nonlinear damping becomes important. Simple Derivation We derive the equation of motion for a simple two dimensional object to demonstrate the overall procedure. Let us assume that we wish to derive the motion of a two dimensional cylinder subject to wave excitation, allowed to move in heave only (Figure 1.3). 15 The incoming wave of amplitude ao and frequency w 0 will cause a force on the cylinder, and, therefore, heave motion. Due to the linearity of the problem, the following decomposition can be used, which simplifies the problem considerably. (a) Consider the sea calm and the ship forced to move sinusoidally with unit heave amplitude, and frequency wo, and find the resulting force. (b) Consider the ship motionless and find the force on the cylinder due to the incoming waves and the diffraction effects (diffraction problem). (c) In order to find the heave amplitude, within linear theory, we equate the force found in amplitude, with the force The force in (a) times the found in (yet unknown) heave (b). (Figure 1.4) (b) can be decomposed further for modeling purposes, again due to linearity: One part is due to the un- disturbed incoming waves and the other part due to the diffracted waves. The first is called the Froude-Krylov force and the second the diffraction force. excitation force The force in The total force is called the [15]. (a) due to linearity can be also decomposed: The first part is simply the hydrostatic force. The second part is the dissipative force, caused by the fact that the refraction waves carry energy from the ship to infinity. For this reason, we define a damping coefficient B so that the dissipative force -Bx where x is the heave velocity will be -Bx where x is the heave velocity. The third part is 16 is an inertia force, caused by the fact that the heaving ship causes the fluid particles to move in an unsteady motion so that we define and "added" mass A and the inertia force becomes -Ax If we denote the undisturbed in- with x the heave acceleration. coming wave elevation amidships as n (t) n(t): (l} · = aoci Where the real part of all complex quantities is meant, here Then the excitation force will be and in the sequel. . F =· Feiwot (2) Where F o is complex (to take into account the phase difference with respect to the wave elevation), 'and the equation of motion becomes: (3) Mx = F - Ax - Bx - Cx Where M the mass' of the cylinder; the motion is also sinusoidal' so with x, complex: (4) x(t) = xoeit A very important remark is that r,A,B depend on the frequency of the incoming wave w0 . This can be easily understood by the fact that at-various frequencies the heaving cylinder will produce waves with different wavelength. We rewrite, therefore, equation (3) as: -Mx, eiwt - By dropping e i (-[M + A( FO(w )a~+ [A(w)w ° t, ]o + - i B(w) - C] x 0 }e iw° t we can rewrite equation (3a) as: WoB() + C} X = Fo ( o ) ao (3b) 17 The motion of a cylinder in water, therefore, results in an increase in the mass and damping term. Equation (3b) is used because of its similarity to a second order system, it is strictly valid, though, only for a monochromatic wave. Ultimately, we wish to obtain the response in a random sea, so equation (3b) must be extended for a random sea. This can be done by obtaining the inverse Fourier transform of (3a), i.e. co co K (t - T) x (T)dT+ f K (t - t) x (T) dt + ~ . Y ; -~~~~~~~~~~ + C x(t) = I Kf(t- Where K ,KX,Kf iwB(m) and F, (X) T)n(t)dT the inverse Fourier transform of -w 2 [M + A(X)], respectively. is denoted by n(t). The random undisturbed wave elevation Equation (5) is not popular with hydrodynamicists, because the effort required to evaluate the kernels K a, Ku, Kf is by far greater than that required to find the added mass, damping ana excitation force. (5) For this reason, equation (5) is rewritten in a hybrid form as follows: - EM + A(w)] X(t) + B(i) x(t) + C x(t) = r(,w)(t) (6) This is an integro-differential equation (or differential equation with frequency dependent coefficients), whose meaning is in the sense of equation (5). 18 Strip Theory The evaluation of A(w), B(w), complex geometries, such as F(w) is not an easy task for the hull of a ship. The hydrodynamic particulars can be found in a later section, but we can give a simple description here of a technique used to simplify the derivations: [15], [17] The ship can be divided in many transverse strips as shown in Figure (1.5) Due to its elongated form and for high frequencies, each strip has small interactions with the other strips, except near the ends. Usually these end effects are smll, so that instead of solving the overall three dimensional problem, we can solve many two dimensional problems sum up all -the partial results. (one for each strip) and For the case of heave, for example, if A(m,x) B(w,x) are the added mass and damping in heave of a strip at location x, then .2I A(w) A(w,x) =f dx (7) -L/2 L/2 B(w) - / B(w,x) dx (8) The strip theory has larger errors at smaller frequencies. It so happens, though, that at small frequencies the hydrostatic forces are predominant, so that the motion error is quite small. Comparison with experiments has shown that for slender ship 19 configurations, the strip theory provides very good predictions [15],[17]. Relation Between Added Mass And Damping The added mass and damping coefficients are not independent of each other, because their frequency dependence is caused by the same refraction waves. If we define T(w) = m 2 [A(w) - B(w) 1 (9) Then T(w) is an analytic function [16 ]. As a result, A(w), B( ), which are real, are related by the Kramers-Kronig in order to describe a causal system. relations, This fact will be used later to obtain a single approximation for T(w) instead of two separate approximations for A(w), B(w). Speed Effects As it can be seen in Figure 1.6 when the ship is heaving with a small angle 9 and at the same time is moving forward with speed 0 U, then a heave velocity results, which is x = Ue. The effect of the forward speed, therefore, is to couple the various motions by speed dependent coefficients. As it can be found in Appendix 1, there are simplified expressions for the added mass, damping and exciting force with a parametric dependence on the speed U. Then 20 expressions greatly facilitate the evaluation of the ship motions. Frequency Of Encounter An additional effect of the ship speed is the change in the If the incident wave has a frequency w frequency of encounter. and a wave number k, then the frequency of encounter we is m e = w + k U cos ( (10) Where ( is the angle between the x axis of the ship and the direction of wave propagation (Figure 1.7). In deep water, the dispersion relation for water waves is 2 = kg (11) so that we can rewrite (10) as We = +- U cos ( (i2) A very important consideration in the difference between frequency of encounter and wave frequency is the following: The ship motions due to linearity will be of frequency we so that the refraction waves are of frequency wA and the added mass and damping can be written as A(w ), B(w ). e e The amplitude of the exciting force though, consists of the Froude Krylov part which depends on w and the diffraction and speed dependent parts which depend on we. is again We t, i.e. ~ - --- ~ ~--_e The time dependence 21 F(t) = aoF(w, we,v) ei W t (13) with a o the incident wave amplitude. This is a very crucial observation and can cause significant errors if not taken into account. 22 Equation of Motion Following the notation of Appendix 1, we write the equations of motion. It should be noted that, due to the slenderness of the ship, the surge motion is left out as a second order motion. is in agreement with experiments [171. This Within linear theory and using the ship symmetry, the heave and pitch motions are not couDled with the group of sway, roll, yaw motions. This is not to imply that the motions are independent, because they are excited by the same wave, so there is a definite relation both in amplitude and phase. (1) Heave - Pitch Motions {[Mo ] o I [+A A 33 5 3 A 3 s ]} X -v Ass + r B 3 , B3 5 1 x Bsa [C 3 C 3 s ] C 5 3 Css5 -v (2) s5 F 3 F5 + ' -V ] n (14) Sway - Roll - Yaw Motion M24 2M M26 -x B22 B2 TI 4x+ -xz X B2j 6 B42 B 4 B62 B4 B6 ! A 22 A 2 f A2 6 A42 A44 A |.O O O x s66 -u + -u (15) O Fe LF6_ 46 O O Where A ij,B ij,Ci.. efficient matrices; the added mass, dampina, F. V X the exciting forces; n = {X3 X5} {X 2 X ,4} hydrostatic cothe wave elevation; (16) (17) T U The frequency and velocity dependence is not written explictly, but is understood, as described in the previous sections. 24 Heave - Pitch Approximation We start with the heave and pitch motions approximation. As it is obvious from equation (14), it involves two stages: (a) Approximation of the exciting force (b) Approximation of the added mass and damping coefficients Data are provided by the hydrodynamic theory for both components and within the wave frequency range. A. Exciting Force Approximation Figure 1.8 shows the exciting heave force on a box-like ship [25]. This information is important to demonstrate several zeros of the amplitude of the heaving force. Figure 1.9 shows the ampli- tude and phase of the exciting force on a destroyer, where, again, the same zeros appear, accompanied in the phase plot by jumps in the phase. The transfer function between the wave elevation and the heave force cannot be represented as a ratio of polynomials of finite degree as evidenced by Figure 1.9. for the pitch moment. Similar plots can be obtained Within the wave frequency range, though, only the first zero is important, while the remaining peaks are of minor significance. This is not true for other types of vehicles such as the semi-submersible, but for ships it is valid for both heave force and pitch moment, so it will be used to simplify considerably the modeling procedure. As it was mentioned before, the exciting force changes with frequency w e , but its amplitude is determined on the basis of the frequency w. The following variables must be included in an 25 appropriate modeling of the exciting forces F 3 (t, as, Fs5 Ct, a., (1) frequency X (2) speed V (3) wave angle Q O, U) = F 3 ,U) = (w, 4)a o eiwet {Fs5 ( ) + i (18) f3W ) e e where a o the wave amplitude, f 3 the heave diffraction force. tions (18), (19) Equa- (19) shIcw that the heave force does not depend on the ship speed, whereas the pitch moment does, in a linear fashion. In order to approximate F3 (w,), Fs(C,c),-f3 (we,) we use the plots in Figure 1.9 as well as Figure 1.10, which show the approximate influence of the wave angle on the excitation force. In order to model the DD963 destroyer, the M.I.T. five degrees of freedom seakeeping program [27 j was used to derive hydrodynamic results. The following model was derived to model shape of the heave force at V = 0 and p = 0 (no speed, head seas) F3(s) = a n [1 + 2J s . + 2 2 (20) Where J = 0.707,al a constant to be determined from hydrodynamic data, n the wave elevation and wa the corner frequency. Remembering the analysis above concerning the dependence of the force on w and Figure 1.10, we can derive 26 /w 2irg = ca + B + cos U cos 2 Lcoso + B 2 1) where L is the ship length, B the beam. Before we establish a relation similar to (20) above, we have to discuss Figure 11, where it is shown that for long waves, the heave force and the pitching movement are 900 out of phase. This means that the transfer function between heave and pitch is a non-minimum phase one, because the amplitude is constant, while the phase is 90 ° . We choose to attribute the non-minimum phase to pitch. Also, the pitch angle tends to the wave slope for large wavelengths, so the pitching moment F5 can be written as a = 2 1i- s/wO ~s 1 + s/w (22) cos4 [ 1 + 2J X + 2T a Wa Where a 2 a constant to be determined, wa is the same (for simplicity) as in equation (21) and Wo is an artificial frequency to model the non-minimum phase. It will be chosen to be eaual to the wave spectrum modal frequency, so we defer the discussion until the corresponding section. B. Added Mass and Damoing By using equation (9), we can rewrite the equations of motion as _? s2 + "S2 + C33 T35 + UsT 3 3 + C 3 T53 s2 - UsT 3 3 + C3 5 Is2 + T 55 5 s-- U 2 T 3 3 + C55 r (23) 27 Here we construct a simplified model where i, Tij= Aijwith Aij,Bij to (24) be evaluated from the hydrodynamic data. Then eaua- tion (23) can be written, after we define O Y1 = X3 y3 = xS Y = X3 Y2 y4 = xs (25) = {y1,Y2,Y3,Y4 } in the form: FY] LY rA 33 LAs53 A3 C33 B ]5 A5s C 33 B 3] C' 5 B s5 B 53 3 B35F 5 Y + F (23a) (23a) F where C3 5 = C 3 5 + U A 3 3 C53 = C3 C's = Css - 5 - U A3 3 2 X A (26) 33 28 Sway-Roll-Yaw Approximation Next we approximate the sway-roll-yaw group of motions, which is uncoupled to first order from the heave-pitch group of motions. Again a two stage approximation is required, i.e.: (a) Approximation of the exciting force (b) Approximation of the added mass and damping coefficients Data are obtained by using the sea-keeping program. A. Exciting force approximation The same infinite-dimensional form is obtained for the exciting force as seen in Figure 1.12 (17abc) for all three motions, as in the case of pitch and heave. Again, within the wave frequency range, a finite dimensional approximation can be achieved, and of reasonably small order. The important fact that the exciting force depends on the wave frequency rather than the frequency of en- counter, is used, while the following three quantities define the exciting force amplitude and phase (1) Wave frequency (2) Speed V (3) Wave angle In Appendix 1 the strip theory approximation of the sway, roll,yaw forces can be found. Using the M.I.T. five degrees of freedom seakeeping program the following finite dimensional approximation was found in case of V = 0, ~ = 900 for sway, roll and yaw 29 2 A2 F2 S (S) = (27) j 5 + 2J 2 + 1 2 F4 (S) S A4 S = A S 2+ 2 2J4 + 1 S2 A6 F6 (28) (S)= (29) + 2J6 S + 1 w where = 0.65 2 J 2 = .5 m4 = 0.85 J4 = .3 W6 = 0.85 J6 = .3 (30) and A 2 ,A 4 ,A 6 are obtained from hydrodynamic data. We redefine the value of z, W4,W 6 such that 2 it will be valid for angles speed other than 0: W. where = (W ~ other than 900, and V + -2 V cos4) sing (31) j = 2,4,6 and w. is given above. It should be noted that the sway,roll and yaw forces 0 are proportional to the wave slope, i.e. 90 out of phase with respect to the wave amplitude. This means that they belong to the same group with pitch, and the same nonminimum phase transfer function: s - s + Wo Wo must be used for all three of them, when the total system (all 5 degrees of motion) is considered. 30 B. Added Mass and Damping The amplitude of the transfer function between the wave elevation and the rolling motion has a very narrow peak so that the coefficients can be approximated as constant [17]. Using Appendix 1: A4 A 0 A A =A44 42 44 B AO 24 42 0 44 Bo 24 V A B 46 46 using the value of A0 +- B O46 -VAO24 46 X B 2 24 at the roll peak. (32) (32) It should be noted that roll involves a significant nonlinear (viscous) damping, which is approximated by introducing an additional "equivalent" damping coefficient B* [ 3 44 31 £irmialar_-, we calculate the sway, yaw coefficients at the same frequency: A A 26 B 0 0 B B 26 A B A B V A22 22 26 AO 62 2 22 Bo +VA AO +Y 26 66 66 V BD A 26 62 66 _ . 0 66 Be 22 26 = = 22 A0 2 22 + Y BO 22 (33) 0 S2 {Aij+Mij} + S{Bij} + i {Cij} {s2{Aij+M. (F 6 ) 2,4,6 = 2,4,6 (34) where Cij = 0 except for C 4 4, which is the roll hydrostatic constant, i.e. C 4 4 = A-(GM) with A the ship displacement and the (GM) metacentric height. Due to the special form of the matrix C, a zero-pole cancellation results from a direct state space representation of the equation above. After some manipulations, the following representation can be obtained which avoids zero-pole cancellation problems: 32 X = T X +UF 0 where {X2 X = F= T ,X4X 6} X4 {fF 2 , F 2 , , F 6 , F 6} (35) where fF indicates the time integral of F and T = U = {uij , with: ti~ r _!2 = =p _p r- -t r2 2 P211 r t4i _p -2p r 22 21 ,ti2-t .r ,t12 1 31 r2 --P21tll- P23t41 t. =-P = t- Pt 13 -P21 23 tl3 .r = - r22 r3 2 , t42.: _ t21 t23 =t, t43= 2 t22 22 2 .r 32p _ p r2 22 2 -P21t12 32 - 12 t r .r , t44 = _2 P23t42- r2 1. 2 P-P r22 14 2 23 23 p 13 33 P22 rC 43- 22C44 t24 24i =-P 21 t 14 -P2 23 t 4 4 - t3 i =0 except for t32 =1 U.. except for U = 0 2 P-P r2 {tij }and (36) -r1 2 r r r 21 22 .r r. r U14 14 =r- 1 r 13 r U U 44 22 U2 5 where =r 33 = r, 21 U 23 22 =r 41 r 32 : 31 r 21 22 r 32 r U U2 1 = 23 P 22 23 = r , U 22 24 =-P 21 21 U 11 -P 23 U -P U 14 23 U 41 44 = r2 3 R (37) {ri.= [A+M] 1 P M T R=[A+M] =' {rij } {P }= R B (38) __ 33 Figures 1.13 through 1.15 show the actual force and moments versus frequency and the achieved finite dimensional for zero speed and beam seas (4=90°). approximation It should be noted again that the approximation is not as good outside the wave frequency range. Figures 1.16 through 1.18 show the same quantities for speed U= 15.5 ft/sec and 450 angle of incidence. From these figures it can be seen that above 1 rad/sec the approximation is poor, nonetheless no significant wave energy is contained in that range, so the approximation is acceptable. Figures 1.19 through 1.21 show the overall transfer function between the corresponding motion and the sea elevation. 34 z1 heave f yaw ~ Yl z \5way i h I x~~x oY X, Figure 1.1 Ship Reference Systems 35 rrxecwn wafer s- rface N/: >>1i < 4/7 Figure 1.2 Wave Profile cyinActer posiEion fec Problem Hydrodynamic Hydrodynamic Problem 36 ;aN~ %race r~~~s Figure 1.4a Radiation Problem ,ninGm - \ cyt inder posi ion Figure 1.4b Diffraction Problem caiwoWave daceiffr= 37 Figure 1. 5 Strip Subdivision of Ship Hull ..A... or Figure 1. 6 Effects of Ship Speed = a 38 U X-- i I Figure 1.7 Angle of Incoming Waves 1.0 0. 8 0.4 0.2 -0.2 0.05 0.1 02 0.3 1.0 2.0 3.04.0 L/X Fiaure 1.8 Heave Exciting Force on a Rectangular Barge of Length L, Beam B incoming waves 39 1.2 0 - Q X -05 0.25 0 \ - 0 0.25 So 50 T* - S lXzp-- . S: - - F'"I' ~ 0.7S 1.00 1.25 t0 iI - . T ......t\t~o t0 IS e O 0 025 050 075 V" 100 1.25 S 1.25 1. 0.75 050 -- 00a-t \ H e. ... s. '_ 0.25 0.50 075 ... ........ 0.25 0 o 1.00 10-... 1.25 1.50 0 0.25 o50 0.75 . ...... 1.00 1.2S Figure 1.9 Effect of Ship Speed on Heave Force and Pitch Moment at Head Seas 1.50 40 t Figure 1.10 The Equivalent Wavelength for Heave and Pitch is X/cosO 90 ., -------- ° IX zero heave f//maximum pitch mean water maximum heave zero pitch Figure 1.11 Phase Difference Between Heave and Pitch 41 SWAY 1, I j~s o 1.Q 0.25 _ 1 02_ · ... 00 ao 9 1 t.0 0.75 - 125 - TS0 0 ,i,. 1 1 1 i 1 - 2S -7 o50 .GS 1.00 1.25 tS0 w.7IS L. , . 02~0751 I II 1I 0.10 1 _tsoi ,, 11 __ I ,I I i i' I -so o 00' I I o .25 I So0 I .S L 1 100 t2S 50 0 I I I 025 0.50 I I 07S 100 I 1.25 1SO Figure 1.12 Wave Force and Moments in Sway, Roll, and Yaw at 600 Angle and Zero Forward Speed 42 SWAIY FORCE pprox 50.00 45.00 48.00 m 3 3500- a: 30,.8- 25.00 I 28.0 I I X e Xu I , I I"0 ro U1 a -1 'n en n> nr X .n cn Omega (dB) Figure 1.13 Sway Force Versus Frequency for the DD963 Finite Dimensional Approximation U = 15.5 ft/sec and 450 Destroyer and Its (dotted line) Angle of Incidence for Speed 43 ROLL MOMENT 65.08 8.80 45.88 G 0 50..00 . 35.00 ·n Ec9 '( a, 19 ' ' a:_. Omega - (a a - (dB) igure Roll Moment for DD963. Same conditions as .14 in 1.13. 44 YAW MOMENT 85.00 - 80.0/ - 7500 - 70.00 e0.0 0 - f 55.00 - 5 .00 e n , i en to m 6J XA X Omega n , (dB) Figure 1.15 Yaw Moment for DD963. Same conditions as in 1.13. X A 45 SWfRY FORCE 37,50 F 35.00 - 32.50 30,00 27.50 -/ 25.00 22.50 20.08 - 17.50 15 . 00 _ I "I-_ en nJ I 6a I n Omega . I n; n9 cn (dB) Figure 1.16 Sway Force Versus Frequency for the DD963 Destroyer and Its Finite Dimensional Approximation (dotted line), U = 15.5 ft/sec and 450 Angle of Incidence for Speed 46 ROLL MOMENT 68.88 55.880 45 s00 : 35.08 / / 48.008 - 30 .00 ---en1 I nO (Ai~ 63 '" J - " ' I en (A I X (A _/ ro en - (A 6 (An Omega (dB) Figure 1.17 Roll Moment for DD963. Same conditions as in 1.16. 47 YAW MOMENT 85.00 - 80.00 - 75.00 - E 70.00 / -/ C Omega (dB) Figure 1.18 Yaw Moment for DD963. Same conditions as in 1.16. 48 SNAY TRANSFER FUNCTION -180.00 -280.00 -30.00 - 3- -50.00 -0 -60.80 . 00 |_ -.- i a_ XCr l_ a I Cn I I X _ COme(dB) Ome (dB) Figure 1.19 Sway Transfer Function for the DD963 Destroyer and Its Finite Dimensional Approximation (dotted line), for Speed U = 15.5 ft/sec and 450 Angle of Incidence 49 ROLL TRANSFER FUNCTION 30,8800 208,00 10. 00 V1 0.00 0 - -10.80 -280.0 1. -30.00 -480 n 00 1 Ul Xn _ I ~ en " I I I f - u' n( X e t n . no u Un -, e Ome (dB) Figure 1.20 Roll Transfer Function for DD963. Same conditions as in 1.19. 50 YRAW TRANSFER FUNCTION I.00 - 0.00 - -380.00 - -48000 -50.8 0b0 ------'*--_-" - * . 0 X1 ; . * L I J1 _ . . . _, . ...... Ome(dB) Figure 1.21 Yaw Transfer Function for DD963. Same conditions as in 1.19. 51 Chapter 2: SEA MODELING The sea waves are generated by the wind, except for very few cases. The process of wave generation is of importance in model- ing, so we will outline, briefly, a simple theoretical model: When the wind starts blowing over a calm water surface, it contains gust components of high frequency, which cause wavelets on the surface. This is due to the inherent instability of the wave air interface. As soon as the surface becomes rough, a significant drag force develops betweer air and water, which becomes zero only if the average wind speed (which causes the major part of the drag) equals the phase wave velocity. As a result, the steady-state condition of the sea develops slowly by creating waveswhose phase velocity is close to the wind speed. Since the process starts with high frequencies, we conclude that a young storm will contain a peak at high frequency. We usually distinguish between a developing storm and a fully developed storm. As soon as the wind stops blowing, then the water viscosity dissipates the high frequency waves so that the so called swell (decaying seas) forms, which consists of long waves (low frequency content), which travel away from the storm that originates them. For this reason, swell can be found together with another local storm (Figure 2.1). 52 A storm usually contains one peak (except if swell is present when it contains two peaks) and the peak frequency Xm is called the modal frequency (Figure 2.2). Also, the intensity of the storm Beaufort is required, which can be described in a number of waves: Scale, Sea State, Wind Average Velocity, Significant Wave Height. The best is the significant wave height H defined as the statistical average of the 1/3 highest waveheight. For a narrow band spectrum of area Mo IH - 4 F~TM (39) From our discussion on sea storm generation, we conclude that it is important to model a storm by both H (intensity) and w (duration of storm). For this reason, the Bretschneider Spectrum will be used defined as: SM(U) 1.25 4 H2 " M exp{ (40) 1.25 (------) -- 1.25 (40) ( The spectrum was developed by Bretschneider for the North Atlantic, for unidirectional seas, with unlimited fetch, infinite depth and no swell. It was developed to satisfy asymptotic theoretical predictions and to fit North Atlantic data. found to fit reasonably well in any sea location. It was Also, by combining two such spectra, we can model the swell as well. main limitations are unidirectionality and unlimited fetch. Its 53 It was felt, however, that it could provide an adequate description for the present application for open sea. As it has already been mentioned, the forward speed of the vessel causes a shift in the wave frequency to the frequency of now, can be defined for ship coordinates The spectrum, encounter. as follows S S(We) = (41) m W= L~e~ f(W e where -1 t = f(W ) + A -1V/ + Ucos 4m = 9 e (42) 2 u cos4 g A rational approximation was found to (29) subject to (30) in the following form S (W Sa(e) 1.25 4 m-+ Be/o)e ( 3 B( e (43) Wo where Sa () a e the approximate spectrum g B(c), m y(a) functions given in table 2.1 (44) 54 Now a transfer function can be deriven from provide an output with the spectrum in (43) such as to (40) when driven by white It is easy to see that noise. (s/ H o) 2 (s) = s ~~wa + ( + -2J ) 5) ''(4 3 where 1..25 2.25 H2 (46) B(ca) m (47) , = y(ca) Wm J = 0.707 A plot of the spectrum for various wind speeds is given in Figure 2.3, while the spectrum and its rational approximation is plotted in Figure 2.4 = for fully developed seas and H 1 / 3 3 m. Important Remark It is customary to define the power spectral density as the Fourier transform of the autocorrelation R(T) S(w) = f R(r)e (48) dT In wave theory, the spectrum is defined one sided (for positive frequencies only) as follows S (W) 1 (0 ! R(T)e -iLT drT > 0 (49) (49) 55 as For this reason, the relation between the spectrum S(w) required for the present application and the Bretschneider spectrum S (X) is: S(M) = ir S(W) X > O Iar X (50) S1(-W) < 0 Therefore, the intensity of the white noise required for driving the transfer function (33) is multiply the transfer function by X i-T). (or equivalently we can 56 TABLE 2.1 Sea Spectrum Coefficients _ _ .00 .10 20 30 y S (c) (a) .95 38 1.0902 1. 1809 1.8. 61 1,6110 1, 3 ':, 1·271-7- 1 ,2116 .40 1.3626 1 .0765 . 1 I. 1.45:39 , . 168 .70 1.5448 1.6 601 .8845 . 8116 . 80 1.7272 .7498 .90 1.00 1.10 1.8182 1.9095 2. 0008 .6968 .6509 . 0C 1.20 . 30 1. 1 2.0918 2.1833 .2744 5-750 .5434 515 1.60 1. 7 1.80 I . 90 :. 0!0. . 4567 2 .5481 '. 63. 954 2.7306 4I4 ,4454 .4262 ,4085 .8 21 1 718 .3923 57 Pe 0.1 ,. o OA Figure 2.1 sC'') eA s,.as !seae9 (SensedC) t a te i - 04 OA tCp'AmA 2\ \/ 1 em F-igure 2_.2 0.i 4.D(r ca/.) 58 £8.08Q r ;o.7 L 3 )ii _.. L. ] 3 ~Ii~ I 1 It2.,30 .-m~~~~~~~~~~~F~.-Ul$-'08 Ls' j Ct,~ DQ.Q. -2F oN I1 Ct' s.. 88 .- ph in 1/____ .: ? ,, j( Fiue23 Beshede=pcrma aiu d sa) p l ( u l rev sp e s win wind speeds (fully developed seas). 59 BRETSCHNEIDER SPECTRUM AND APPROX. HI/3 = 3.(cm) 1.25 - Approximation d .75/ VE/ >- z *2 **/ a. ul .2500 Up At Cs s C1 FREQUENCY Figure 2.4 t na tA (rad/sec) Bretschneider Spectrum and Approximation uT m 60 Chapter 3: DERIVATION OF THE STATE-- SPACE EQUATION We proceed to derive a state-space form of the equations of motion. Starting with the sea, we can construct the following representation (three cascaded second order systems) 0 1 0 0 -_2 -2J, o 0 0 O 0 0 0 =02 0 0 0< 0 0 0 0 0O O O 0 0x 1 +0 0 0 0 0 0 0 1 0 0 _2 -2Jo -2 n=[E .s/ (51) 1 oo 00o ] x5 or 0 x --s n= = As x C s s + Bw x -s Heave-Pitch Model The following model is derived for the force (some algebra was involved to reduce the dimension of the state). (52) 61 o 1 0 2o -m o -2J o o 0 Xf = 0 0 81 82 0 O O 0 0 1 0 xf '+ 1 - 0 1I a2 cos 0 0 0 0 x 0 8 03 :F in 0 -2Jwo -wo _ 0 or written 0 short 0 Xf = Af xf +Bf (54) = Xf Cf The inertia model can be written as follows: 0 1 1 $2 0 70 0 3 0 X 4 Xm = D2 3 m 0 5s5 0 0 1 $6 B7 I3 x= 88 O 0 D D4 3 (55) 1 0 o 0 1 1 0 Xm 62 where D1 D2 .1= D3 rB1 s2 3 F:I :~ ::a7 |$Bs B a4 D2 DI - = s8 6 D A3 3 A35 s Ass 1 (56) C 33 3 B3 IL(57) D C B D3 4 3 3 5 3 C35 B C5 B 5 35 and we can write in short = A x x -m m -m B + m (58) r~x 3 Xs -C m Xm Now the total model can be constructed as follows: 0 t = At xt + B t W (59) XS= C t X t where A A t s s O t= A B.C zAs O 'B 0 m O0 B C mf A Lo f T 5s B. t O -J (60abc) 63 Data were obtained for the DD-963 destroyer from the M.I.T. Seakeeping Program and are given in Appendix 2. Appendix 4 lists a computer program that produces the At, B t , C t matrices once the ship speed, wave heading, significant wave height and modal frequency were specified. The output can be used directly as input to the LIDS control and filter design package. Table 3.1 provides the numerical values of matrix A for speed 20 knots, angle 0 ° (head seas), significant wave height 10 ft. and modal frequency 0.72 rad/sec (sea state 5). 64 Sway-Roll-Yaw Model The sway, roll, yaw exciting forces and moments are essentially driven by the slope of the sea elevation, which for regular waves equals the wavenumber times the amplitude, or for deep water we can write 2 Slope = -a i.e. in the time domain: 1 Slope = d2 2 dt (61) As outlined in Chapter 1, the fact that only roll has a spring constant causes zero pole cancellation problems, which can be avoided by introducing the matrices T and U described by equations (36) through (38). Il The original equation is in the form I I (M+A) -1 + Bxi + C = F (62) where x1 is of dimension 3: = roll = [ X (63) yaw M is the mass matrix, A and B the added mass and damping matrices respectively, C is the hydrostatic matrix and F the vector of exciting force and moments. S2 X -l = - SCx s( -l Then [M+A] 1 [M+A] -1F -l+A] (64) 65 By letting x T FT = [x2 , X4 , X 4 = [IF 2, , x 6] F2 , F 4 , SF4 , F 6 ] (65ab) and using the T and U matrices of equations (36),(38) we obtain a state space description of (64) without zero-pole cancellation, in the form x = T x + U F (66) The state space representation of equations (27) through (29) is in the form = -2 j j (67) where j = 2,4,6, while fFl Aj-2 J l Fj O j.2.j -- l j (68) For this reason we build a force matrix AF: AF = A2 0 0 0 A4 A (69) A6 and a matrix BF driving the force dynamics with n , i.e. using the sea model, which is exactly the same as in the case of heave, pitch (6 states), so 66 BF(i,j) = 0 except B F (2,2) = -A22 B F (4,2) = -A (70abcd) 4 W4 F B F (6,2) = -A6W Then using the same sea model described in equations (52ab) we obtain the overall model as (16x16) x -t = A (16xl) x t -t + BtW t (71ab) x2 (3x16) ix 4 C x where As At = T BF AF +U TE T BT t =[ BT s Ct = ] I (72abc) 1000 0 0010 ' 0001 Data for all quantities involved are given in Appendix 2, while Appendix 4 lists the computer programs that can produce the matrices At, Bt, Ct once the ship speed, wave heading, significant wave height and modal frequency are specified. 67 Table 3.1 provides the A t matrix for speed 15.5 knots, heading 450, significant wave height 10 ft. and modal frequency 0.72 rad/sec (sea state 5). 68 SW II - ! ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ f~~~~~~~~~~S a . 1~ WW h II) sr~~~~~~~~~~~~~~~L 1~; 9 WI I La W If I1 W I B oCD C) siairdaim~~~~~~~~~~~~~~~~~~~~'n S~~~~~~~~~~~~~~~~~~~~~W W fMnS xWSWWCDSWWQ~WW W ' W L W S~~~~~~~~~~~SI WaW I di NO) C) ' CP CP a ) ~~ ~ ~ ~ ~ I: 0*) C ;' + V~~~~~~~~~~~~~~~~~~~~~~~~~~~O S-S' 5'a x ft W W55 C) . . LfS . . L C)U . . ~~~~~~~~~~~~~+ 4. +W dcd~Mddd~Mddddd CDSSCDS-.-.S W W U · tm Lac , La L CD L. LL iii~~~cce4 LA ~iiiQe22amddMdddd a; 0;M 0 a -II: iaa;;a;cMMv x3 C ~i~ii~ii LI La -~ o CDC 55 a ')AA C') LaLaLa CD S ~ ~~ ~ ~ ~ ~ ~ ~ C8 4 CD -~iiiiii La La N 'C-- .4.)55ff55 D5 S..55 C La ~ 55 ·- N C W S C')S ft-U)iiii CD .LU m WC 4.~~~~~~~~~~ La~i~~ CD A' SS S - SSfSS)t.i D" 0) ----- C La 4. S- A~~ii~i 33s~i~d La SW 0~~~~~~~~~~~~~~~~~LA~W L 0O) 2Mdiddddiddddddr (CC) - ft DCD55 Ph cu L 3, 'r La CD 0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ W E9 E-l~~~~~~~~~~~. La LaW I~~~~~I M N a;~~~~~~~~~~ 53~m toL L~~~~~~~~L La a~~~ 4.1 * 1 .LO X -- U)- s~as~ N~~ 69 Chapter 4: KALMAN FILTER AND SIMULATION The heave-pitch approximation resulted in a 15 state system and the sway-roll-yaw approximation in a 16 state system. Given that 6 states describe the sea, the total system required for 5 degree of freedom motion studies would contain 25 states. If the sea spectrum contains two peaks then a 31 state model is required. The heave-pitch group is not coupled with the sway-roll-yaw group so that the study of each group can be independent. This is not to indicate that in a total design the two groups must remain independent, since they are excited by the same sea. Heave-Pitch Motions It is assumed that the heave and pitch motions are measured. The gyroscopes can provide accurate measurements of angles, up to about 1/10 degree. The noise therefore is due to structural vibrations, which in the longitudinal direction can be significant due to the beam-like response of the vessel. As a result the measurement noise was estimated based on data from ship vibrations. The same applies to the heave measurement noise. A Kalman filter was designed for speed V = 21 ft/sec and waves coming at 00 (head seas) with significant wave height H = 10 ft. and modal frequency im = 0.73 rad/sec (sea state 5). The measure- ment noise intensity matrix was selected from ship vibration 70 data to be 0.75 0 0 0.0003 The model poles are shown in table 4.1, while the filter poles are within a radius of 1.3 rad/sec as seen in table 4.2. Typical simulation results are shown in Figures 4.1 and 4.2. In these figures exact knowledge is assumed for the significant wave height and modal frequency. The accuracy of the filter is very good both for heave and pitch. Subsequently, the same filter was used combined with a ship and sea model different than the nominal one, to investigate the sensitivity to the following parameters: The influence of the significant wave height is very small when the modal frequency is accurately known. On the contrary, the influence of the modal frequency is quite critical; particularly for pitch (see Figures 4.6 and 4.7). The same conclusion is reached when a double peak sea spectrum is used [22],[23]. The effect of the forward speed and wave direction was found to be unimportant particularly for heave, while for small changes in wave angle (+ 150) the pitch prediction error was not affected significantly [22]. 71 Sway-Roll-Yaw Motions As in the case of heave and pitch, the measurement noise consists primarily of structural vibrations rather than instrument noise. For roll such vibrations are quite small for a destroyer vessel and similarly for sway and yaw the vibrations are smaller than in the case of heave and pitch. The noise intensity used was nonetheless similar to the heave, pitch noise, so as to bound the filter eigenvalues below 2 rad/sec, which is the typical wave bandwidth. A specific example has been worked out for a forward ship speed of 15.5 ft/sec and waves at 450 and sea state 5 (significant wave height of 10 ft. and modal freqeuncy of 0.72 rad/sec). The measurement noise intensity matrix was diag {0.l ft2 , 2*10-4 (rad) 2 , 2*10-4 (rad )2} The simulation shows very good estimation as seen in figures 4.4, and 4.5. Yaw is very small and the measurement noise 4.3, is large relative to the yaw motion, nonetheless the yaw estimation based primarily on the roll, sway measurements is very good. Table 4.2 presents the results of a sensitivity study of the influence of the various parameters involved. parameter is again the modal frequency. The most critical The ship speed and the wave direction are not critical for the estimation error. This is a very important conclusion as far as the wave direction is concerned, because in reality, seas are directional and very difficult to measure, or even model appropriately. 72 The influence of systematic measurement errors was studied This factor is defined to be the by using a calibration factor. ratio of the measurement fed to the filter over the actual measurement, thus introducing a systematic error. If C is the calibration factor,then the systematic error as a percentage of the actual measurement is 100 (l-C). In the case of a 10% error, the most significant change was found in the case of the roll motion. In the case of a calibration factor 0 (indicating a dis- connected measurement) significant errors resulted, especially for roll in the case of disconnected roll measurements (Table 4.3). Table 4.4 presents the poles of the model used, while Table 4.5 shows the poles of the Kalman Filter derived for the nominal condition as described above. Figures 4.8 through 4.10 are sim- ulation results and show the significant effect of the modal frequency on the estimation error. Finally, Figure 4.11 shows the simulation of the sway motion estimation when the roll measurement is disconnected. 73 Ship Speed: U = 21 ft/sec SEA: H=10 ft, w m = 0.72 ra a sec Heading angle: 00 -d H=10 ft, w m = P1 1,2 = -0.199 + 1.111 same P = -0.286 + 1.016 i same 3,4 -0.576 + 0.576 i same P8,9 = -0.576 + 0.576 i 8,9 same P = 10,11 = -0.863 + 0.863 P10,11 10,11 -0.571 + 0.571i 2 ,1 3 = -0.571 + 0.571t P12,13 = -0.863 + 0.863 P1 P14,15 = -0.863 + 0.863 P 1 4 ,1 5 TABLE 4.1: sec P 5 = -0.696 P 5 = -1.058 P6, 6,7 0.52 rad = -0.571 + 0.571i Poles of the heave, pitch model 74 Ship Speed: U = 21 ft/sec, SEA: H = 10 Heading angle: 0° W m = 0.72 rad/sec ft P1,2 = -1.289 + 0.540i P34 = -1.134 + 1.033i 3,4 P 5 = -1.340 P 67 = -0.903 + 0.7280i 6,7 P 8,9 P1, = -0.777 + 1.376 i = -0.273 + 1.572 i P 2,13 = -0.248 + 0.9181i P12,13 P14,5 = -0.0936 + 0.0940i P14,15oles of the heave, pitch Kalman Filter TABLE 4.2 Poles of the heave, pitch Kalman Filter -& TABLE 4.3 Sensitivity of the RtMS error of sway, roll, yaw motion, to changes in the parameters of the ship-sea model. C (sway) indicates a calibration coefficient in the sway measurement (and similarly for the other motions), to handle systematic errors in the measurements. Parameter Changed Error Sway (ft) Error Roll (deg Error Yaw (deg) Basic Case 0.241 0.56 0.0776 U=20 ft/sec 0.245 0.568 0.0963 wm=0.52 rad/sec 0.314 0.91 0.0858 4=60 ° 0.296 0.624 0.112 C(sway)=0.9 0.255 0.586 0.081 C(roll)=0.9 0.247 0.708 0.0808 C(yaw)=0.9 0.242 0.56 0.0777 C(sway)=0.0 0.518 1.21 0.1408 C(roll)=0.0 0.376 4.08 0.158 C(yaw)=0.0 0.242 0.563 0.0785 RMS values of the motions (nominal case) 0.60 4.56 0.227 (.810)2 (.810)2 Measurement noise intensity (.316 ft:)z 76 Ship speed : SEA H=102ft 0 Ship model poles p = -0.754 3 4 = -0.754 : =0.72rad/s --A I=2t : + 0.754 0.754 .i P56 = -0.754 ± 0.754 P7 , 8 = -0.223 + 0.873 i P -0.335 + 0.588 i P11,12 = -0.260 11,12 + ' 0.440 i p 1 3 1 4 = -0.00983 + 13,14 0.484 L p1516 = -0.0204 T 15,16 0.05971 TABLE 4.4 _o C : an m=O.=0.-8G7rad/s : 1,2 3,4 e inc e. ?ea U=15.5ft/s + P1 2 '3,4 P5 6 = -0.470 + 0.470 i = -0.470 _ 0.470 i = -0.470 + 0.470 SA'ME : Poles of the sway, roll, yaw model £ 77 . 72rad,/s =t Kal.-.an filte-r ooles : ?1,2 = -1.067 ± 1.086 P3,4 = -0.457 ± 1.312 i -5,6 P5,6 == -1.279 ± 0.477 p7,8 = -0. 210 P9,1 0 = -0.365 i 0.934 i _ 0.523i ?p1,12 = -0.087 11,12 + 0.446 V = -0.159 ± 0.165 { P13,14 15,16 = -0.0203 + TABLE 4.5: 0.0595 i Poles of the sway, roll, yaw Kalman filter 78 7.50 4 5.00 I 2.50 Li 0 .00 6 -7.50 -52~. U1 CD. * m + m + , * m + + I * mlrr + * * m m + + id m + TIME (sec) Figure 4.1 Results of Heave Simulation and Its IKalam Filter Estimate (dotted line), Using Accurate Model at U=21 ft/sec and 0° Angle of Incidence and in sea state 5. r + 4.~ 7~~ 3.008 2.0880 1.800 1.00 - - -4.t''I FI 1 ' FI ffi rrl m TIME F m (see) Figure 4.2 Results of Pitch Simulation and Its Kalman Filter Estimate. Same conditions as in 4.1. rq rrl 80 3.00 r ~~actual F~~~~~~~ A estimated - I .001.84 - 8n.8 n_ *20 * * * mr m ~ m + + + s * -+ m m m + TIME (sec) Figure 4.3 Results of Sway Simulation and Its Kalman Filter Estimate (dotted line), Using Accurate Model at U=15.5 ft/sec and 450 Angle of Incidence, and in sea state 5. 81 I5.00 r actual estimated I0.00 5.00 L -) J -5. 00F -1A.00 C m m + -*- * - *+ CD -+ TIME + I w + (sec) Figure 4.4 Results of Roll Simulation and Its Kalman Filter Estimate. Same conditions as in 4.3. + 82 1.00 , actual estimated i! i i 25 - .00 ' 1 '1 a en 1 I I CDn-i X X 3 t3 ml + m X m + m + TIME m + X X 6U m + m + (see) Figure 4.5 Results of Yaw Simulation and Its Kalman Filter Estimate. Same condition as in 4.3. 83 7.50 5.00 i- 7 i58 11 A9 _.0 8k/\ 7LiaV -7.50 61 0) + + + CD CD _ - + + + + TIME + (sec) Figure 4.6 Results of Heave Simulation and Its Kalman Filter Estimate (dotted line). The actual wave spectrum model frequency is 0.52 rad/sec, while the value used in the Kalman Filter is 0.72 rad/sec. All other parameters as in 4.1. + 84 4.00 r 3.00 -2.0088 - -3.00 V I -4.000 . 0 . + · '. X + X0 9 . . · - m · X X + + + + + + + TIME (sec) Figure 4.7 Results of Pitch Simulation and Its Kalman Filter Estimate. Same condition as in 4.6. 85 'I. cn AX -3.00 rI -i .800 -8.F8 Results Sway Actual -4(dotted Simulation line). and rad/sec, Its W=.52 Filter while Estimate in Kalman Filter W= 7of + + Kalman + + TIME (see ) Figure 4.8 Results of Sway Simulation and Its Kalman Filter Estimate (dotted line). Wm=0.72 rad/sec. Actual Wm=0.52 rad/sec, while in Kalman Filter All other parameters as in 4.3. 86 2830.00 0 =10.00 · . . Same conditions as in 4.8. - - * U m r m Results of Roll Simulation Figure and Its 4.9 Kalman Filter Estimate. Same conditions as in 4.8. I 87 1.50 1.00 1.580 .// 1 yr -. 50 -I 0 .0 -I .50 Same conditions as in 4.8. - 88 2.00 i- n s. i 0.00 -1-.00 - i ! 2 IZ * (1 \1 ji s ~~~~-4.~00~ TIME (sec) Figure 4.11 Results of Yaw Simulation and Its Kalman Filter Estimate (dotted line), using noisy measurements (light line) when the roll measurement is disconnected. in 4.3. Same other conditions as 89 3.00 F 2.00P -2 _1.00 'iA en m'~ X :~'i ( uiv * Cn . Jl ! LIII m m r iiX ui se)i . .. mn m m Figure 4.12 Results of Sway Simulation and Its Kalman Filter Estimate (dotted line), using measurements (light line) when the actual angle of incedence is 600 and the value used in the Kalman Filter is 450. Same other conditions as in 4.3. 90 Chapter 5: SHIP MOTION PREDICTION It is of interest to use the models developed in the previous sections to forecast the behavior of the vessel within a few seconds. The feasibility to predict the motions could assist significantly the pilot in committing the aircraft to landing under favorable conditions. An automatic landing does not require within LOG theory such information since the predictable part of the motions is included in the state and therefore used directly. Nonethless the prediction is of primary importance for pilot landing or semi-automatic landing. Similarly for offshore operations a display of a prediction of the most critical vessel motions could reduce the operation risk significantly. The operator could choose a time window of minimal motion or acceleration and then transfer cargo or personnel. The subject has been considered in the literature [6], [11], [24] using both frequency and time domain techniques. Theoretical Background The first to treat the subject of developing a predicting filter was Wiener [29]. If a random process has power spectrum S(w) a spectral factorization is first required, i.e. A(w) = ( (C) 4* (w) where * denotes complex conjugation and f(s) (73) is an analytic function of s with the exception of a finite number of poles in the left half plane. Then the transfer function of the optimal predictor K(w), in the sense of minimizing the expected value of the error, is 91 given by the expression: K() Y(t+a)e - i = t dt where Y(t) is the inverse Fourier transform of prediction time. (74) (w~)and a is the The importance of this result is to provide a number of intuitive results, such as the fact that a narrow band process is predictable while at the extreme a wide band process is unpredictable. The disadvantage of this approach is that it may require differentiators in its implementation, depending on the form of the spectrum [29]. The alternative is to use state space models where no such problems appear. form. In fact the predictive filter has a very simple If the system has a state space description 0o x = Ax + W (75ab) Y = Cx + W where W 1 and W 2 are white noise signals, it is not hard to see that if the state x is perfectly known at t then the predictable part of x(t+T), denoted by xp, is AT x (t+T) = e -P x(t) (76) If the state is not available the Kalman filter estimate is used instead [5],[8]: AT x (t+T) = e -p ^ x (t) (77) 92 i.e. the propagation of the equation 0 -p(t+T) = A x (t+T) (78ab) x (t) -p x (t) For a stable matrix A x (t+t) - 0 as T +-D, reflecting the fact -P that the influence of the driving white noise completely alters the state of the system, once the homogeneous solution has died out. The covariance of the error e x -p - x denoted as Pp is governed by the equation T 0 Pp () = A Pp(T) + Pp(T) A + V1 (79) When the state is perfectly known the initial condition is: Pp (o) = 0 (80a) While in the case of using the Kalman filter estimate P (o) = P(t) (80b) where P(t) is the error covariance of the Kalman filter at the "present" time t. The two models developed for the vertical and horizonatl motions have been used to study the predictability of the ship motions. Figures 5.1 through 5.5 show simulation results assuming perfect state knowledge. Similarly perfect state knowledge has been assumed and the covariance has been propagated using equations (80). (79), Figure 5.6 is a plot of the heave, pitch motions rms error versus prediction time, while Figure 5.7 depicts the sway, roll, 93 yaw rms erros versus prediction time. The error has been non- dimensionalized with respect to the corresponding rms motions. As expected, the error tends to 100% for large prediction times. Roll is a narrow band process and as expected it is the most predictable motion, up to ten seconds ahead. The remaining four motions are predictable up to five seconds ahead. These results hold in the ideal case. The actual performance will be lower due to the presence of noise in the measurements, fewer measurements than states and modeling errors. To assess the effect of measurement noise some simulations were made, whose results are shown in Figures 5.8 through 5.10. extreme case was considered: A rather In the case of the vertical motions only two measurements were available, heave and pitch, and in the case of the horizontal motions, only the sway, roll, and yaw motions were available. As seen in the figures, the same noise used to derive the Kalman filter gains was used, which is quite significant. As expected, the performance deteriorated although roll is still predictable up to eight seconds. The other four motions are pre- dictable up to about two seconds. In the case of modeling errors, let the correct model be 0o x = Ax + W (81) while the prediction model is o x with A* = A + 6A. = A* x (82) Then the error e = xp - x is governed by the equation e = A*e + 6A x - W (83) 94 or by letting xl = o {eT FA* A 0 A X1 F-WI x + W (84) so by denoting: P = E{e T , = E{ x e ,V = E{ x xT } (85) the following equations are obtained P = A o V = A V P +* P +A*TA-V + T-AT V T T V A * T + U 6A O U = A.U + U A + V1 V1 T + V! UAU~~t ~(86abc) +V~ 1 U, the covariance of the vessel model, can be assumed to be in steady state so the first two equations can be used to propagate V and P. Figures 5.12, 5.13, and 5.14 depict the error covariance of the vessel motions when the model used is different than the actual one. Since one of the most critical parameters is the modal frequency, its effect has been studied: The nominal value is 0.52 rad/sec, while the value used in the prediction filter is 0.72 rad/ sec. The covariance at the initial time is assumed to be zero (perfect state knowledge). As can be seen from these figures, the effect of modeling error is important in the case of the modal frequency, providing a reduction of about 30-50% in the prediction time within prescribed confidence limits. 95 The sinusoidal behavior of the covariance propagation at about twice the motion natural frequency roll, (as seen in the case of for example) can be explained by the form of the covariance equation P = A P + P AT - V (87) For example the unforced equation 0 T P = A P + P A (88ab) P(O) = Q has the solution At A't P(t) = e Qe (89) which is composed of exponentials in the form (hi + Xj)t where Xi are the eigenvalues of A. As a result 2Xi will appear and in the case of roll it is obvious that twice the roll natural frequency dominates the response. 96 Reduced Number of States Another aspect of interest is omitting states. In such a case we denote by x the nominal state of dimension n and by x* the implemented state of dimension m, n>m. We assume that x* is obtained by simply omitting some states, so that x* = Wx (90) where W = [I m 0] with Im the unit (mxm) matrix. (91) The nominal system equation is 0 x = A x + BW (92ab) z = Cx + W2 while the prediction filter is 0 x =A* x (93ab) zp = C* xp where x has dimension m, A* is the mxm reduced system matrix and -p C* the reduced observation matrix. Then we define: A = W T A* W -A (94) and we obtain the covariance equations [8]: P = T-A*-W-P + P-WT.A*T-W + 6A-V + VT 1 V = A-V + v.wT A*T.W + U-.AT - V 1 U = AU + UAT + V A + V (95abc) 97 where P, V, and U are the same quantities defined before, except that the error is defined as: e = WT x -p - x (96) Note that we may eliminate any row from the original system, so its is convenient to set all the rows to be eliminated at the end by performing a row permutation. The matrix that interchanges rows has the form S 0O ... O 1 0 ... 0 (97) where each row has only one nonzero entry equal to one, and no row is the same as any other one. The inverse of S is ST thus minimizing the computational effort. Once the appropriate permutations have been established, it is easy to construct S and equation x = Ax + BW (98ab) y = Cx becomes z= A z + BpW (99ab) y =Cp z where = Sx A = S A ST B = S B C =C p (100abcd) ST 98 Then we proceed to determine the error covariance as explained above using Ap, Bp, and Cp as the system matrices. The inclusion of non-minimum phase zeros was considered to be an important part of the overall modeling. This was confirmed by studying the effect of omitting these zeros on the prediction error covariance. This is seen in Figure 5.15, where the heave and pitch rms error, non-dimensionalized over the corresponsing rms motion, is plotted versus prediction time. As expected, pitch error increases substantially, since pitch lags heave at low frequencies by 90° . Because heave and pitch are coupled, the error in heave is also affected, resulting in poor prediction of both heave and pitch. 99 3.880 2.880 1 1.80 a I -2.0088 -3. 88 W a X^W 1 aX Ln n A Lat X 8 63 CA 01 W s 0a n a TIME (sec) Figure 5.1 Heave Simulation Results and Its Prediction (dotted line starting at t=40 sec) for U=21 ft/sec and ~=0O, and in sea state 5. Perfect state knowledge is assumed. 100 2.88 1.58 -2.00 .50 . -28 8 · CA LA CR 'Ca TIME (sec) Figure 5.2 Pitch Simulation Results and Its Prediction. as in 5.1. Same conditions 101 .75 ' i I J ii .50i ",1 1X t / ,', \ I ir PerfectXcstate knowledge~ is assumed. c r / , . i ~ us TIME (e:Xn ) Sway Simulation and and Its Prediction linestate starting at t=40 se ) for Results U=15.5 ft/sec =45t, an (dotted d in sea 5. Perfect state knowledge is assumed. 102 7.50 ;f' 250 U I ' a. 00 U/\ X |X! i -2.50 ' /\/ I / -5.00 -7 .50 Lo. _ -W 7.5 X tW a ~Ln (,0 - ± L3n1n ,,X an 0X Ln 5} .. TIME (sec) Figure 5.4 Roll Simulation Results and Its Prediction. as in 5.3. Same conditions Ln "(. 103 .30 .20 A .\ / I I (A i i! ha -·. L in 5 [a h) A A . . . u . u O\i .. TIME (sec) Figure 5.5 Yaw Simulation Results and Its Prediction. as in 5.3. Same conditions ... 104 PREDICTION ERROR 12.SeE+01 Fs Z 18.88E01 W 0 Z O 58.88E+ n:/ - HEAVE , Ix 0 Ixy- o 25.8E++8/ PITCH - _ a a6 (it _ 6a (t 0 X TIME (SEC) Figure 5.6 RMS Prediction Error Over RMS Motion Versus Prediction Time for Heave and Pitch, U=21 ft/sec, Wm=0.72 rad/sec, sea state 5. m 105 PREDICTION ERROR 12.50E+01 - I- Z 10.0E+801 W Z 75OOE+00 / 0 H2 aO 25.0OE+/_ /,rS/ CA l OL ROLL z L Cc SWAY TIME SEC) Figure 5.7 RMS Prediction Error Over RMS Motion Versus Prediction Time for Sway, Roll, 00.00£-01 Yaw. U=15.1 ft/sec, =450, sea state 5. a · ~~----------------~~~~~~ - ~ ~ __ 106 HEAVE SIMULATION, PREDICTION 7.50 5.00 I 2.50 -5.00 1| ACTUAL 1MO0TION i - -8 j1 -7.50 -10.00 NO I SY MEASUREMENT PREDICTION . p B, Bt TsM U1 (s e Bt ) TIME (sec) Figure 5.8 Heave Simulation Results, Its Kalman Filter Estimate (up to 40 sec) and Its Prediction Using the Kalman Filter Estimate (after t=40 sec). Same conditions as in 5.1. 107 SWAY SIMULATION, PREDICTION 3.00 2.00 1 1.0088 -2.00 - 1 IEASURED - C- -3.00 . Ca I I I (.a C) CA CA PREDICTED i CAl CA c CA TIME (sec) Figure 5.9 Sway Simulation Results, Its Kalman Filter Estimate (up to 40 sec) and Its Prediction Using the Kalman Filter Estimate (after t=40 sec). Sam conditions as in 5.3. 108 ROLL SIMULATION, PREDICTION 2.50 88 -7.5 PREDICTED 5.881 -18.00 jFigure , TIME (see) 5.10 , I- ---- Roll Simulation and8.IPrediction. Same conditions as in5.9. -5. ~88 1 1i ACTUAL ___ Roll Simulation and Prediction. Same conditions as in 5.9. ' 109 .75 actual .- - estimated/predicted 50 .25 -.50 ll -. 75 L *·¢ Ln -a CD X ·9 Ln 0s 06 TIME Li m) Ln m 0 Ln · ·s (sec) Figure 5.11 Yaw Simulation and Prediction. Same conditions as 5.9. CD 110 .75 actual .... - estimated/predicted .508 .25 I' I I 0 00 P"3 i/ - - \ -. 25 -. 50 -. 75 I T ·1 a) TIME i T 0 n z , , CD 0 (J (sec) Figure 5.11: Yaw Simulation and Prediction. Same conditions as 5.9. 0(A PREDICTION ERROR (SWAY) 15. E+81 12.58E+81 Uti d 18.8+81 z. / = O- 0 0 75.88E+8I a88. BEs- TIME(SEC) Figure 5.12 RMS Prediction Error Versus RMS Motion Versus Prediction Time for Sway. Actual Wm=0.5 2 rad/sec, used Wm=0.7 2 other conditions as in 5.3. rad/sec. All 112 PREDICTION ERROR (ROLL) 8. 88 z /c 7 8.8 - 0.08 48.88 ~Figure 18.88 8.88 4i en 1\/ 5. 1388 ea 6sw TIME (soc) Figure 5.13 RMS Error Versus Prediction Time Roll. in 5.12. Same condition as 113 PREDICTION ERROR (YyW) 12.58+81 ItL I- I II I / CDn t a: =2S.EWFtOot . TIME ( sec ) Figure 5.14 RMS Error Versus Prediction Time for Yaw. as in Ero P4 essPeicinTm Same conditions 5.12. o aw aecniin 114 PITCH HEAVE Figure 5.15 RMS Prediction Error Over RMS Motion for Heave and Pitch. U=21 ft/sec, 4=0 ° , sea state 5. In the prediction model the nonminimum phase zero have been omitted. 115 CONCLUSIONS A satisfactory approximation of the ship motion equations as provided by hydrodynamic theory has been achieved. The approximation is valid within the wave frequency range and for seas described by the Bretschneider spectrum, whose major limitations are (a) uni-directional seas (b) unlimited fetch, deep water The resulting two groups of motions, i.e., heave-pitch and sway-roll-yaw can be approximated separately requiring 15 and 16 states respectively. If both must be used a 25 state system is required. The model depends parametrically on the ship speed, the wave angle and the significant wave height and modal frequency. The Kalman filter is designed using as measurement noise intensity, values from ship vibration amplitudes. For sway- roll-yaw the vibration levels are small, nonetheless, to bound the filter eigenvalues below 2.0 rad/sec similar values were used, as for heave-pitch. It should be remembered that heave and pitch are related by a non-minimum phase transfer function, resulting in reduced filter accuracy. Actually, heave is 90 ° out of phase for low frequencies with respect to all the other motions. 116 A sensitivity analysis of the filter performance indicates that the most critical parameter is the spectrum modal frequency. It should be remembered that a sea spectrum may contain more than one peak, in which case it is essential to obtain an accurate estimate of both peak frequenices. Of particular interest is the fact that the wave direction does not have a significant influence on the estimation error. This means that although our modeling used a uni-directional Bretschneider spectrum, it can be applied in its present form for directional seas. The models derived herein can be used to predict the ship's motions up to 5 seconds ahead in time for all motions and 10 seconds for roll. When modeling errors and noise are taken into account, a more realistic estimate of 2-3 seconds for all motions and 6-8 seconds for roll is obtained. Again, the modal frequency of the sea spectrum is the most critical parameter. Also, the nonminimum phase zeros can deterio- rate the performance of the predictor significantly if omitted. 117 ACKNOWLEDGEMENTS Erik Tiemroth, graduate student in the Ocean Engineering Department, has prepared some of the programs presented in this study and assisted in the debugging and running of several others. Material and information was provided by Prof. Gunter Stein of the EECS Department and Chris McMuldroch, a graduate student in the EECS Department. The work was carried out at the Laboratory for Information and Decision Systems with support from the National Aeronautics and Space Administration Ames Research Center under Grant NGL - 22 - 009 - 124. 118 REFERENCES (1) Athans, M., "The Role and Use of the Stochastic LinearQuadratic - Gaussian Problem in Control System Design", IEEE Trans. on Automatic Control, Vol. AC-16, Dec. 1971. (2) Balchen, J.G., Jenssen, N.A., Saelid, S., "Dynamic Positioning Using Kalman Filtering and Optimal Control Theory," Automation in Offshore Oil Field Operations, 1976, pp. 183-188. (3) Barr, R.A., Ankudinov, V., "Ship Rolling, Its Prediction and Reduction Using Roll Stabilization," Marine Technology, Jan. 1977. (4) Bodson, M., "Lateral Control System Design for VTOL Landing on a DD963 in High Sea States," S.M. Thesis, Dept. of Electrical Engineering and Computer Science, MIT, Cambridge, Mass. 1982. (5) Bryson, A., Ho, Y.C., "Applied Optimal Control," Publishing, Co., New York, 1975. (6) Dalzell, J.F., "A Note on Short-Time Prediction of Ship Motions," Journal of Ship Research, Vol. 9, No. 2, 1965. (7) Doyle, J.C., Stein, G., "Robustness with Observers," IEEE Trans. on Automatic Control, Vol. AC-24, August 1979. (8) Gelb, A., editor, "Applied Optimal Estimation," MIT Press, Cambridge, MA. 1974. (9) Grimble, M.J., Patton, R.J., Wise, D.A., "The Design of a Dynamic Ship Positioning Control System Using Extended Kalman Filtering Technique," Oceans '79, San Diego, CA, Sept. 1979. Hemisphere (10) Kallstrom, C.G., "Simulation of Ship Steering," Dept. of Automatic Control, Lund Institute of Technology, Lund, Sweden, Coden: LUTFD2/(TFRT-7109)/1-353/1976. (11) Kaplan, P., Kaff, A.I., "Evaluation and Verification of Computer Calculations of Wave-Induced Ship Structural Loads," Ship Structure Committee, Report No. SSC-229, 1972. (12) Kaplan, P., "A Study of Prediction Techniques for Aircraft Carrier Motion at Sea," AIAA Sixth Aerospace Science Meeting, Paper No. 68-123, 1968. (13) Kim, C.H., Chou, F.S., Tien, D., "Motions and Hydrodynamic Loads of a Ship Advancing in Oblique Waves," Trans. SNAME, Vol. 88, 1980. 119 (14) Kwakernaak, H., Sivan, R., Systems," Wiley, 1972. "Linear Optimal Control (15) McMuldroch, C.G., "VTOL Controls for Shipboard Landing," Laboratory for Information and Decision Systems, Report LIDS-TH-928, 1979, MIT, Cambridge, MA. (16) Newman, J.N., "Marine Hydrodynamics," Cambridge, MA. (17) Ogilvie, T.F., "Recent Progress Toward the Understanding and Prediction of Ship Motions," Fifth Symposium on Naval Hydrodynamics, Bergen, Norway, 1964. (18) Salvesen, N., Tuck, E.O., Faltinsen, O., "Ship Motions and Sea Loads," Trans. SNAME, Vol. 78, 1970. (19) Sidar, M., Doolin, B.F., "On the Feasibility of Real Time Prediction of Aircraft Carrier Motion at Sea," NASA Techn. Memo X-6245, 1975. (20) Smith, O.J.M., New York 1958. (21) St. Denis, M., Pierson, W.J., "On the Motions of Ships in Confused Seas," Trans. SNAME, 1953. (22) Tick, L.J., Differential Equations with Frequency Dependent Coefficients," J. of Ship Research, Vol. 3, No. 2., 1959. (23) Triantafyllou, M., Athans, M., "Real Time Estimation of the Motions of a Destroyer Using Kalman Filtering Techniques," Laboratory for Information and Decision Systems Report, MIT, Cambridge, MA 1981. (24) Triantafyllou, M., Athans, M., "Real Time Estimation of the Heaving and Pitching Motions of a Ship Using a Kalman Filter," Proc. Oceans '81, Boston, MA, Sept. 1981. (25) Triantafyllou, M., Bodson, M., "Real Time Prediction of Marine Vessel Motion Using Kalman Filtering Techniques," Proc. OTC, Houston, Texas, 1982. (26) Vugt, J.H., "The Hydrodynamic Forces and Ship Motions in Oblique Waves," Netherlands,Ship Research Center TNO Report No. 150S, Dec., 1971. MIT Press, 1978, "Feedback Control Systems," McGraw Hill, 120 (27) Yumori, I.R., "Real Time Prediction of Ship Response to Ocean Waves,Using Time Series Analysis," Oceans '81, Boston, MA., Sept. 1981. (28) "5-Degree of Freedom Seakeeping Program Manual," Design Laboratory, Ocean Engineering Dept., MIT, Cambridge, MA, 1979. (29) "Dynamic Control Systems Software - User's Manual," Laboratory for Information and Decision Systems, MIT, Cambridge, MA. 1980. (30) Wiener, N., "Time Series," MIT Press, Cambridge, MA. 1949. 121 APPENDIX 1 Hydrodynamic Theory In the text, the simple one dimensional equation of heave motion was derived to demonstrate the principles involved. Here we will proceed to write the overall equations of motion. We will avoid extensive hydrodynamic theory developments since [15 [17 ] provide an in depth coverage. ] and Within linear theory, we intend to write the added mass, damping and exciting force terms. The equations of motion can be written as {_-2 [M + A] + ix = B + C} x Fx n (1) where M = {Mi} B = {Bij} , C = {Cij} A = {Aij, (2) We can find the various matrices from hydrodynamic theory [15 ], [ 17] and dynamics. We omit surge as a second order quantity so that x is a vector of dimension five: m 0 O m 0 0 I 0 0 -mz c O -mz c xx 0 O ° O O 0 xz I YY 0 0 with m the mass of the ship, I z O' (3) O I the distance of the center of C 122 gravity, vertically, from the origin. A2 A = 2 0 A2 4 0 0 A3 A2 6 O A i A 42 0 A44 0 A4 0 A 53 0 As5 0 A62 0 A 64 0 A60 B2 0 B2 0 B2 6 B 33 0 B35 0 O B4 0 B4 6 . O Bs 3 0 Bss 0 B 62 0 B 64 0 B6 0 0 0 0 0 O OpgA O O O 0 0 A(GM) 0 O C5 3 0 0 0 0 33 5 0 6 (4) fI 2 0 B C BB 4 = 2 4 4 (5) 6 O0 (6) A(GM)L 0 0 0 With Ap the waterplane area of the ship, A the displacement, (G'M) the transverse metacentric height and (GM)L the longitudinal metacentric height. The following relations hold true within strip theory [17] 123 A B 33 5 3 B = B° + U a33 A=3 s A° AO B - + 2U UA Bo + UA3 3 - 33 35 3 5 35 As 3 U xbA b3 3 Ux 2 A W Ux a A 33 0 U B3 3 + U + B5 - UA 33 = 35 - UxA a Ass As3 = = A 5ss As + + U2 -22 B = Bs5 + U2o U B 3 3 + UX The superscript = ~ A w 33 - A xA b 33 A o B5 3 U2 U2 A 7 b 33 A o *3 2 U 0 A3 33 - 0 denotes quantities 33 2 U 2XA U2 A 3 + at A ,2xA a 3 2 b 33 3 + U zero xA A 3 speed. XA is the distance to the aftermost cross-section of the ship and a.j, 0A.are its sectional added mass characteristics. These last two 13 quantities are important only for cruiser sternships. B?. can be found using the M.I.T. Seakeeping Program 13 Similarly, we can find: U A2 2 - 2 A 2 x2 B2 2 + U 2 0 A2 2 B 22 = X2 The Aj, [271. 124 A2 . = A4 2 = A2 B2 4 = B42 = B2 4 + A2 6 = A26 + B2 6 = B2 6 - A4 4 = A4 B44 = B44 b -- A O A 4 4 U 0 =U 2C2 U B22 - with B44 A U- + A U2 a22+ A C2 2 A b2 A U 4 + UX 2 + XA b 2 2 o0 UA 2 U2 A bA * + B44 A U a44 the equivalent nonlinear damping A4 6 U o = A4 6 + 2 B2 4 - U2 A A o o B62 = B26 A6 4 = A4 B64 = B46 UA 2 U 6 "- + U + UXA a2 , 2 U 0 2 U2 A xA b24 + B2 ° A XA A b2A4 -2 UA24 + Ux A a 2 4 -2 2 A A a2 125 A6 5 AG= = B66 P f F4 = a p P = 2 U U A22 - BB 6 6 + U2 B2 2 F 2 = ao F3 + o p a FS = -P (f 2 x X 2 b A 2+ + + UXA C22 h2 ) d + U U2 2 U iU + co. p iW 2x A A XA b A2 A hA tfP U- A+ ¾ [[ + hdE + ai-mh2 2) + (f4 (f +2 h4) p iU di hA o f U (f 3 + h 3) dS + o f [[(f 3 p ao U + h 3) + iW h a22 i U hAA 3 ] dE - P p s' -- U XA hA Xh hA where fj is the sectional Froude-Kryloff force and h. the sectional fraction force diffraction force. 126 VCG = Vertical Center of Gravity = - 4.6 ft. (from waterplane) GM = Metacentric Height XCG = Longitudinal Center of Gravity = 1.07 ft. AFT = 4.16 ft. ( from admidships) A = Displacement Cb = Block Coefficient b = 6,800 ton = 0.461 127 APPENDIX 2 The M.I.T. Sea-keeping Program [27] was used to derive the hydrodynamic data. The values used to develop the simple models for heave-pitch and sway-roll-yaw (as already mentioned, to first order heave and pitch are uncoupled from sway, roll and yaw). All units are consistent such that the forces are obtained in tons, the moments in ton-ft., the linear motions in ft. and the angular motions in radians. The programs change the angular motions and express them (output) stage. in degrees only in the final Heave-Pitch Characteristics = 214 C33 587 A°33 = 281 BO 260 A°35 = 15500 C 35 = 260 M 33 B I55 = 3.76*10 A B5A = 3.8*10 C55 5 5~C5 Al = A2 =120000 15500 35 55 = 4.20*106 9.53*10 550 The principal ship characteristics have as follows LBP = Length between perpendiculars B = Beam = 55 ft. T = Draft = 18 ft. = 529 ft. 128 Sway - Roll - A44 = I44 = 104,000 B 44 = C4 4 = A A46 A 22,800 800 + B44 28,800 -760 -760 24 = 181,000 220 22 Yaw Characteristics B24 = -50 B46 = 5,600 B 22 = 10 66 I6 = 3.8 * 106 B = 130,000 66 66 A26 26 = 14,500 A = 380 A, 4 A6 6 = = 2,400 23,000 23,000 B 26 =370 129 APPENDIX 3 Simulation The simulation of a continuous system x = Ax + BW 1 y = Cx + W2 where the white noise signals W 1, W 2 have intensity V 1, V 2 respectively, is performed by constructing the equivalent discrete system x(T+6t) = Ak x(t) y(t) + Bk Wkl = C x(t) + Wk 2 A2 where A k = e Aft =I + A-6t _BW 1 At+6teA(t-T)BW BkWk = £t+6t e BW k kl t 2/ l( +. T)dd An approximation would be x(t+6t) = (I+A-'t)x(t)+(B-6t)Wkl y(t) = Cx(t) + Wk 2 where Wkl, Wk 2 discrete white noise of intensity Vkl, Vk 2 respectively, given by: 130 Vkl = V 1 /6t Vk 2 = V 2 /6t If a random number generator is provided with a range between 0 and 1, the following relation provides Wkl (and simi- larly Wk2) Wkl = (RND-0.5) · 12-Vkl where RND is a random number between 0 and 1. Higher order approximations to Ak may be necessary in some cases. Problems appear in particular with lightly damped systems, i.e., when the matrix A possesses eigenvalues -a + ib where O<a<<b. Such for example, is typically around 10% is the case of roll whose damping of the critical value. The approximation: e (-a+ib)6t (l-a-6t) is valid provided a6 -t<<l b3 t<<l The next order term is: (-a+ib) 2 6t 2 so for accuracy we require, (a2 -b2)-6 I ab6t 2 <<b6t <<a6t + ib-6t 131 which results in a single additional condition 6t<<i 2a 1 a -b It is easy to see that when a<<b then 6t<< min [1 ' 2 (2)a This requirement may be very-demanding for simulation by imposing an extremely small time step. If an appropriate step is not chosen then the real part of the (neglected) second order term reduces the first order real part, thus resulting in a reduction of the already small damping. By using a higher order approximation, in the form e A.6t I + A-6t + such problems are resolved. . .+ Am 6 t m In the present case, the sway- roll-yaw model requires such treatment for efficient simulation. 132 Appendix 4 Computer Program Listing 133 Cd a, a, O O 0 4-- CL0)~~~~~U 0 o J O 0 ~ 4 Cd - cc~~~~~~~~~~~~~ 0 06O O~a 0 0O0 0OOGOOO O OO 0OOO0OOO o o O OC, o OO OOOOZOO OOO oooo.oooo U~~~~ t oVI- w krdoooooo.ooooooooooooooooooooooooooooo rl~ :4 tu o:o c o o o o o ~ CC U"I 0~~~ "~' -Z C3QO OQ 4 if~~~~~~y000000000000000000000000000000000000 cc ,4 U~~~~~0U U) 0n zf o ~ Q nn0~ o O o o'o o o o o o QSn~ o o o 00 o o o 00 ~~~~U) ", r O 00C)Q000000600000000000Q00000600000000ZDOOOOi3a 0 - Q, o 0OOO N .J w 4 0~~~~ w ·r UJI C., - -~: ~C'J -. * r , C~~~~~~C Z U. 0 4 . dr Z 4 -)( r. 03 CJ* U. (Cl -4 > 4I -. 3 - (_9 +<. X U-. T.-. , C - '- rC 0* a, . 0 ,U _.nU 0 -" W, ·- r '. 00 Z-O O~.0 0 0I- O~ -C')*--c 3- ?LU ~ .i~~ ~ ~ cr U- s-- Ir OAL 0 Li . -c vi-- 0 ' ~~ w Z, OC 0 U)~~~~~ O -~ -Z· 0* i- 0-4 *0 -u) (0-5 u.--Z LL ii I'3 ,X '--~ i (NO 0 d->0-Z-!n ICC -) -- (" 1 * (1D 0 I- 0 - --- C.·-- Z) 0 In-- * 0% .k =tfl *C4 I: ~ 0+-0: ' 002 0 O 0 i , ~ ) '- 't nL -'.* C4 ( U0~ l',--. ,.·" 4C40 m· I C) 00 *0.-00 ~ .- 0 C)r 0 in Cui (OKciCJ- ; - ~ 3L~ cdi 3~)' 010 0 ji ~ Cf~ )a n *~~~~~~~~~~~~~~~~W n M r ) 0)0ai-,c na~L77~ t'- * L In , .- , ono .) L*0 n .n 0 II C~ .'*")-k+ ·. 4.- + i- + + + >-*0C 00000000 (N 3 *oo *000,i, K-0u1K-(0C O I' 1:~(II U) (0IC1 C) co ~ 0· I hA~Jh c0 0 ·,I - l - 4c~0c U- * 4T7 . - t-U.<0I---0(.4"--0i *ri-c:: C n- -a 0 _ 0-+ r- M W U: 4 \ C - I <-. -· O-+ - -- 04 0 *-(aJ 0 4 -- 00 00J) C O f--3 W W 3+ *C1 0u)0-'. )00 ':I .~ !I ' - .0 IC 1 (Z t++*+.J 3 L )-.--I: . ~*-( T r-+++C1C· C· 0C1 K, K I C' r 134 w 0000000000000000000000000000000000000000000000000000000 LI) OOC'0000000000000000000000000000000000000000---------- OOOOO0000000000000000000000000000000000000000000000000 C000000000000000000000000000000000000000000000000000000 0 * * cO U) + ,o O 0 0 - . .J _ O Z+ t -L 0 L X O----*LZ 'tO- sL-t O ** I LLOsCb91* 5 L vi LL*+ 1 - r;i tn) J u n rI .-t 0II^-t-J 'i1 -Fii )C-o. i0 1^i-a Nt' AEf 1-iO*OO**O****SS7Z-1 CiV 11 l1 L 41 1--OCS. - F1l1 '. ," 7 it siC ti1 l1 -l 1>>i - ; ^sFr S :ro<C l* S _S - U ) _ : *r-_..n ..***C'Ci*O Z , 0 * Q DO It l ;I I II tln I u .,0 2U.'.-U. O+v;s -- 1) L 0 U* I II I *' In0 -W' It N -: L N Z 1 II C - ZN---) _ , _ _ 1 ^ -N u .. I II xt ;X - O I 0 i , 11 - 0 - 1i - Lf) - C- 0 i I I- 0 I C I +fIn-~fCA IN('3J0 O O h C C *! '--. O_ > . ', C'9 C . CO - - - - t- - - --"' It - - N C)j 11 i 135 0 0 .- C 000000000000 0000000000000 c¢ I0 z O 0 z ft: L. oo Cooooo *-* 0 - N --U. ~. r ' LL. O! LL o '? Lf)' C-, -00(' -, .' ,--CN Lu Z: . N 0~uOC 0'30 > <~% 136 O 0 Ul oH o 0 Cd 3 r1 o o a4 4 r4 0000000000000000000000000000000000000000000000000000000 000000000o o o 0888000000008 o o-----------o o s o o in cs tn cl v)X v) cs cs cl o Co w v c n v ) et )o c 0 w tn tn c 1 v0 v) v) 8800000000000000000000000000000000000000000 0000000000000000000000000000000000000000 in- rxCS ~~ z ct ox tA0( LAU Z ctx 0 1 0 In LU C. (0 - C* -- - C - x -(OD L aM x j >- x -x 1aL in I 3 V)Z CZ 0 >- 00 LUL hi) 1-CxO Lw Z x>Cx 0 C Q C in a >- uO U1 4 J LUc4I0>In 1 QO 2(x.j 00000 _ )-C -3-Q ZZZZZ 00000 - - -ac - 4 w -a 00000 0 Z inininn -0(ZC-t --- iLU -J'- wLULLUULw -j-j-j-J O mcOo O O a Z 00000 cOc00 * 0 0- 0 L lU 0-v in 000 * 0 NJ- * C U) Q--- 0- N 3UCx m in aL 0 0 x LU* *N + -N i V) 0 00 4 - .0 zo *Cu m w c-+ 1tC N C) *1. .'.r -a* .'OinO * ' ZZ > IIIUco-U)MO) *It <<4-< .00 4 II O III Z : _jwwxwwOU.0 4NQOmU)n>> N o - -C4 II - *- - 0 h - .- t -C-I---ii CZ t 0. i ni Z - 0* * Ol *C . * C'(CD CD U 4 ZI 1- t it 1 --- n-C4 -C, t II it if ' 00.i 0C'4C ,) U-O 0 + *aa0iOU, C) Lo OCx WQ+ *I41NNN000 D S) 1 - 1*- O--C 0 LI Q *xhIUIII ZV* 00'-' 0II U a 0 La- I 0 a in n- 00000 11 0 > a rCNI t * 0 00 * 4 *>> 0 a In ao 0 0 I-in**in*>++ in * * -<1-'in0 ZZ'-Z000 -ZZ - 4 *IIa _ - . c * a.'' * * -*J Cr-.' LL _ ZZ + * C'.JIO J+ ** _-- *4 Z* 0 ci * ejnu ,4* 1I 0.LUM Z xO ->cw -jN -X J *)i <-Cx0 -( I 0 CII * -LflD-0 ;-y 1>'IIL - -LU4 ' LU w I--coI W O -aQ C it-CxZ 00 1- 0 LU 3xow0cx04 Cx - -- i L*, IN C4. . *. 10LUu.* ZC 0- Q* C'4 . a C :3 U. i.j . t4 _ ZO' I Qn n CD7 *C'J * N e 3 L i- 4 LUV I- -I L - c-IL N 0* I ,a)nL Z N L. > 0 CL -- LULLULU w -t J _)l*N xx Cxa1 0..C 0 r-* Z 4 - ('.C')3 _0 * 1 * IL* U. - II 0 ' -a - Lu *3 3 . . vC _ 'CD _ 0 Z.O '-i ,a 0 Z 4 (UJ 0 - 0 4 0 CD *-. _-(D CD 0 05 4 LU C) I ** LU C _ - . 000 * tc 1 .. <<)M If 0 137 0 0 2 w 0000000000000ooooooooooooo000000008oo000ooooo00000000000000000000ooooooo LZ U: 0 DI0 f l Q ) r m m 0 0 o o o~ 0000000000000000000000000800000000000000000o;0 000000000000000000000000000 00000000000000000000000000000 c0 0 zw, _ ~~0I 0 0 W1) w 2.O- OOWW0 -- * * 4 U. L 5 O0 . > , mm <t In z * ID C<4CIJC.- - N N N- - U t co o * 0 1 0^ Z CN N Q CK Q C .w1 Q 4C-_ C'-- - - * *-II CY. -- . - IN - N - *if0 11 i0f N-N- - - - - - - I r _____ - - - 0 1 -N- - - - - - - - - - i N * '- - - - _ - _ N _*- _ 0 N znwwo - + . U. N i1N Q * 0 N >. ZZtJ O O O O-0000wi+v*0vJ - NZ N Q C C - * 3 I- It'J:)-- Q - ^ -I IN Z- 4 - ^ - - - 111 138 Cl 0 0 00000000000000000000000000000000000000000000oooooooooooooooooooooo0000000000 w 0000000000000000000000000000000000000000000000000000000 V) z t z LL o o C.) D U) _ w 4 Cr -'~-QC-CM --(N 0 U. O 0C CSI - 0 I- -^^f( N- v *T tn- Wle - > -(NC.4C.4(N 0Ur:~ U-~ ~ Z : *0 O o It 11 1 --I 1t -C D o t ** o 00N 0 -J * .. U. 't* s *C-0 (N*U C-) -- (0Q0 - 1 1 1 11 w O0 . 0 ~d~' ~ ~~,c 0X< C.4 0 ..'- * *t - aq- v T 1.D <Cc ~ 0 * x 0. . 111 ~ O< -- -- - o.0 W M r- ( '~ CM ** . 4:I - C IS O* 5L * O-, -W _ c 1i It 11 1 - - - - -W * 0 * *1 OCN * -C CI 1 111 1 if 11I 11 - -CIJ .0.'- * l( cr I( ^ I + *C =C.) tt r- Z 0 + - o + -0 *C- a.** 4Z*a. *v * CI *a.* 00 0 t 1 rI I I hW~dlI hIAl~ tol (O v , C'SC0tt O0. .* -* *C I Ir - - - - M ^ O * (Dw C it - 'u - D -o-o Ir * *-"cs --o o o0 0O I) Ih -0 IT 139 0 0 w CD w _ -- -0 0 0000000000000000000000000 -0 - - -0 0 UX _ Z o _ nz~~~~~~~~~~~~~~~~~~~~z ct t-t n I z X n, Z o XI o z w ,I0~~C I , I, w- II I I I I _ _ 3 - 0. a . 1 t , <U - c t Q . - -^ Z ^vC~vQ + . _ _ _ _c D D. 0* 0. -S- - ._ . tn _ _ _ _ 0. * -vvCv - A - ·- · - -. vva ^^- - (flN -CI 0* . *_ tw h- ,_, _ -ma _ _ _, -4 - t4'C -C _00. _, '-- * C')C' -- t* I£n *- t X - ^ I O O - , tJ I I * O O W -W . 4. II tN I3 O . , 0- 1- ,. X C3- Ew - tI4 , 'I - - I X t I C' tO- 1DtY _- . ~ Z t,' t L * Z- 1x X -. - .tN 1-4 U WZ II IA eI W U,~ I - I- _ I ), , 4D l ,_,, X >4 -I_ I~ _ O - 0, z 0 O - 140 0 0 ooooooooooooooooooooo0000000000000000 0000000000000000000000000000000000000 tn LU z 0 4 U, _ I z o § a. U I _. I * > c I C) _ Il4. _ >Z 0 u U Z t _ cl.QU -oO U. - CI - < < --< C---------- - <* . 4-Q- _ I .C 0 <O z >oo ,- F ''N' I 11 . - 1 I1 11 I I| ,'Q)t 0 iI V * - W .. _ -N4 + - -C4 C ' * a----------- W a4U n 0I 0 < _< - 4 *^ < -C _ CS- - <: w -W NC if _CN . -O\OuO\ .e wJzO - 11 - 11 11)1 - 11 '.4u zm 11- 11 11- 1 -' -*- It 1 It 11 ) C) o0 .) .. ..CO i F '-) - 0 , 11 0 11