Dept. of Naval Architecture and Marine Engineering Publication No. 233 January 1981 SURFACE SHIP PATH CONTROL USING MULTIVARIABLE INTEGRAL CONTROL Michael G. Parsons Hua Tu Cuong Department of Naval Architecture and Marine Engineering College of Engineering The University of Michigan Ann Arbor, Michigan 48109 This research was sponsored by NOAA Michigan Sea Grant #NA79AA-D-00093 Michigan Sea Grant Publications Office 2200 Bonisteel Blvd. Ann Arbor, Michigan 48109 Single copies free Bulk orders $2.00 ea. copy

Abstract The surface ship path control problem is formulated as a multivariable, linear state variable control problem subjected to measurement noise and nonzero mean disturbances. A multivariable generalization of integral control is presented and then specialized to the surface ship path control problem. The controller has the property of providing zero steady-state error to a constant commanded set point. It is insensitive to errors in the knowledge of the system characteristics. The controller has a nonzero steady-state error to a ramp commanded set point (nonzero heading straight path). This error is established analytically which allows its calculation in advance. The effect of the error can, therefore, be eliminated by simply shifting the time at which turns are initiated. The performance of the controller in straight steaming, lane changing maneuvers, passing maneuvers, and a series of turns and straight path segments is illustrated by digital simulations. The multivariable integral controller shows promise as an effective and practical surface ship path control concept. -i

Table of Content page Abstract i List of Figures iii List of Tables iv Nomenclature v 1. Introduction 1 2. Problem Formulation 4 2.1. Equations of Motion 4 2.2. Measurement Selection 8 2.3. Design Process Disturbances 10 3. Derivation of Integral Controller 14 3.1. General Derivation 14 3.2. Application to Ship Path Control 18 3.3. Steady-State Error to a Ramp Commanded Set Point 24 4. Controller Design and Evaluation 28 4.1. Controller Design 28 4.2. Verification of Need for C1Kyno Term 31 4.3. Performance with Bias Disturbances 35 4.4. Lane Changing Maneuvers 39 4.5. St. Marys River Turns 46 4.6. Cubic Turn Command 50 5. Turning by Coordinate Rotations 54 5.1. Derivation 54 5.2. Simulation Results 58 6. Conclusions 61 7. References 62

List of Figures Figure page 1. Coordinate System for Path Control 4 2. Design Passing Ship Disturbance 11 3. Design Lateral Current Disturbance 13 4. Schematic of Multivariable Integral Controller 19 5. Commanded and Actual Ship Paths in Nonzero Set Point Startup with Incomplete Control Law 32 6. Rudder Angle in Nonzero Set Point Startup with Incomplete Control Law 33 7. Commanded and Actual Ship Paths in Nonzero Set Point Startup with Complete Control Law 34 8. Rudder Angle in Nonzero Set Point Startup with Complete Control Law 34 9. Commanded and Actual Ship Paths with Design Lateral Current - Design for H/T = 1.89 Operating at H/T = X 36 10. Rudder Angle with Design Lateral Current - Design for H/T = 1.89 Operating at H/T = X 36 11. Commanded and Actual Ship Paths in Nonzero Set Point Startup with Step Disturbance 38 12. Rudder Angle in Nonzero Set Point Startup with Step Disturbance 38 13. Commanded and Actual Ship Paths in Lane Change with Passing Ship Disturbance 41 14. Path Error in Lane Change with Passing Ship Disturbance 41 15. Rudder Angle in Lane Change with Passing Ship Disturbance 42 16. Commanded and Actual Ship Paths in Lane Change with Design Lateral Current 43 17. Rudder Angle in Lane Change with Design Lateral Current 43 18. Commanded and Actual Ship Path in Single Turn with Design Lateral Current 45 19. Path Error in Single Turn with Design Lateral Current 45 20. Commanded and Actual Ship Path in Passing with Bias Disturbance 46 21. West Neebish Channel (from NOAA Chart 14883) 48 22. Commanded and Actual Ship Paths in St. Marys River Turns 49 23. Rudder Angle in St. Marys River Turns 50 24. Commanded and Actual Ship Paths with 370 Cubic Turn Command 52 25. Path Error with 370 Cubic Turn Command 52 26. Rudder Angle with 370 Cubic Turn Command 53 27. Rotation of Local Coordinate System 55 28. Global Coordinate Transformation 57 29. Global Path in 370 Turn by Coordinate Rotations 59 30. Local Path in 370 Turn by Coordinate Rotations 59 31. Rudder Angle in 370 Turn by Coordinate Rotations 60 -iii

List of Tables Table 1. Characteristics of Tokyo Maru Model and Prototype 7 2. Coefficients of Tokyo Maru versus H/T at Fn = 0.116 (12 knots full-scale) 7 3. Reference Measurement Noise Characteristics 10 4. Optimal Gains for Tokyo Maru at H/T = 1.89 and Fn 0.116 47 6. Sample Results in St. Marys River Turns Simulation 50 -iv

Nomenclature A state weighting matrix eq. (88) or system dynamics matrix eq. (73) Ai coordinate rotation point, Fig. 27 a slope in ramp lateral offset command aij element of matrix A B control weighting matrix eq. (88) or control distribution matrix eq. (73) bij elements of control weighting matrix B Ci element of Cx Cx optimal state feedback gain matrix Cy desired output feedback gain matrix di distance from Hi to Ai, Fig. 27 di+1 distance from Hi+1 to Ai, Fig. 27 E[...] expected value operator e error due to ramp input ess steady-state error due to ramp input F system open-loop dynamics matrix Fn=U/gW Froude number based on ship length fij element of system open-loop dynamics matrix F G system control distribution matrix or ship center of gravity, Fig. 25 g acceleration of gravity H water depth (m) or system measurement scaling matrix Hi projection of ship center of gravity on (i-axis, Fig. 27 Izz ship yaw mass moment of inertia [kgm2] I' =2I /pL5 nondimensional yaw mass moment of inertia zz zz J optimal control cost function Jzz yaw added mass moment of inertia [kgm2] J' =2J /pL5 nondimensional yaw added mass moment of inertia zz zz kij element of Kx ki element of Kv Kv output integral error state estimator gain matrix Kx Kalman-Bucy state estimator gain matrix Ky output due to steady disturbance estimator gain matrix L ship length between perpendiculars tm] or eq. (36)

element of L M matrix defined by eq. (56) m control vector dimension or ship mass [kg] or power in eq. (75) mij element of M m'=2m/pL3 nondimensional ship mass my sway added mass [kg] my'=2my/pL3 nondimensional sway added mass N total yaw moment or yaw moment disturbance [Nm] N'=2N/pL3U2 nondimensional N Nr derivative of yaw moment w.r.t. yaw angular velocity [Nms/rad] N'=2Nr/pL4U nondimensional Nr rv derivative of yaw moment w.r.t. lateral velocity [s] Nv derivative of yaw moment w.r.t. lateral velocity [Ns] N. derivative of yaw moment w.r.t. ldrift angle rad]tion Ns2] v NB derivative of yaw moment w.r.t. drift angle [Nm/rad] N'=2NO/pL3U2 nondimensional No Nj derivative of yaw moment w.r.t. drift angle rate [Nms/rad] N;'= 2N;/pL4U nondimensional Ni N6 derivative of yaw moment w.r.t. rudder angle rate [Nms/rad] Na=2N6/pL3U2 nondimensional N6 N- ~ derivative of yaw moment w.r.t. rudder angle rate [Nms/rad] n dimension of state vector or power in eq. (75) p dimension of measurement vector Q process disturbance power spectral density matrix q dimension of disturbance vector or disturbance spectral density qii diagonal element of Q R measurement noise power spectral density matrix r=d*/dt yaw angular velocity [rad/s] r'=rL/U nondimensional yaw angular velocity yaw angular acceleration [rad/s2] i'=rL2/U2 nondimensional yaw angular acceleration rii diagonal element of R s Laplace transform variable T output selection matrix T' output selection matrix for augmented state F Tr nondimensional rudder control time constant ~~t ~time [s]

ti time at coordinate rotation i t'=tU/L nondimensional time At integration time step U ship speed [m/s] u longitudinal component of ship speed [m/s] or scalar control u control vector steady-state control vector u' perturbation control vector v lateral component of ship speed [m/s] or scalar integral error state vc' nondimensional current velocity v integral error state, eq. (44) w process disturbance vector steady process disturbance vector w' zero-mean process disturbance vector X global coordinate XG global coordinate of ship center of gravity XH global coordinate of projection point H, Fig. 28 x longitudinal axis of ship x state vector Xi state vector of i-th system in coordinate rotations IXs steady-state state vector x' perturbation state vector Y total sway force or sway force disturbance [N] or global coordinate YG global coordinate of ship center of gravity global coordinate of projection point H, Fig. 28 Y'=2Y/pL2U2 nondimensional Y Yd Laplace transform of Yd Yr derivative of sway force w.r.t. yaw angular velocity [Ns/rad] Y'=2Yr/pL3U nondimensional Yr Y. derivative of sway force w.r.t. yaw angular acceleration r [Ns2/rad] Y.'=2Yr/pL4 nondimensional Y. r r Ys~ ~ Laplace transform of Ys Y,~~ ~ derivative of sway force w.r.t. lateral velocity [Ns/m] YB~ ~ derivative of sway force w.r.t. drift angle [N/rad] -vii

Y~=2Y8/pL2U2 nondimensional Yg uY derivative of sway force w.r.t. rudder angle [N/rad] YI=2Y6/pL2U2 nondimensional Y6 Y. derivative of sway force w.r.t. rudder angle rate [Ns/rad] y transverse axis of ship or scalar system output ]Y system output vector Ic constant of integration, eq. (46) Xd desired or commanded system output system output vector of i-th system in coordinate rotations difference between Yd and Yw, eq. (28) ~Ys=[s ss steady-state system output system output due to steady disturbance z measurement vector 4z measurement vector of i-th system in coordinate rotations coefficient in denominator of transfer function, eq. (75) or amount of coordinate rotation, Fig. 27 B=B' |drift angle [rad] relative to earth B' effective drift angle relative to water Bi coefficient in numerator of transfer function, eq. (75) f3 drift angular velocity [rad/s] o'=BL/U nondimensional drift angular velocity r system disturbance distribution matrix Yij element of matrix r 6=6' rudder angle (rad] 6c=6c' commanded rudder angle [rad] n lateral offset or transverse local coordinate [m] nondimensional lateral offset ni lateral local coordinate in system i ei,8i+1 angles defined in Fig. 27 Xi nondimensional eigenvalue v measurement noise vector Vi element of V longitudinal local coordinate _ augmented state vector S~i ~ longitudinal local coordinate in system i p density of water [kgm3] -viii

(aj RMS measurement j noise aj' measurement j noise standard deviation used in simulations T correlation time for measurement j noise angle between global X-axis and local ~i-axis, Fig. 27 heading angle Irad] (...) vector quantity derivative w.r.t. time (...) estimate of quantity (.7.) mean value of quantity initial value of quantity or acceptable value of quantity used in forming A and B weighting matrices [...]T transpose of a matrix [...]-1 inverse of a matrix -ix

1. Introduction The problem of controlling surface ships along prescribed paths in restricted waters is important from operational, safety, and environmental viewpoints. In the Great Lakes system, the difficulty associated with the safe movement of large bulk carriers through restricted waters such as the St. Marys River below the Soo Locks is a controlling factor in the evolution of larger, more economical vessels. The use of ships larger than the present 305 m (1000 ft) vessels may well be limited by a lack of maneuvering safety and/or excessive dredging costs. Since these larger bulk carriers would be a small, dedicated fleet, it might be practical for them to use onboard, microcomputer based automatic path controllers in the most restricted channels such as the St. Marys River. Precise, reliable, automated control might exceed the expected day-to-day performance of human operators and thus permit the use of larger, more economical vessels in more restricted channels with resulting reduced dredging costs and/or increased safety. These reduced system costs could more than offset the automatic control system costs. In this report, we investigate the feasibility and effectiveness of using a multivariable integral control law for the path control of a surface ship in restricted waters. This type of control would be a potential candidate for use onboard Great Lakes bulk carriers. These ships are subjected to short-term, essentially zero-mean disturbances due to passing ships, current and wind variations, waves, and bank and bottom changes. They are also subject to more long-term, non zero mean disturbances due to current, wind, second-order wave forces, and banks. The dynamic characteristics of the ships also change significantly depending on depth-under-keel, draft, trim, and speed. Maneuvering situations can place severe demands on the helmsmen and occur often in the Great Lakes system due to the high percentage of the voyage time spent in restricted waters. In previous work, we have investigated the feasibility and effectiveness of other control schemes for the path control of surface ships in restricted waters. Our earliest workl investigated the use of nonadaptive, optimal stochastic controllers for this purpose. These control systems consisted of a steady-state Kalman filter and a steady-state optimal state feedback controller. The Kalman filter uses noisy measurements to generate an unbiased estimate of the state which is then used by the controller to generate the -1

-2rudder command. In that work, the yawing moment and lateral force disturbances acting on the ship were modeled using first-order shaping filters. These controllers were shown to provide effective control when a ship is subjected to short-term, essentially zero-mean disturbances. That work did not produce controllers which could accommodate more long-term disturbances without a mean offset from the desired path. We also showed the desirability of an adaptive system which could automatically account for changes in the dynamic characteristics of the ship due to changes in depth-under-keel or other operating conditions. Our more recent work2'3 investigated surface ship path controllers which can accommodate long-term disturbances with a zero-mean offset from the path and which can adapt for changes in the dynamic characteristics of the ship. These control systems consisted of four major components arranged into two loops. The inner or control loop consisted of a steady-state Kalman filter and a steady-state optimal state feedback controller as used previously1 except that we modeled the yawing moment and lateral force disturbances using a brownian motion approach. The brownian motion approach assumes that the rate of change of the unknown disturbances is white noise of known spectral density. Standard optimal stochastic control techniques are then used to design a Kalman state estimator and state feedback controller. With the use of the brownian motion disturbance models, the Kalman filter can effectively estimate both the essentially constant and the stochastic disturbances which are acting on the ship at any time. The controller using the brownian motion disturbance model was shown to be very effective provided the ship hydrodynamic characteristics were well known. Since the dynamic characteristics of a ship will generally not be well known in restricted waters, our more recent work2'3 added a second, outer or gain update loop consisting of an on-line parameter estimator and a second function which recalculates the Kalman filter and controller gains using the latest estimates of the ship characteristics. A minimum variance parameter estimation scheme was utilized. This approach has shown some promise but has also shown the design concept to include seriously conflicting requirements. For the parameter estimator to be fully effective it is necessary to cause the ship to move dynamically about its desired path so that its rudder command and resulting motion histories can be used to estimate the hydrodynamic character

-3istics. This motion is, however, in direct conflict with the objective of precise path control. Compromise is therefore needed between the accuracy of the parameter estimates and the effectiveness of keeping the vessel on the desired track. A preferred approach for surface ship path control in restricted waters would be a scheme which would be effective with known ship characteristics but also "robust" or insensitive to changes in the ship characteristics from those upon which the design was based. The multivariable integral control approach studied here has these desired characteristics. Holley and Bryson4 recently completed a survey and evaluation of multivariable control techniques applicable to the automatic loading approach control of aircraft. This work has served as the starting point for our current work. Holley and Bryson concluded that a multivariable generalization of integral control provided effective control which allowed the zeroing of steady offsets due to essentially constant disturbances. Further, the resulting systems were insensitive to model (dynamic characteristics) errors as are present in the ship path control problem due to water depth, bank and speed changes if a nonadaptive controller is to be used. This approach does sacrifice some performance compared with the optimal stochastic controller using brownian motion disturbance models when the ship hydrodynamic characteristics are correctly known. This report is presented in four principal parts. First, the surface ship path control problem is formulated as a linear multivariable control problem. The selection of measurements and the definition of design process disturbances, which are used in the evaluation of system performance, are discussed. Second, the derivation of the multivariable integral controller is presented in general and then more specifically for the ship path control problem. The steady-state error for this controller to a ramp commanded set point is estimated using a deterministic approximation. Third, a multivariable integral path controller is designed for a specific ship and its performance is evaluated using digital simulation. Finally, a revised approach using coordinate system rotations is introduced and its performance is evaluated using digital simulation. The report closes with conclusions based upon this work.

2. Problem Formulation In this section, we formulate the surface ship path control problem as a linear, state-variable control problem. The selection of measurements and typical process disturbances are also discussed. 2.1 Equations of Motion. The development of the linearized, state-variable equations of motion for a surface ship moving in the horizontal plane presented here is based on the formulation by FujinoS and is presented in more detail in our earler work.1 The coordinate system for the problem is shown in Fig. 1. The 0 - Sn system is fixed in space with the desired ship path predominantly along the (-axis so that the prescribed lateral offset nd could be programmed as a function of This approach is typical of many maneuvering situations where the ship is to follow a series of straight paths or leading line segments along a general direction. A more general approach without restriction on the prescribed path will be utilized in Section 5. The G-xy system is fixed at the center of gravity of the ship. The positive sense of the drift angle B, heading angle, yaw rate r, and rudder angle 6 are shown. Neglecting the effects of pitch and roll, the ship motion can be described by coordinates x, y, and 4 desired path x,X"u U u = dx/dt v = dy/dt nd () / P U = (u2+v2)1/2 G r = d#/dt r,N yY,v Figure 1. Coordinate System for Path Control The exact equations of motion of the ship are integro-differential equations in which convolution integrals represent the memory effect of the fluid to previous motion.6 An alternative formulation yields differential equations -4

-5with frequency dependent coefficients. Fujino7 has shown that for the maneuvers of interest here the frequency dependence is negligible and constant-coefficient differential equations can be utilized. This assumption becomes less and less valid as the water depth to ship draft ratio H/T + 1. When the equations of motion are linearized about the nominal path, the equation in the x-coordinate decouples so that the ship motion can be given by, (m+my) d = YvV + (-mU+Yr)r + Yr + Y + Y, (1) dr (Izz+Jzz) dt = NvV + Nrr + N- v + N66 + N, (2) dt = U(S-)) (3) dt which are valid for small deviations from the nominal path, O0, and a constant speed U condition. An external sway force Y and an external yawing moment N are included to account for disturbances which act on the ship. It is common and convenient to utilize drift angle 8 instead of the lateral velocity v so we can use, v = -UsinO -=u_, (4) to express eq. (1) and (2) in terms of the drift angle. These equations can then be nondimensionalized as shown in the Nomenclature to yield, — l =" r', (5) dt' d~' -(m'+my')dt,: Y,'B + (-m'+Yr,)r' + Y r + Yd6d' + Y' (6) (I +J ) -Na'B' +Nrr' +N +' N'.6' + N' + N(7) zz zz dt' dt F=B' -' (8) dt' (6' - ) (9) dt' T c r where we harve now included a first-order model for the steering gear dynamics. The control is the comanded rudder angle dc'. The unit of nondimensional time t' is the time it takes the ship to travel one ship length.

-6Equations (5) through (9) can be transformed into state-variable form; i.e., 4,' 0 1 0 ~' 0 0 0 r' 0 f22 f23 0f25 r' 0 Y21 Y22 dt'' 0 f32 f33 0 f35 0'+ Y31 Y32, (10) dt'=' c nl' 1 0 -1 0 0 no o0 0. 0 0 0 0 -~1/Tr 1/T O 0 6'OL' I~~~1T or, nxl mx1 qxl = Fx + G u + r w. (11) The coefficients of the open loop dynamics matrix fij and the disturbance distribution matrix Yij are algebraic combinations of the stability derivatives and mass and inertia terms in eq. (6) and (7)*1 The multivariable integral controller of interest here can handle disturbances with a nonzero mean. We will therefore model these disturbances as the sum of two components; i.e., an unknown but constant part ws and an additive white noise disturbance w'. Equation (11) then becomes, = = Fx + Gu + rw + rw' (12) The problem thus has five states, one control, and four disturbance components. For a particular example, we will utilize the data obtained by Fujino5'8'9 for a model of the 290m (951 ft) tanker Tokyo Maru. This is obviously not a Great Lakes vessel but is utilized here due to the availability of data and to allow a direct comparison of the performance of the multivariable integral controller with that of the other approaches we have considered previously.l12 Systemmatic data for Great Lakes bulk carriers is only now becoming available through a Maritime Administration sponsored project at Stevens Institute and ARCTEC. Fujino conducted planar motion mechanism (PMM) and oblique tow tests of the model at various water depth-to-draft ratios H/T. Selected characteristics for this vessel are shown in Table 1. The coefficients fij and yij obtained for the Tokyo Maru at 12 knots full-scale at H/T values of 1.30, 1.50, 1.30, 1.50, 1.89, 2.50, and X are given in Table 2. As shown by FujinoS this vessel is course unstable for the intermediate depth-to-draft ratios from about 3.0 down to 1.75

-7as is typical of many large vessels. Fu ino' s characteristic model -prototype linear scale ratio, X 145.0 _ length between perpendiculars, m 2.000 290 breadth, m.3276 47.5 draft, m.1103 16.0 displacement 58.4 kg 179,100 LT block coefficient 0.8054 0.8054 rudder area 3,390.0 mm2 71.29 m2 propeller diameter 53.8 mm 7.80 m P/D 0.740 0.740 expanded area ratio 0.619 0.619 number of blades 5 5 Table 1. Characteristics of Tokyo Maru Model and Prototype. H/T 1.30 1.50 1.89 2.50 f22 l -1.6508 -1.7136 -1.7657 -1.8177 -1.9515 f23 9.3157 6.6235 5.7359 4.6112 3.1591 f25 -0.55543 -0.79235 -0.88074 -1.0416 -1.0410 Y21 346.69 385.98 477.68 536.00 567.13 Y22 4.8040 -2.2145 -5.0043 -5.8625 2.3365 { f32 1 0.02974 0.13890 0.17199 0.23621 0.31507 -1.0388 -0.71895 -0.52766 -0.54560 -0.63651 -0.09995 -0.12092 -0.15607 -0.'16639 -0.16163 vYi1 11.825 14.230 21.141 21.942 16.844 Y32 | -19.216 -23.123 -28.233 -31.490 -37.384 Table 2. Coefficients of Tokyo Maru versus H/T at Fn=0. 116 (12 knots full-scale) The output of the system, eq. (12), will be given by, ~~~~~~y =-~~~~~~~ T~x, (13) where T is the output selection matrix. We will take the lateral offset as

-8the output here so y is just a scalar and T becomes the row vector, T = [0, 0, 0, 1, 0] * (14) 2.2 Measurement Selection. All of the states in the ship path control problem as formulated in eq. (12) are available for measurement. The heading i' can be obtained from a compass; the yaw rate r' can be obtained from a rate gyro; the drift angle B' = -v' can be obtained from a doppler sonar; the rudder angle 6' can be obtained from the rudder stock or less accurately from the steering gear rams. The lateral offset from the desired path r' must be obtained using navigation aids such as DECCA Hi-Fix or radar. Each of these measurements may be subject to bias and zero-mean measurement errors and system transmission noise. In the presence of this measurement "noise" and with the measurement of only selected states, the complete state vector can be estimated using a Kalman filterll10'11 provided all of the states are observable with the chosen measurements. The authors have previously studied the observability of the ship path control problem.1 It was shown that it is necessary to measure the lateral offset n'. Additional measurements improve the ability of a Kalman filter to estimate all the states and thus improve the effectiveness of an optimal state feedback controller. The yaw rate r' is the next most effective measurement. The heading 4' is readily available and is the next most effective measurement. The drift angle 8' measurement was shown to add little to the effectiveness of a ship path controller which already measures nl' r', and 1'. With the steering gear model used in Section 2.1, there is little need to measure the rudder angle since the state is known exactly given any initial condition 6'(to) and the subsequent rudder command history 6'(t), t ) to. c In any practical application, the steering gear would have its own, separate feedback system; the first-order model included in eq. (12) is just a means of introducing a realistic rudder time response into our study. For the purposes of the controller design, it is reasonable to assume a measurement vector consisting of measurements of,' r', and n' each contaminated by Gaussian, white noise; i.e., pxl 1 0 0 0 0 v1 z = o 1 00 x + v2 Hx + v (15) 0 1 L 3

-9The final part of measurement definition is to establish reasonable levels for the measurement noise v. The white noise power spectral density needed in our continuous system design approach can be estimated by assuming the noise to be exponentially correlated with an RMS noise level aj and a correlation time Tj. The Tj should be much faster than the time constants of the ship and less than the system sampling time for the white noise model to be valid. The power spectral density can then be estimated by, rjj 2(j)2Tj. (16) To evaluate the control system effectiveness in this study, we use digital simulation with a fixed-stepsize Euler integration scheme. This has the effect of approximating the continuous Gauss-Markov process, eq. (12), eq. (14) and eq. (15), by a discrete Gauss-Markov process. In these simulations, the covariance of the computer generated random measurement noise must be selected to be consistent with the design noise power spectral density. To provide equivalent state estimate error covariances, it is necessary that the simulation measurement noise variance given by, aj'2 = L2,(17) At where At is the integration stepsize.l112 This can also be considered from a more direct viewpoint. If the controller is implemented digitally in an onboard computer with the system sampled at each At, the measurement noise will be a white sequence with variance aj'2. The reference measurement noise levels used in this work are shown in Table 3. In view of our earlier comments about the rudder model included in eq. (12), we assume exact knowledge of the rudder angle. Astrom and Kallstroml3 note that all sensors have dynamics with time constants less than 1 sec. and that the measurement errors are about 0.10 in 4, 0.020/s in r. Millers14 uses RMS errors of 0.2~ in 4, 0.010/s in r and 10 m. in n. Canner15 states that DECCA Hi-Fix crosstrack errors are as low as 1 m. when the baseline is along the desired path as is done at the entrance to Europoort. Astrom and Kallstroml3 and Bystrom and Kallstrom16 have found errors in r of less than 0.002~/s in systems identification of full-scale experiments. In view of this data, the reference levels assumed in Table 3 seem reasonable. The values for rjj and aj' are nondimensional. The aj' are calculated by eq. (17)

-10assuming At' = 0.005 (At -.24 s) which we use in our simulations. A noise correlation time of 0.1 s. and a sampling of.24 s. imply some correlation but the resulting values for aj' yield reasonable covariances for a white measurement noise sequence. Dimensionally, the aj' are about 93 percent of the assumed values for aj. measurement source aj Tj rjjaj compass 0.10 0. ls 1.298x10-8 1.611x10-3 r' rate gyro 0.001~/s 0.1s 2.860x10-7 7.563x10-3 DECCA Hi-Fix 3m 0.18s 4.559x10-7 9.549x10-3 Table 3. Reference Measurement Noise Characteristics 2.3 Design Process Disturbances. While operating in restricted waters, a ship can be subjected to a wide range of disturbances. Many of these can be characterized as being short-term relative to the time constants of the ship and as having essentially a zero mean value. First-order wave forces, wind gusts, and passing ships can be included in this category. Other disturbances remain long enough relative to the time constants of the ship that they must be considered to have nonzero mean value. Second-order wave forces and the effect of a lateral current, bank, or steady wind are included in this category. For the purposes of this study, we utilize two typical or design process disturbances in digital simulations to evaluate the performance of the path controller. These design disturbances were defined in our previous work.2 These definitions are repeated here for completeness. Passing Ship. The lateral force Y' and yawing moment N' due to a passing ship were selected as a typical short-term, essentially zero-mean disturbance. The assumed design disturbance is shown in Fig. 2. This disturbance is based on results originally presented by Newton17 for two Mariner vessels passing in deep water. These results are considered to be representative forces and moment histories and thus are reasonable for use in comparisons here. Yung18 and Abkowitz, Ashe, and Fortson19 show that the magnitude of the disturbances increase in shallow water as H/T + 1 so the magnitudes in Fig. 2 are known to be low at the shallower depths. In Fig. 2, the nondimensional

time scale is in ship lengths and the ships are beam-to-beam at t' = 0. These lateral force and yawing moment histories are assumed to be independent of depth under keel in our simulations. N'x105 y'x105 50 40 30 y' -2 -1 20. -02 Figure 2. Design Passing Ship Disturbance. Lateral Current. The effect of a lateral current was selected as a typical long-term, nonzero-mean disturbance for use in our ship path controller simulations. In a steady current, the ship assumes an equilibrium condition with 6' = 0 and,' = 8' so that the effective drift angle relative to the water Be' is zero. In this condition, there is no external hydrodynamic lateral force or yawing moment on the ship. In our ship maneuvering equations, we have assumed to this point that the drift angle B' is with respect to the earth. In shallow water, a doppler sonar would actually measure lateral velocity relative to the bottom. In a lateral current vc' without an additional disturbance, eq. (6), (7), and (8) should properly be written, -(m' + my') — = Y'Se' + (-m' + Yr,)r' + Y.' + Y6'6' (18) dt'

-12dr' (Iz + Jz dt' = NB'Be' + Nrwr' + N Be' + Np6', (19) dt'=' - Be' + VCI ~ Now if the drift angle relative to the earth 8' is introduced, Be' = B' + vc', (21) eq. (18), (19), and (20) become as follows in a steady current: (m y')d = YBB' + (-m' + Yr,)r' + Y~,r + Y66' + YBVc' (22) dt (I' + J' )- = NB'B' + Nr'r' + N. A' + N6,6' + NaBvc' f (23) zz zz dt' (24) dt' Thus when using the drift angle with respect to the earth in eq. (12), a steady current has the effect of applying an external lateral force and yawing moment given by, V = YfB'vc (25) and, N' = NB'vc'. (26) For design evaluation purposes, we have used eq. (25) and (26) to establish the lateral force and yawing moment produced by a 1 knot lateral current on the Toyko Maru moving at 12 knots in an intermediate water depth H/T = 1.89. This disturbance was assumed to be constant for 15 ship lengths and then to reduce linearly to one half this value at 20 ship lengths. This design disturbance is shown in Fig. 3. These lateral force and yawing moment histories are assumed to be independent of depth under keel in our simulations.

-13Y'/Y' Yo.0023277 N'/N o 0 N' =.0010262 5 10 15 20 25 30 t' Figure 3. Design Lateral Current Disturbance

3. Derivation of Integral Controller In this section, we derive the multivariable integral controller which is evaluated in subsequent sections. The general development follows that given by Holley and Bryson4 except that it is more general in that it can accommodate a nonzero initial command without a startup transient. The general development is followed by a specialization to the ship path control problem. The section closes with an estimation of the steady-state error of the controller to a ramp commanded set point. 3.1 General Derivation. The multivariable integral control law can be taken as a state variable feedback plus a feedback on a fictitious output yo; i.e., _ = CxX + Cy Yo ~ (27) At this point, o can be defined as the difference between the desired steady-state output Xd and Yw the output due to the constant disturbance Zo d Yw ~ (28) Substituting this expression into eq. (27), we arrive at an alternative expression for the control law; U = Cxx + Cy (Yd - Yw) ~ (29) Note that when the constant disturbance component _w is zero, its corresponding output Xw should also be zero, yo becomes identical to Xd and the control becomes, U(ws= 0) = CxX + CyYd. This provides a heuristic justification for the presence of the second term in the control law. The steady-state condition of the system in the presence of the constant portion of the disturbance ws can be used to derive an expression for the feedback gain matrix Cy. In the steady-state, we have x 3 0. We can designate the steady-state values of the state and control by the subscript s and the perturbations from these values by a prime; i.e., -14

-15x xS + xI t -- -s (31) us + u' ~ Substituting eq. (31) into eq. (12), we obtain the following system:' = Fx' + Gu' + rw', (32) O F + G + rws The steady-state solutions, if any exist, must then satisfy, (F + GCx)Xs + GCy Yo + rb = o. (33) If the system (F,G) is controllable, (F + GCX) is negative definite so, s = -(F + GCx)-(GCy Xo + rWs), (34) and the steady-state output is, Ys = TXs = -T(F + GCx) 1(GCy Yo + ris) (35) Defining, L = -T(F + GCx)' (36) this becomes, s = LGCy Xo + Lrws (37) If we now require that the steady-state output ys be equal to the desired output yd, comparison of eq. (37) and eq. (28) yields, Yw =Lr, (38) and LGCy= I, (39) where I is the identity matrix. Equations (36) and (39) allow the calculation of the gain matrix Cy from F, G, Cx, and T. Estimating Yw. In practice ys the nonzero bias component of the disturbance is rarely known. Therefore both the state x and Yw are unknowns in the control law, eq. (29). A recursive scheme can be utilized to provide an on-line estimate of Yw. Here we will present an argument which will show heuristically how one can arrive at a particular form for this estimation

scheme. Starting with a linear estimate approach: A A A w Ky(yw - Yw) = KyL(r - ws), (40) A A where ws is an estimate of w and Yw is an estimate of yw. From eq. (12), rFs is approximately equal to, rWs - x- Fx - Gu, or using eq. (27), rs! x - (F + GCx)x - GCy o ~(41) Substituting eq. (41) into eq. (40) yields, yw = - KyLC + KyL[(F + GCx)x + GCy Yo + rs], or using the definition of L, eq. (36), and eq. (39), A w KyL - Ky(Tx - Yo Lrws). (42) Recalling from eq. (37) that, Yd Yo + L rr, we can replace Is by ws in eq. (42) and it can then be rewritten as, Yw! - KyLk - Ky(Tx d).(43) We can now define a new state variable v by the differential equation, v= Tx yd ~ (44) This new state is thus the integral of the output error. Substituting eq. (44) into eq. (43), we obtain the following linear estimation scheme for Yw: Yw = Ky(Lx + v).(45) Equation (45) can be integrated to give, Yw = Ky(Lx + v) + Yc,(46)

-17where Yc is the constant of integration which can be obtained from the initial conditions of the system. In eq. (40), it can be seen that if the gain matrix Ky is selected to be negative definite with eigenvalues to the left in the complex plane compared with the other closed-loop eigenvalues, this estimator will provide a rapid estimate with yw + yw. Substituting eq. (46) for yw in the control law, eq. (29), yields, u'= Cxx + Cy Xd + Cv(Lx + v) - Cy Yc, (47) where, Cv = CyKy. (48) This final control law is comparable with that obtained by Holley and Bryson4 except that they omitted the final term. This additional term resulted from the constant of integration in the integration of eq. (45). It is necessary to allow the startup of the system with a nonzero desired output without an undesirable startup transient. This will be illustrated by simulation results presented in Section 4. Augmented System. In equations (12), (13), and (44), the system states x and the integral error states v are only available through the noisy measurements of the states, eq. (15). The system is also subjected to the process disturbances s and w'. In this situation, we can utilize the Separation Theoreml0 and estimate the states of the augmented system by the system of Kalman filters, 5x= Fx + Gu + Kx(z - Hx) (49) v = Tx Xd + KV(z - Hx) and then utilize these estimates in the control law, u Cx, + Cyyd + Cv(Lx + v) - Cyc. (50) From eq. (40) it can be seen that the eigenvalues of Ky will determine the dynamics of the estimate of Xw independent of x. These will also be eigenvalues of the closed-loop system. The other eigenvalues willbe those of the closed-loop controller (F + GCx) and those of the state estimator (F - KxH).

-18The remaining Kalman filter gains Kv can be determined to produce a zero mean output error Tx - Xd even when the system is subjected to measurement noise and process disturbances. Holley and Bryson4 show that this condition will result if, (T - KvH)(F - KxH)r 0. (51) This can be used to obtain Kv. Summary. The complete multivariable integral controller is defined by equations (49) and (50). The control gain Cx and state estimator gain Kx can be obtained using optimal control methods such as eigenvalue decomposition20,21,22 as used in our previous work. The gain matrix y can be obtained by pole placement. The matrix Cy is defined by eq. (39) with matrix L from eq. (36); i.e., the solution to, -T(F+GCx)-lGCy = I. (52) The matrix Cv can then be obtained from eq. (48). The matrix Kv can be obtained from eq. (51). This control concept will produce a zero steady-state output error with respect to a nonzero desired output when subjected to a constant disturbance. This result is independent of errors in the knowledge of the system matrices F, G, and r as will be present in the ship path control problem. The system can accommodate zero mean measurement noise and process disturbances. A schematic block diagram of the complete system is shown in Fig. 4. 3.2 Application to Ship Path Control The general form of the multivariable integral controller has been derived above. This can now be specialized to the surface ship path control problem as represented by equations (10), (14), and (15). We have the following system from above: = Fx + Gu + r + rw', z = Hx + v, (53) Y = Tw where,

-19disturbances W + W' +, actual output |F m' H desired output Yd controlled measurement+ v + z measurements plant noise control u C C C L _T controller initial condition y, integration constant state estimators Figure 4. Schematic of Multivariable Integral Controller

-200. 1 0 0 0 O 0f22 f22 3 f25 0 F = 0 f32 f33 0 f35, G= 0 1 0 -1 0 0 O o0 0 0 0 -l/T,/T'Y21 Y22 1 0 0 02 Y= 731 Y32, H = 0 1 0 0 0 0 0 00010 O O and, T = [ 0 0 1 0] Output Integral Error Estimator Gain Ky. In this problem, the output y.= n, the desired output d = nd, and therefore the output integral error v are all scalar quantities. With three measurements, the gain Kv is then a vector of dimension three which drives the estimator, A v = Tx- Yd + Kv(Z + Hx) * (54) This gain vector can be obtained from a system of two equations produced by, (T - KvH)(F - KxH)r -. (5i5) The solution to this equation is a one-dimensional subspace because we have three measurements of a system which is subject to two unknown disturbances. The state estimator gain Kx is in general a full matrix and we can designate its elements as kij. We have shown previously1 2 that for the system model used here the rudder angle is known exactly from knowledge of its initial value and the commanded rudder history. The last row of the state estimator gain matrix Kx is, therefore, zero for the statistical steady-state filter; i.e., k51 = k52 = k53 = 0. We can define the matrix, 5x2 M = [mij] = (F - XH)r, (56)

-21which has the elements, mil = (1 - k12)21 m12 ( - k2)Y22 m21 (f22 - k22)Y21 + f23Y31, m22 = (f22 - k22)Y22 + f23Y32 m31 = (f32 - k32)Y21 + f33Y31, m32 C (f32 - k32)Y22 + f33Y32 m41 = -k42Y21 - Y31, m41 -k42y22 - Y32 m51 = -k52y21 =, i52 = -k52Y22 = 0 If we now let the gain matrix Kv = [kl, k2, k3], eq. (55) becomes, [-k1, -k2, 0, 1 - k3, O]M =, which yields the following system of equations: -klm11 - k2m21 + (1 - k3)m41 = 0 (57) -klm12 - k2m22 + (1 - k3)m42 = 0 The resulting one-dimensional solution is, m4 1m22-m42m21 kI- (1 - k3) mi lm22-m 12m21 (58) ~ m1 m42-m 12m41 1 k2 =..-, -(1 - kI3) m1 m22-m1 2m2 1 where k3 is free to be selected by the designer. This indicates that an optimization process could be utilized to determine the value of k3 which would minimize the ITAE (integral of time multiplied by the absolute value of the error) or a similar transient response performance index.23 In this report, we will consider the simple special case where k3 = 1 so that eq. (58) yields, kl= k2 = 0 and, K= ('[0, 0, 1]. (59)

-22In this case, integral error state v has the definition, v n- nld, from eq. (44) and the estimator for this state, eq. (49), becomes simply v z3 nd * (60) Thus v is the integral of difference or error between measured lateral offset and the desired lateral offset. We therefore have a proportional plus integral control as desired. Feedback gains Cy and Cv. These two gains are used in the feedback control law, eq. (50). The gain Cy which multiplies the desired output Yd is a scalar. From eq. (50), the gain Cv will also be a scalar since u and v are scalars. Using eq. (39), we have, CY LG (61) where L and GT are each (1x5) vectors. The state feedback gain matrix Cx can be obtained using an optimal control technique such as eigenvalue decomposition which yields the statistical steady-state controller given F, G, a design cost function, and assumptions for the spectral densities of the measurement noise and process disturbances. This gain matrix can be represented here as Cx = [C1, C2, C3, C4, C5] ~ Using equations (36) and (53), the gain Cy is given by, 1 -[o 0 0 1 ol(F+GCx)-l O o. 1/Tr We therefore need only the (4,5) element of be shown to equal Tr/C4. The gain Cy th y - C4,q 1(62) and eq. (48) yields, Cv C4Ky. (63)

-23Integration Constant Yc. The scalar Yc appeared as the constant of integration in the integration of eq. (45). Substituting Yd = nd eq. (62), and eq. (63) into the control law, eq. (50), we have, ~~A ~A A u = Cx d - C4nd - C4Ky(Lx + v) + C4Yc. (64) A A A A A A Using, x = [(, r, O, n, 6]T and, L = [-1, ~2, ~3, ~4, ] this becomes, A A u = (C1 - C4yKyI1)I + (C2 - C4KyX2)r + (C3 C4Ky-3)S A A A A + (C5 - C4KKy5)5 + C4( - - - C4Ky~4n + C4Yc - C4KyV ~ (65) This equation permits a solution for Yc given any set of initial conditions on the states x(O) = x(O) and the control u(O). For the typical case where the ship is to begin on a straight course with an initial offset rn from the a-axis, we have, Yd = nd = no and, x(0) = x(0) = [0, 0, 0, no, O]T v(O) = v(0) = 0 Without a disturbance the commanded rudder angle control will also be zero, i.e., u(0) = 6c(0) = 0 Substituting these initial conditions into eq. (65) yields, Yc = Ky4Tno, (66) so that the subsequent control law becomes, A A A u = C x - C4nd - C4Ky(Lx + v) + C4Kyk4no. (67) From eq. (36) and eq. (53), ~4 can be seen to be just the negative of the (4,4) element of the matrix (F + GCX)-' which can be shown to yield,

-24C1 94 =-0- I C4 so that eq. (67) becomes, A A A u = Cx - C4nd - C4Ky(LX ) + v) + CKyo. (68) We noted earlier that Holley and Bryson4 omitted the last term in their development of this multivariable integral controller. Without this term, a startup of the system in an equilibrium position at no will produce a transient away from the equilibrium such that v will generate a value which will produce the same effect as the final term; i.e., Cl ro v=. (69) C4 We will show below that this startup transient can be quite large. Including the last term in eq. (68) avoids this undesirable startup transient. This completes the development of the multivariable integral path controller for a surface ship. In Section 4, this controller will be designed for the Tokyo Maru and its performance will be evaluated using digital simulation. 3.3 Steady-State Error to a Ramp Commanded Set Point In most practical situations, the prescribed ship path nd is defined by a series of straight lines or leading line segments. In general, most of these will not be parallel to the -axis and thus the commanded set point yd = nd will be a ramp function and not a constant. The controller developed above will produce a zero mean error for a constant commanded set point but its performance will deteriorate when yd is time-dependent. It is therefore of interest to study the common case where the commanded set point yd is a ramp yd = at', (70) where a is a constant. To establish the behavior of the output error of this controller when subjected to eq. (70), the closed-loop equations can be rearranged so that Yd appears as an input and then the steady-state error to a ramp input can be established.

-25The problem of interest here is stochastic with random measurement noise v and process disturbances w'. To simplify this section, however, we will treat only the deterministic problem v - w' = O. This greatly simplifies the evaluation of the error and reveals the essential character of this type of controller. For the deterministic case without a constant disturbance, we have, x Fx + Gu v= T~x - Yd, (71) y = Tx with the control law, u = C 4Yd - C4Ky(Lx + v) when no is taken as zero. Substituting this control into eq. (71) yields, x = (F + GCx - CG4KyL)x - GC4Kyv - GC4Yd v=Tx- Yd or, d [x][ (F + GCx - GC4KyL)(-GC4Ky)] [x] [ -GC4] T (72) where the prime for nondimensional time is omitted here. If we now define an augmented state variable 5 = [x, v]T and let T' = [0, 0, 0, 1, 0, 0], we have the system, ~ =A+By d, (73) y T', where A and B are defined by comparison with eq. (72). Equation (73) is a single-input, single-output system for which the transfer function can be obtained. If we let Y(s) and Yd(s) be the Laplace transforms of y(t) and Yd(t), respectively, eq. (73) can be expressed as, Y(s) = {T'[sI6 - A]'lB}Yd(s) = H(s)Yd(s), (74)

-26where I6 is the 6x6 identity matrix. The transfer function H(s) can be further expressed as, B0 + 1is +.oo + Bmsm H(s) =. (75) a0 + a1l +.. + ann The steady-state error e(t) due to a ramp input yd 3 at can be defined as, Yd(t)-Y(t) e(t) = lim t+00 a or, e(t) = Yd(t)-yss(t) (76) Chen24 shows that the steady-state response of a system with the transfer function eq. (75) to a ramp input is given by, o0 a0B1 - alB0 ySS(t) =-at + a a0 ao2 so that the steady-state error becomes, e(t) = 0 t - 1. - (77) cao aO Thus, if ao * Bo, the steady-state error to the ramp input is unbounded. If 0 = 0 * 0, the steady-state error reduces to the bounded result, e(t) = | - a1. (78) a0 Using eq. (74), the transfer function H(s) can be evaluated. This yields the following results for the first two coefficients in the denominator and numerator, respectively: a0 = Tr(f33f25 - f23f35)a56, (79) a1 = Tr[(f33f25 - f23f35)a54 + (f32f25 - f22f35 - f25)a56], (80) B0 = Tr(f33f25 - f23f35)a56, (81) B1 = C4(f33f25 - f23f35) + Tr(f32f25 - f22f35 - f25)a56, (82)

-27where a54 and a56 are the elements of matrix A in eq. (73); i.e., C4-C1Ky a54 = (83) Tr C4Ky a56 2 - - (84) Tr Note that as = 8o * 0 so the steady-state error to the ramp input is bounded. Substituting equations (79), (80), (82), (83), and (84) into eq. (78) yields, e(t) =. (85) C4 This, simple result is extremely important. This indicates that even though the multivariable integral controller has a steady-state error with a ramp commanded set point, this error is bounded and known in advance. The controller can, therefore, be programmed to compensate for this error by simply changing the point at which a turn is initiated. Notice also that a control law which minimizes C1, the feedback gain on the heading, will have a minimum steadystate error. The error can be eliminated with C1 = 0.

4. Controller Design and Evaluation In this section, we design a multivariable integral path controller for the Tokyo Maru using the characteristics defined in Section 2. The performance of this controller is then evaluated through a series of digital simulations. As noted in Section 2, the tanker Tokyo Maru is used because of the availability of data and to allow a more direct comparison of the performance of this controller with that of alternative concepts studied earlier. 4.1 Controller Design The characteristics of the 290 m tanker Tokyo Maru are given in Tables 1 and 2. The only undefined parameter is the rudder time constant Tr which we have taken as 10 seconds. In the nondimensional form used here this becomes, Tr = 0.21287, (86) at Fn = 0.116 or 12 knots full-scale. The first step in the design process is the calculation of the state feedback gain matrix Cx and the state estimation Kalman filter gain matrix Kx. These can be taken as the optimal steadystate solutions to the stochastic Linear Quadratic Gaussian (LQG) problem,2,10,25 k= Fx + Gu + rw (87) z = Hx + v where w and v are vector white noise processes with power spectral density matrices Q and R, respectively. The design cost functional is defined as the expected value of the integral, ft J = E[if (xTAx + uTBu)dt] (88) to where A and B are now weighting matrices which can be initially established by the designer to reflect the relative acceptability of errors in the various states and the use of the various controls. All four matrices Q, R, A, and B are diagonal with some nonzero elements. We have previously2'3 defined these matrices for the design of other optimal, stochastic controllers for the Tokyo Maru and use the same values -28

-29here. Detailed discussion of the choice of these quantities can be found in the earlier work so only brief comments will be included here. The power spectral density for the process disturbance Q is based upon the passing ship disturbance shown in Fig. 2. We calculated the root mean square values of N' and Y' between t' = -2 and t' = 1.4 and using an assumed first-order-process correlation time of one ship length for each disturbance established the diagonal terms, q11 = 1.548 x 10-8 (89) q22 = 8.970 x 10-8 The diagonal terms of the measurement noise power spectral density matrix R are given in the fifth column of Table 3. The nonzero diagonal terms of the A and B matrices were taken as, a44 = (~o)-2 = 772.5, a55 = (6o)-2 = 131.3, (90) B = b1l = (6co)-2 = 131.3, based on a dimensional use of 50 of rudder when the lateral offset error becomes 10.43 m (slightly less than one-quarter beam). The solution to this optimal, stochastic control problem requires the closed-form solution of two matrix Riccati equations2 which can be obtained from Potter's algorithm using eigenvector decomposition.20, 21 This technique was developed into a practical design tool by Bryson and Hall22 in their OPTSYS computer program. The Michigan Terminal System (MTS) version of the OPTSYS program has been used here to produce the state feedback gains C, and the Kalman filter gains Kx listed in Table 4. This design was developed for the characteristics of the Toyko Maru at the water depth to ship draft ratio H/T = 1.89, which is the least course stable depth for this ship. We have previously shownl that if a nonadaptive optimal, stochastic path controller is to be used, the best overall performance is obtained if the controller is designed for the ship's least course stable water depth.

-30controller gains CxT Kalman filter gains Kx 5.5421 4.6883 0.9507 0.0035 2.6601 20.9479 109.7887 -0.4755 6.3895 2.7730 9.0086 -8.6949 2.4252 0.1239 -0.7579 4.1275 -0.8499 0.0000 0.0000 0.0000 Table 4. Optimal Gains for Tokyo Maru at H/T = 1.89 and Fn = 0.116 The closed-loop eigenvalues for the Tokyo Maru at H/T = 1.89 when controlled by the optimal, stochastic controller given in Table 4 are as follows: X1, X2 = -0.52137 + 0.87033i A3 = -6.64361!4 = -0.97623 X5 = -2.32090 Recall from Section 3 that the gain Ky can be chosen by pole placement with respect to the eigenvalues of the closed-loop system (F + GCX). Here, we place Ky at the eigenvalue furthest to the left in the complex plane; i.e., Ky = -6.64361 ~ (91) This will cause the estimate of Yw (through v) to converge rapidly compared to the time response of the system. The output integral error gain Kv was taken as the simple form, Kv = [0, 0, 1], (92) as developed in Section 3.2. The value for the gain Cy can be obtained from eq. (62) and Table 4; i.e., Cy = -C4 = -2.4252. (93) Equation (63) then gives the gain Cv as, Cv = -C4Ky = 16.1121. (94) Finally, the steady-state error of the system to a ramp commanded set point can be evaluated numerically using eq. (85) and Table 4; i.e.,

-31C1 5.5421 es= - 1 - 1- - = 1 2.285 (95) This completes the design of the multivariable integral path controller for the Tokyo Maru. The performance of this design is evaluated in a series of digital simulations which follow. 4.2 Verification of Need for C1_Ky Term We noted in Sections 3.1 and 3.2 that the development of the multivariable integral controller presented by Holley and Bryson4 omitted the last term in the control law derived here. We derived, u =CxE - C4nd - C4Ky(Lx + v) + CiKyno, (68) for the special case where the ship begins controlled operation in equilibrium on a straight path offset no from the C-axis. The simulations described in this section were performed to show the general performance of the controller and to verify the need for the final term in u. We noted above that without this term, the ship would undergo an undesirable startup transient which could be prevented if the initial state were known. In both these simulations, the Tokyo Maru controlled by the multivariable integral controller designed in Section 4.1 is operating in a water depth H/T = 2.50. Thus, the ship is not operating in the water depth for which the controller was design and the simulations show the robustness of the design to errors in the knowledge of the dynamics of the system. We use an incorrect water depth here only as a mechanism to introduce a rational set of errors in the system dynamics. The errors in the system model could actually be due to any cause. At the water depth H/T = 2.50 the ship is still course unstable; the change in system coefficients between the design depth H/T = 1.89 and H/T = 2.50 can be seen in Table 2. The desired path is a straight line offset onehalf beam yd = no = nd = 0.0819 from the (-axis. The ship begins in equilibrium for this condition; i.e., x(O) = [0, 0, 0, 0.0819, O]T There is no initial disturbance but the ship is subjected to the design passing ship disturbance defined in Fig. 2 with the ships beam-to-beam at t' = 7. The measurements are contaminated with the measurement noise defined in Table 3 throughout the simulations.

-32Startup without the C1Kyno Term. The Tokyo Maru was simulated as described above using the control law eq. (68) but without the C1 Kno final term. The results of this simulation are illustrated by Figures 5 and 6. Without the final term the commanded rudder angle u' 6c at time t - 0+ was 3.016 rad. The actual rudder angle was driven to a maximum of 0.664 rad. which would exceed full rudder on ocean going ships and certainly exceed the validity of the linear modeling upon which the design is based. The resulting commanded and actual ship paths are shown in Fig. 5. (Note that the plotter routine has rounded the labels on the vertical scale; the line spacing is 0.008.) The startup transient caused by the lack of the ClKyno term in the control law causes the ship to deviate 0.076 or 22 m from the commanded course before returning to the commanded course. This transient is needed for the variable v to develop a value, eq. (69), which will compensate for the truncation of the control law. The influence of the passing ship disturbance on the path is lost in the final part of this transient. The rudder history during this simulation is shown in Fig. 6. The rudder activity due to the passing ship is evident at t'= 7. The rudder activity due to the measurement noise is about 10 after the ship returns to the commanded course..a I- - 4Am 12A0 *AD 24A 2U 36.00 TIME!. J 1==........__o~ i i XX.X \d L_____.__.._.0 W XXX/'W Incomplete Control Law

-33--— W -- ~ ~Control Law startup with Com ete Control Law. The Tokyo efaru was simulated as desecribed above using the complete control law eq. (68) to show the effect of and 8 The coanded 1200 00actual ship paths24. 20 2.in Fig 7 The ship beTIME Fig urinse. Rudder Angle in Nonzero Set Point Startup wi th Incomplete Control Law Startup with Complete Control Law. The Tokyo Maru was simulated as desecribed above using the complete control law eq. (68) to show the effect of the final term. The results of this simulation are illustrated by Figures 7 and 8. The commanded and actual ship paths are shown in Fig. 7. The ship beThe maximum deviation occurs at about t' = 8 due to time omparedssing ship disturbance. This deviation is less than 1 m. In comparing Fig. 7 with Fig. 5 note The integral controller provides effective control with zero mean disturbances such as the passing ship. This performance is at H/T = 2.50 and thus with incorrect knowledge of the dynamics of the ship. The final term in the control law correctly handles the effect of the nonzero initial set point. The rudder history during this simulation is shown in Fig. 8. The maximum rudder angle is about 7.50 in response to the passing ship. The rudder activity level in response to the measurement noise is again about +1~0 as would be expected.

-34S 1 4.0 4W 12.00 4.0 20.00 24.00 2.0 32.00 36. TIME Figure 7. Commanded and Actual Ship Paths in Nonzero Set Point Startup with Complete Control Law 0 S~~,fltt I I Figure8. Rudder Angle in Nonzero Set Point Startup with Complete Control Law

-354.3 Performance with a Bias Disturbances To evaluate the performance of the multivariable integral controller with a bias or nonzero mean disturbance, the Tokyo Maru was simulated under the control of the controller designed in Section 4.1 while operating in deep H/T = X water. The controller was designed for H/T = 1.89 so again the ship is operating with errors in the knowledge of the ship dynamics. The change in ship characteristics between H/T = 1.89 and the stable depth H/T = can be seen in Table 2. The simulation begins with the ship in a no lateral current equilibrium condition on the commanded straight path at nd = 0. The ship is subjected to the "design lateral current" disturbance shown in Fig. 3. This disturbance is constant for the first 15 ship lengths and then reduces linearly to half this value by 20 ship lengths. The disturbance then remains constant again after 20 ship lengths. The initial magnitude is established to be a one knot lateral current when the ship is operating in a water depth H/T = 1.89. The disturbance does not represent a true lateral current disturbance at any other water depth. When the simulation begins, the ship and controller are in effect subjected to a step change in yawing moment and lateral force. This simulation, therefore, represents a severe startup test for the controller. The results of this simulation are illustrated by Figures 9 and 10. Figure 9 shows the commanded (nd = 0) and actual ship paths. The maximum deviation due to the step change in disturbance is about 60.9 m or 1.3 beams. The controller then returns the ship to the commanded path. The maximum deviation due to the ramp change in disturbance beginning at t' = 15 is about 17.4 m or less than 0.4 beam. The controller returns the ship to the commanded path after the disturbance has stabilized. The rudder activity associated with this maneuver is shown in Fig. 10. The maximum rudder angle in the startup transient is about 330 which stretches the validity of the linear model but this stepchange in disturbance is a design test which is unrealistically severe. The equilibrium, mean rudder angles are about 6 = 0.131 prior to t' = 15 and 6 = 0.0655 after t' = 20. If the simulation had been conducted at H/T = 1.89 so that the disturbance would correspond to a true lateral current, these equilibrium, mean rudder angles would then have been zero and the equilibrium, mean state would have been E = f8, r = 0. The controller prowvides effective control when the ship is subjected to constant disturbances even when operating with errors in the'knowledge of the ship dynamics.

-36-; i-!....-tt —-— t — + —! l.. -l. ztz i <j*Iir - - I I I w ndz= =g............- - -. =- - 1 0.00 4.00.00 12.00 6.00 20.00 24.00 26.00 32.00 36.00 TIlE Figure 9. Commanded and Actual Ship Paths with Design Lateral Current - Design for H/T 1.89 Operating at H/T =. - - 1 T 11 \ 1., I I 1 1. / 1' —1 1 1. 1 |I oZO0 4.00 Ii.00 12.00 16.00 20.00 24.00 21.OQ 32.00 36.00 TIME Figure 1.o Rudder Angle with Design Lateral Current - Design for H/T 1 1.89 Operating at H/T = ~. -I 00 0.00/' 4.0 80'20.0 2.0 40 80 20 60 T~,l

-37In our previous work,2 the best performance for the simulation shown in Figures 9 and 10 was achieved with a controller which used brownian motion processes to model the unknown yawing moment and lateral force disturbances. This design was restricted to straight path nd = 0 operation. Comparison of Figures 18 and 19 of reference 2 with Figures 9 and 10, respectively, show comparable responses with the multivariable integral controller described here having somewhat superior performance. The maximum deviation at the startup are 60.9 m and 63.5 m and the maximum deviation due to the ramp change in disturbance are 17.4 m and 18.5 m for the multivariable integral controller and brownian motion disturbance model controller, respectively. This improvement in performance is achieved without any attempt at this point to optimize the transient response of the multivariable integral controller through the selection of k3 * 1 in eq. (58). This controller shows considerable promise and greater general capability than the brownian motion disturbance model controller studied earlier. To ensure that there is not an unacceptable change in the performance of the multivariable integral controller when subjected to a bias disturbance if the commanded set point is nonzero, the simulation described in Section 4.2 was repeated using a constant "lateral current" disturbance at a value corresponding to 1 knot at H/T = 1.89. As in Section 4.2, the ship was operating in a water depth H/T = 2.50; the complete control law eq. (68) was utilized. The commanded set point is a half beam offset, 1nd - 0.0819. The constant disturbance is applied in a step at the start of the simulation; i.e. the initial state and control are for the no disturbance condition. The results of this simulation are illustrated in Figures 11 and 12. Figure 11 shows the commanded and actual ship paths. The maximum deviation from the commanded path is about 42.9 m or about 0.9 beam soon after the startup. The rudder angle history for this simulation is shown in Fig. 12. The equilibrium, mean rudder angle in Fig. 12 is nonzero because again this design disturbance is not a true lateral current at the simulation water depth of H/T = 2.50. The controller provides effective control with no noticeable deterioration of performance with a nonzero commanded set point.

d 4l' - - - - - - ~l - _ —---- _....-... ~ w - - ~ —! 0".0o0 4.00 6.00 12.00 16.00 20.00 24.00 20.00 32.00 36.00 TIME Figure 11. Commanded and Actual Ship Paths in Nonzero Set Point Startup with Step Disturbance....~...!..........!1.....~~~~~~~,/ 0.0 1.0 1.0 MO 40 9 20 60 TIM

-394.4 Lane Changing Maneuvers A common maneuver in restricted waters is the change from one straight path to a second, parallel straight path. This can be extended to a second change back to the original path as would occur in a passing or overtaking maneuver. A series of lane changing maneuvers were simulated with various disturbances to evaluate the performance of the multivariable integral ship path controller. Lane Change with Passing Ship. In this simulation, the Tokyo Maru under the control of the multivariable integral controller designed in Section 4.1 was simulated to be operating at a water depth H/T = 1.30. The controller was designed for H/T = 1.89 so this simulation is with "incorrect" knowledge of the ship dynamics. The change in ship characteristics from H/T = 1.89 to H/T 1.30 can be seen in Table 21 the ship is very course stable at H/T = 1.30. The commanded ship path is nd = 0 for the first 10 ship lengths, varies linearly to 4 beams, nd =.655, at 20 ship lengths, and then remains constant at nd =.655 after t' = 20. The central portion of the maneuver represents a heading change of only 3.750. More extreme heading changes will be illustrated below. The simulation was initiated with the ship in equilibrium on the commanded path. The ship was subjected to the passing ship disturbance shown in Fig. 2 with the ships beam-to-beam at t' = 7. The results of this simulation are illustrated by Figures 13, 14, and 15. The commanded and actual ship paths are shown in Fig. 13. A small perturbation due to the passing ship is evident at about t' = 8; this only a little over 1 m. The ship turns to the transition course smoothly and completes the lane change with an overshoot of less than 2 m. The dominant feature of the response is the continuous error in the path during the transition phase. This is the steady-state error to a ramp commanded set point property of the controller studied in Section 3.3. The actual path lags the commanded path by about.15 or 43.5 m. These results compare exactly with the analytical results from Section 3.3. In Section 3.3, we found that for a ramp commanded offset, Yd = at' the steady-state error will be given by, C1 eSS = -4.

-40For the design developed in Section 4.1, the steady-state error was found to be, ess = 2.285 For the simulation shown in Fig. 13, the slope a is 4 beams in 10 ship lengths or a = 0.06552. Using the definition of the steady-state error, eq. (76), the steady offset error will be, Yd - y = 2.285a = 0.1497,(96) for this particular maneuver. These results correspond almost exactly with the simulation results. With this knowledge available in advance, the maneuver could simply be initiated 2.285 ship lengths "earlier" than shown in Fig. 13 and thus the steady-error in this maneuver can be eliminated as a practical concern. This can be illustrated further by Fig. 14 which shows the path error during the lane change shown in Fig. 13. The error during the transition can be seen to be about 0.15. Shown on Fig. 14 as a dashed line is the revised zero point for the path error if the turn were actually programmed to be initiated 2.285 ship lengths earlier than shown in Fig. 13. The zero point would then effectively be 0.1497 between t' = 10 and t' = 20 as shown on Fig. 14. (Actually the whole maneuver would be shifted 2.285 to the left.) There would be farily large errors as the ship completed the two heading changes smoothly over three ship lengths after t' = 10 and t' = 20. The error between t' = 13 and t' = 20 during the transition would then be 1.3 m or less. Recalling that the controller is operating with incorrect knowledge of the ship dynamics, this is highly effective path control.

-41d d w o0oo 4.00 o.00 12.00 16.00 20,00 24.00 2800 32.00 38.00 TIME Figure 13. Commanded and Actual Ship Paths in Lane Change with Passing Ship Disturbance ~. I I I I111 1 111111 ____''~-.....~'~ ~ — \ 3teady error.1497 4.00,00a 12.00 16.00 20.00 2400 2.00 32.00 ", TIME Figure 14. Path Error in Lane Change with Passing Ship Disturbance

The rudder activity in the maneuver shown in Fig. 13 is shown in Fig. 15. The maximum rudder is due to the passing ship disturbance and reaches a magnitude of about 8.30 just after the ships are beam-to-beam at t' = 7. The rudder angles to initiate the two heading changes reach magnitudes of about 5.20 and 4.00. The general level due to the measurement noise is within +1.50. 3 Ht ill, I -JIM,.[ A 1.f JA I A TIME Figure 15. Rudder Angle in Lane Change with Passing Ship Disturbance Lane Change with Bias Disturbances. This simulation is similar to that just described except that (1) the ship is subjected to the design lateral current disturbance shown in Fig. 3, (2) the ship is operating in a water depth H/T = 2.50, and (3) the lane change is for a 0.5 ship length offset in 10 ship lengths. In this situation, the commanded set point slope a = 0.0500 during the transition and the resulting offset error corresponding to eq. (96) is 0.1143. The results of this simulation are illustrated in Figures 16 and 17. The commanded and actual ship paths are shown in Fig. 16. The corresponding rudder activity iu shown in Fig. 17. H/.5,ad 3)telnchneifoa0.5.hp.egh.fse.n10si lengths.~ Inti iuto,Itecmaddstpit'lp.50drn the~ ~ ~ ~ ~~~~V trnito and/ the/ reulin ofstero orepnin*o'. 9)i O.14.Th esls fthssiuato ae lusrte n iurs16ad 7 The~~~~~~~~ ~~~ comne an*culsi ah r hw nFg 6 h orsodn rudder~ ~ { acivt,s Vhw inFg 7

-43c -- -, — * CTI -d J I- ~- - -1 0.00 4.00.00 12.00 1600 20.00 24.00 2.00 32.00 36.00 TIME

-44As in the simulation with the design lateral current illustrated in Figures 9 and 10, the ship is subjected to a step change in yawing moment and lateral force disturbance at the start of the simulation. The ship path shown in Fig. 16 therefore, includes a transient from the commanded path which reaches a maximum offset of about 43.5 m or 0.9 beam at about t' = 3. The ship enters the transition effectively and the path lags the commanded path by about 2.285 ship lengths with a path error of about 0.1143 during the constant disturbance prior to t' = 15. The' period of time-varying disturbance between t' = 15 and t' = 20, however, appears to introduce an additional time lag of about 0.7 ship lengths and introduce additional path error. If the turn had been programmed to be initiated 2.285 ship lengths earlier to eliminate the steady-state error prior to t' = 15, the resulting cross-track error at t' = 20 would have been about 11.8 m or one-quarter beam. The additional path error introduced by the time-varying disturbance between t' = 15 and t' = 20 appears in Fig. 16 to be decreasing near the end of the transition. To ensure that the time-varying disturbance does not introduce a change in the steady-state error to a ramp commanded set point, we repeated the simulation shown in Figures 16 and 17 with the second turn at t' = 20 eliminated. The resulting commanded and actual ship paths are shown in Fig. 18; the path error is shown in Fig. 19. The time-varying disturbance can be seen to perturb the ship path but not alter the steady-state error properties. The maximum cross-track error would be only 11.8 m or one-quarter beam at t' =20.4 if the turn were programmed 2.285 ship lengths early to accommodate the known steady-state error. The controller is, therefore, effective with both bias disturbances and large time-varying disturbances as might be expected from banks or current changes.

0.00 4.00 800 2.00 16.00 20.00 24.00 28.00 32.00 36.00 TIME Figure 18. Commanded and Actual Ship Paths in Single Turn with Design Lateral Current cr~ -0-o 12.00 - 0.00 I I.00 2.00 6.00 -- _ ~, - -:..... T. TIME Figure 19. Path Error in Single Turn with Design Lateral Current Lateral Current ~ ~ ~'st II I....11.1..1 1 1 1 -

-46Passing with a Bias Disturbance. As a final illustration of the effectiveness of the multivariable integral controller with a bias disturbance, the Tokyo Maru under the control of the design developed in Section 4.1 was simulated to pass another ship using two 0.5 ship length lateral transfers requiring 10 ship lengths each. The transfers were commanded beginning at t' = 5 and t' = 25. This simulation was conducted in a water depth H/T = 2.50. The ship was subjected to yawing moment and lateral force step disturbances with magnitudes corresponding to a 1 knot lateral current in H/T = 1.89; i.e., Fig. 3 with t' < 15. The commanded and actual ship paths are shown in Fig. 20. Again the performance is very good. The time lag and steady-error properties can be seen to be independent of turn direction and bias disturbance magnitude. Each of the four turns could be initiated 2.258 ship lengths early to eliminate the steady-error during the two transition periods. The overshoot at t' = 19 and t' 39 is about 3 m. 11~SF a e W 00 00 0 28 00 Al3,00 Figure 20. Commanded and Actual Ship Paths in Passing with Bias Disturbance 4.5 St. Marus River Turns I-~- - -~-~\

-47of the 304.8 m (1000 ft) Great Lakes ore carriers,26 the series of three turns in the St. Marys River between Sand Island and Moon Island in the West Neebish Channel was selected as the prototype for this simulation. These turns were identified as some of the most difficult in the Great Lakes system for the large bulk carriers. This path is the downbound lane and is shown in Fig. 21. The second leg of this path is the Rock Cut which has vertical cut stone walls. The channel at this point is roughly three beams wide for the largest ships in the system. For the purposes of the simulation reported here, only the path was utilized from the prototype. The simulation maneuver is defined in Table 5. segment length heading turn at end of segment 1 13.38 4 = 0 370 turn to port 2 14.76 4 _ -37~ 280 turn to starboard 3 13.79 4,= 90 400 turn to port 4 8.07 4 = -490 Table 5. Simulation Maneuver Based Upon St. Marys River For-this simulation, we continued to utilize the Tokyo Maru under the control of the controller developed in Section 4.1. The simulation was performed at a constant depth of H/T = 1.89. Thus, the operating and controller design depths were the same. Bank effects and current were not included in this particular simulation. The passing ship disturbance shown in Fig. 2 with the ship beam-to-beam at t' = 7 was included for test purposes even though the prototype is actually a single direction channel. The simulation was started in equilibrium on the commanded path. The results of this simulation are illustrated by Figures 22 and 23. The commanded and actual ship paths are shown in Fig. 22. Because of the coordinate system defined in Fig. 1, the plot of n versus t' appears as a mirror image of the actual path. For clarity of the results, the turns were not initiated 2.285 ship lengths "early" but this could have been done to eliminate most of the path error in the final three segments when the commanded heading is not zero. The analytical path error Yd - Yss for the second, third and fourth segments are 1.722,.362, and 2.629, respectively, when the turns are not programmed to be initiated 2.285 ship lengths early.

-48-424 4;a's ww D.3~;~1 P 1 2 4 ~ -..2.. 2I':. MII \ D I~ 2 3 ('33 I't ~,,\4P~ OkidtsVI' 6~d CLO 4 2~~~~~~~~~~~K~ 8 IM~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~~~~61 I 3 2 6 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~1 1'~ 3 ~~~~~~~~~~~~~~~~~~~~~~2' C 4.12 2 ~ S~,~-~ I 1 0~~~~~~~~~~~~~~~,,,IMr~RIc,p~ ~IqpkotiI 3 33 4 6' E 31 Co" "36"1~ 1,000/, -~~~~~1,R4 I J~ 0~4 2 A~~~~~~~~~~~~~~~~~~~Ar \/ /=\ "F,:,, M.'. \' _~~~~~~~~~~~~~~~~~~~ 9UnIIIinL OkRdq~c11~jr 3'SSrs/' I 0! \i 2 3 0 2 Ri,6 Oc\\,60A Na.0::..,.., 0 > l:' Alwlp~~~~~~~~~~~~~~~~~~~~~~l ~~~~~~~~~~~~~~~~~~c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c Fiur 21. Wes NeebshCannl (f N Cr 4 %~~~~~~~~~ _ 3 2 2 3~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~I/ \Oq part 0~~~~~~~~~~~~~Rt'-'~x, ~ ~~~~~~~~~~ 17 q~C 5 ~ ~~~ i 2 4I 2 6\ / 8 r~~~~~~~~~~~~~~~~~~~~~~~~~...,,=;,,,,:;~'~ =~~~~~~~~~~~ 2 2 -~. —.-. —,_-=4'''I 1 2,0 3.'d:'"'",'I 4 M 13 3 3 4 6 17 - 3 4 PI 4 &a#" 2 1 3 6 2 26 " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~,,Ok\as, cr 2 1 3 3 1 1 3 5 5~~~~~~~~~~~~ 2 4 2 9 r I 2 2 3 5 C1 4t~~~~~~~~~~~I;, ~,' " 1 ~~r'L1 -r~~~~~r 6 On r \4 I~~9 7 3 31 ~ - 2 3 5 9\A3 4,a' Pk L~-l* ~LI Lnin Wlwl3 P~" ~~~~~~~~~~~i~~~~~~~~~~~~~~~~~~~,i',,b-'.(,.s' I, " ryrlror~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~/".'",.,, 8 ----------— ~ I Fiue2.I W e s N e bs C h ne ( fo NOA Chart f ~~~~~~-;~~~rtt II~1483

-49Sample results from the simulation are shown in Table 6. The resulting path error with the 2.285 shift in the turn commands is shown in this Table. In general, the performance is excellent. The larger error at t' = 50 is still only.15 beam for this ship. This error seems to result from the presence of the transient initiated at the start of the segment since the error at t' = 47, 48, 49, and 50 is +.025, +.004, -.025, -0.25, respectively. Allowing for the shift of 2.285; the maximum overshoot after the first turn is.094 or 27.3 m. This corresponds to a crosstrack overshoot of 21.8 m or.46 beam. The maximum overshoot after the second turn is.082 which corresponds to a crosstrack overshoot of 23.5 m or.49 beam. The maximum overshoot after the third turn is.121 which corresponds to a crosstrack overshoot of 23.0 m or.48 beam. The rudder angle in these turns is shown in Fig. 23. The maximum rudder angle at the start of the first turn is 23.90; maximum rudder angle at the start of the third turn is 30.00. In the next section, we investigate an approach for reducing both the overshoot and maximum rudder angles. 8i'0.00.00 16.00 24.00 32.00 40.00 4 56.00 64.00 n.oo TI1E Figure 22. Commanded and Actual Ship Paths in St. Marys River Turns

-50-i TIME Figure 23. Rudder Angle in St. Marys River Turns simulation analytical path error with Yd nd QYd Yd -Yss 2.285 shift 26 -9.506 -7.787 1.719 1.722.003 = 0.9 m 28 -11.013 -9.295 1.718 1.722.004 = 1.2 m 39 -12.842 -12.48.362~.005.362.000 = 0.0 m 41 -13.158 -12.80.358~.005.362.004 = 1.2 m 48 -20.283 -17.65 2.633~.005 2.629.004 = 1.2 m 50 -22.584 -19.98 2.604~.005 2.629.025 = 7.3 m *printed simulation results have only four significant figures so results are ~+.005 Table 6. Sample Results in St. Marys River Turns Simulation 4.6 Cubic Turn Command In order to reduce the path overshoot and reduce the maximum rudder angles the commanded turn can be a smooth Curve rather than a discrete change of head

-51ing as used in Fig. 22. To investigate the effectiveness of this approach, the first turn in the St. Marys River turns simulation described in Table 5 with the passing ship disturbance shown in Fig. 2 was repeated using a cubic turn command. For simulation convenience, however, the 370 turn was made to starboard. The commanded path was as follows: Yd=0 t' ( 10 Yd -.0086704 (t'-10)3 +.0056385(t'-10)2, 10 ( t' ( 15.17 Yd =.75355(t'-13.38) 15.17 t' This cubic transition path reaches the desired 370 path beginning at t' = 13.38 at t' = 15.17 with the correct 370 heading. The commanded and actual ship paths for this turn are shown in Fig. 24. As in Fig. 22, the turn is not programmed to be initiated 2.258 ship lengths early in order to eliminate the steady-state error. Analysis of the simulation results allowing for the 2.258 shift in-the actual path shows that the response has a maximum overshoot of about.0511 at t' = 18.3. This overshoot can be seen in Fig. 25 which shows the path error for the turn shown in Fig. 24. With the 2.258 ship length shift, the effective path error after the turn is completed would be with respect to the steady-state error value of 1.722. Dimensionally this overshoot is 14.8 m compared with the 27.3 m overshoot in the first turn of Fig. 22; the crosstrack overshoot is 11,8 m or one-quarter beam compared with 21.8 m in Fig. 22. The cubic turn command, therefore, provides a significant reduction in path overshoot. The rudder activity in the turn shown in Fig. 24 is shown in Fig. 26. The maximum rudder angle magnitude in the initation of the turn is 5.60; the maximum rudder angle used in checking the turn is 8.90. These values compare favorably with the 23.90 and 14.70, respectively, shown in Fig. 23 for the first turn. The cubic turn command, therefore, provides a significant reduction in the rudder angle magnitudes used in turns made under the control of the multivariable integral controller. Controller implementation of this control law could include both the 2.258 ship length shift and cubic turn commands. The resulting controller shows considerable promise.

o~~8 - ~5200 4.00 8 12.00 16.00 20. 24.00 28.00 32.00 36.00 TIME Figure 24. Commanded and Actual Ship Paths with 370 Cubic Turn Command ___ _ ->' I -II | 1 1 1 1 1 ~ 1 1 ~1 - 1:;\|'fady errm 1. 7 2 2 0.00 4.00.00 12.8 16.00 20.00 24.8 20.00 32.8 36.00 TIME g~~~~~~~~~.....

-53a 0.00 4.00 0.00 2.00 16.00 20.00 24.00 2. 32.00 36.00 TIME Figure 26. Rudder Angle with 370 Cubic Turn Command

5. Turning by Coordinate Rotations In the previous section, a cubic turn command was utilized to reduce rudder usage and overshoot in turns under the control of the multivariable integral controller derived in Section 3. In this section, we study an alternative approach in which the turn is achieved by a series of coordinate rotations within the control computer without a specific turn command. The necessary coordinate transformations are developed and then the effectiveness of this approach is illustrated by digital simulation. 5.1 Derivation In the results presented in Section 4, turns were at discrete locations and were initiated as discrete changes in heading; i.e., ramps in commanded offset, or were initiated using cubic transitions in the commanded offset followed by a ramp commanded offset. This approach included a steady-state error to the ramp commanded offset which will occur with any constant nonzero heading.'Knowing the steady-state error in advance we noted that this could be taken into account by simply initiating the turns 2.258 ship lengths early for the specific design studied. This approach would be practical for a series of discrete turns with intermediate periods of constant heading as in the St. Marys River turns simulation. In more complex situations involving more closely spaced turns or perhaps continuous turning, this would become difficult to manage. An alternative approach would be to rotate the coordinate system used by the controller so that it is aligned with the desired heading. This could be implemented in small, discrete steps consistent with the cycle time of the control computer. Since the final condition will always have the nd = 0 with respect to the internal coordinate system, there will be no steady-state error. In this approach, we will be dealing with two coordinate systems. The global system (X,Y) is fixed to the earth. It is utilized to keep track of the position of the ship with respect to the navigational charts. The local system (,Bn) is generally aligned with the channel with the E-axis along the centerline or desired track. The control law will be based on the local system since the state vector will be computed with respect to the local system. Rotation of Local Coordinates. The geometry associated with a discrete rotation of the local coordinate system at a time ti is shown in Fig. 27. Point G is the center of gravity of the ship. The desired path can be construc

-55ted as a series of straight segments separated by the rotation points Ai. At the rotation time ti, the local coordinate system ( ni, fi) is at an angle *i with respect to the global X-direction. The local lateral offset is ni = HiG in Fig. 27. At ti, the local coordinate system is rotated through an angle ai to a new local coordinate system (Fi+1, ni+l) which is at an angle *i+1 with respect to the global X-direction. The new local lateral offset is ni+l = Hi+1G ~ The angles Oi and 8i+1 locate the rotation point Ai with respect to the ship using the ni- and ni+l-axes, respectively. A reasonable logic for establishing the rotation time ti would be when 8i = 28i+1 for the particular geometry shown in Fig. 27. Suitable equivalents could be developed for the other combinations of ai and ni. ship path -/ commanded path a, Figure 27. Rotation of Local Coordinate System If we define GHj positive in the direction of positive ni and define di = HiAi and di+i = Hi+jAi positive in the direction of positive Ei and i+1 respectively, we have, di.i = tan1 di(97)'-i i+1 = tan1 di+- (98) and,

-56ai = 8Oi Oi+1 (99) The rotation time ti might then be taken when Oi+1 = ai for the geometry shown in Fig. 27. After the rotation, the new lateral offset is, ni+l = -GAi cos 8i+1 = - /; + di2 cos(8i - i ) (100) The new heading angle *i+l will be related to the prerotation heading angle by, *i+1 =4 i + ci (101) The other states will be unchanged so the state vectors before, x[, and after, Xi+1, the rotation at ti will be given by,,xi| = ([ij r, no, n, iT, 1 (102) Xti XXi+1 (103) Xi+l| = [i+1, r, 8, ni+,1 6]T, ti with the new ni+l and 4i+1 given by eq. (100) and eq. (101), respectively. At ti, the system model is switched from the system i; i.e., Xi = Fxi + Gu + rw!i = Hxi + V, yi = Txi, to the system i+1; i.e., Xi+l = FXi+l + Gu + rw +l - Hxi+l + V Yi+1 = Txi+1 Since the state vector xi(ti) contains all past history of the system prior to ti, xi+l(ti) from eq. (103) can be used as the initial condition for system i+1. The validity of the switching process is guaranteed by the SemiGroup axiom of system theory.27 Finally, the error integral state v in the integral controller should continue undisturbed; i.e.,

-57rti ti+1 v(t) =... + (nli fndi)dt' + (fi +1 ndi+l)dt' +. ti-1 ti Global Coordinate Transformation. With the rotation of the local coordinate system it is necessary to have a transformation from the local to global coordinates to be able to keep track of the position of the vessel. The required geometry is illustrated in Fig. 28. The ship center of gravity G and projection H on the Ei+I-axis are shown for the rotation time ti and a later time t. If we assume that the ship path is never at a large heading with respect to Fi+1 and that speed is constant as in the initial development of the linear model, the nondimensional time t will have units of ship v A ship path G(t}/ Ei+l / \ com manded path G(ti ~~~~~~O K +1 X rli+l Figure 28. Global Coordinate Transformation lengths along the Ei+1 axis. At the time of the rotation ti, the distance di+1 from Hi+l to Ai is established. Thus, the coordinates of Hi+l in the global system (Xi,Yi) can be established. The position of H(t) is then given by, XH(t) = (t - ti)cos *i+1 + Xi, YH(t) = (t - ti)sin qi+1 + Yi The position of the ship in the global coordinates is then given by,

-58XG(t) = XH(t) + ni+lsinji+1 (105) YG(t) = YH(t) - ni+1cosi+1 ~ In a treatment using the full nonlinear equations of motion the F-coordinate equation could be utilized and speed changes could be accounted for. The treatment here is restricted to the linear, constant speed model. The prime has been dropped from the nondimensional time in the above development. 5.2 Simulation Results In this example, the Tokyo Maru under the control of the multivariable integral controller designed in Section 4.1 is simulated through a 370 turn which is initiated using the coordinate rotation approach. The ship is operating at the controller design depth of H/T = 1.89. The commanded set point Yd = nd 0 O throughout since the turn is created by the coordinate rotations. A single rotation of 370 would produce unrealistically large rudder commands. The turn is therefore created in a series of 40 equal rotations at each 0.1 ship length. The turn begins at t' = 11.38 and is thus completed at t' = 15.38. The passing ship disturbance shown in Fig. 2 with the ships beam-to-beam at t' = 7 is included in this simulation. The results are illustrated by Figures 29, 30, and 31. Figure 29 shows the global path for the turn. The 370 heading final track begins at Y = 0 at X = 13.38 as in the first turn of the St. Marys River turns simulation. The local commanded and actual ship paths are shown in Fig. 30. In this approach, the commanded path is always nd = 0 so the local path represents the path error. There is no steady-state error in this approach. The maximum overshoot occurs at t' = 17.5 with an overshoot of.1526 or 44.3 m. This is a crosstrack overshoot of 35.3 m or.74 beam compared with 11.8 m for the cubic transition turn shown in Fig. 24 and 21.8 m for the discrete 370 first turn shown in Fig. 22. The rudder angle in the turn is shown in Fig. 31. Using coordinate rotations, the rudder magnitude reaches about 280 at two points in the turn. This compares unfavorably with the maximum rudder angle magnitude of 8.90 shown in Fig. 26 for the cubic transition turn and 23.90 shown for the first turn in Fig. 23. The coordinate rotation approach does not offer improved turning performance. It could be used periodically, however, to realign the local coordinate system and eliminate the steady-state error.

......Ft..., ~GLOBAL X Figure 29. Global Path in 370 Turn by Coordinate Rotations Cr C.oo 4.00 8. 12.00 16.00 20.00 24.00 28.00 32.00 36.00 TIME Figure 2930. Local Path in 370 Turn by Coordinate Rotations _~ —1~~~~,!,:w............... -,,- - - ~ - i.- - - ~ 1 l... Figure 30. Local Path in 370 Turn by Coordinate Rotations

-600 o ---' 00 4.00 8.0 12.0 16. 20. 24.00 2. 32.00 36.0 TIME Figure 31. Rudder Angle in 370 Turn by Coordinate Rotations

6. Conclusions The following are the principal conclusions based upon this work: * The multivariable generalization of the integral controller presented here follows that of Holley and Bryson4 except that we have obtained an additional term in the control law which results from a constant of integration. This additional term in the control law allows the controller to accommodate nonzero initial setpoint commands without a highly undesirable startup transient. * In the specialization of this controller to ship path control, we have studied at this time only the special case where the integral error variable Kalman filter gain is Kv = [0, 0, 1]. The design of this gain contains the free variable k3 which can be selected to optimize the transient response of the controller with the other elements of Kv = [kl, k2, k3] then obtained using eq. (58). We hope to investigate the effect of optimizing k3 on the performance of this controller in the near future. The multivariable integral controller has the property of zero steadystate error with a constant commanded set point when subjected to disturbances and measurement noise. In ship path control, a common situation is a nonzero heading straight path which corresponds to a ramp set point command. In this case, the multivariable integral controller has a nonzero steady-state error which can be interpreted as a time shift in the turn response. We derive a simple analytical expression for this error, eq. (85), which allows its calculation in advance. The effect of this error can then be eliminated by simply initiating the turn a fixed time earlier than would normally be expected. In simulation results presented in Section 4, we show the value of the additional term in the control law in eliminating unwanted startup transients if the initial commanded set point is not zero. With the complete control law, performance with a constant set point is excellent even when the ship is subjected to large time-varying disturbances and when the design is based upon incorrect knowledge of the characteristics of the ship. The nonoptimized controller studied here shows superior performance to the path controllers we have developed earlier2 using brownian motion disturbance models. The multivariable controller provides effective control in lane changing and passing maneuvers provided the time shift needed to offset the steady-state -61

-62error to a ramp commanded set point is implemented. This performance is not sensitive to errors in the knowledge of the dynamics of the ship. The controller is effective when subjected to large bias or time-varying disturbances. The multivariable controller provides effective control in larger magnitude turns between straight path segments as included in the St. Marys River turns maneuver defined in Table 5. The time shift must be included to offset the steady-state error to ramp commanded set points when the path heading is nonzero. Cubic transition set point commands can be introduced to reduce the overshoot and rudder activity in the turns. The use of coordinate rotations was introduced as an alternative means of eliminating the steady-state error to a ramp commanded set point and possibly reducing the overshoot and rudder activity in the turns. The latter objective was not realized, but this approach could be used periodically to "update" the local coordinate system in the controller and to accommodate continuous turning situations.

7. References 1. Parsons, M.G., and Cuong, H.T., "Optimal Stochastic Path Control of Surface Ships in Shallow Water," The University of Michigan, Department of Naval Architecture and Marine Engineering, Report No. 189/ONR Report ONR-CR215-249-2F, 15 Aug. 1977. 2.- Parsons, M.G., and Cuong, H.T., "Adaptive Path Control of Surface Ships in Restricted Waters," The University of Michigan, Department of Naval Architecture and Marine Engineering, Report No. 211, March, 1980. 3. Cuong, H.T., "Investigation of Methods for Adaptive Path Control of Surface Ships," Ph.D. Dissertation, The University of Michigan, 1980. 4. Holley, W.E., and Bryson, A.E., Jr. "Multi-Input, Multi-Output Regulator Design for Constant Disturbances and Non-Zero Set Points with Application to Automatic Landing in a Crosswind," Stanford University Report SUDAAR No. 465, Aug., 1973. 5. Fujino, M., "Maneuverability in Restricted Waters: State of the Art," The University of Michigan, Department of Naval Architecture and Marine Engineering, Report No. 184, Aug., 1976. 6. Ogilvie, T.F., "Recent Progress Toward the Understanding and Prediction of Ship Motions," Proceedings of the 5th Symposium on Naval Hydrodynamics, Bergen, Norway, Sept. 10-12, 1964, pp. 3-80. 7. Fujino, M., "The Effect of Frequency Dependence of the Stability Derivatives on Maneuvering Motions," International Shipbuilding Progress, Vol. 22, No. 256, Dec., 1975, pp. 416-432. 8. Fujino, M., "Experimental Studies on Ship Manoeuvrability in Restricted Waters - Part I" International Shipbuilding Progress, Vol. 15, No. 168, Aug. 1968, pp. 279-301, and "-Part II," Vol. 17, No. 186, Feb. 1970, pp. 45-65. 9. Fujino, M., "Studies on Manoeuvrability of Ships in Restricted Waters," Selected Papers from the Journal of the Society of Naval Architects of Japan, Vol. 4, 1970, pp. 157-184. 10. Bryson, A.E., Jr., and Ho, Y.C., Applied Optimal Control, Blaisdell, Waltham, Mass, 1969. 11. Gelb, A. (ed), Applied Optimal Estimation, The MIT Press, Cambridge, Mass., 1974. 12. Parsons, M. G., and Greenblatt, J.E., "SHIPSIM/OPTSIM Simulation Program for Stationary, Linear Optimal Stochastic Control Systems," The University of Michigan, Department of Naval Architecture and Marine Engineering, Report No. 188/ONR Report ONR-CR215-249-1, 23 June 1977. 13. Astrom, K.J., and Kallstrom, C.G., "Identification of Ship-steering Dynamics," Automatica, Vol. 12, 1976, pp. 9-22. -63

-6414. Millers, H.F., "Modern Control Theory Applied to Ship Steering," Proceedings, of the IFAC/IFIP Symposium on Ship Operation Automation, Oslo, Norway, July, 1973. 15. Canner, W.H.P., "The Accuracy Requirements of Automatic Path Guidance," Proceedings of the Fourth Ship Control Systems Symposium, Den Helder, The Netherlands, Oct. 27-31, 1975, pp. 1-141 to 1-151. 16. Bystrom, L., and Kallstrom, C.G., "System Identification of Linear and Nonlinear Ship Steering Dynamics," Proceedings of the Fifth Ship Control Systems Symposium, Vol. 3, Annapolis, MD., Oct. 30-Nov. 3, 1978, pp. J2 2-1 to J2 2-21. 17. Newton, R.N., "Interaction Effects Between Ships Close Aboard in Deep Water," DTMB Report 1461, 1960. 18. Yung, T.-W., "Hydrodynamic Interactions Between Ships in Shallow Water," Ph.D. Dissertation, University of Michigan, Department of Naval Architecture and Marine Engineering, Jan., 1977. 19. Abkowitz, M.A., Ashe, G.M., and Fortson, R.M., "Interaction Effects of Ships Operating in Proximity in Deep and Shallow Water," Proceedings of the 11th Symposium on Naval Hydrodynamics, London, UK, March 28 - April 2, 1976, pp. VII.37-VII.57. 20. MacFarlane, A.G.J., "An Eigenvector Solution of the Optimal Linear Regulator Problem," Journal of Electronics and Control, Vol. 14, No. 6, June, 1963, pp. 643-654. 21. Potter, J.E., "Matrix Quadratic Solutions," SIAM Journal of Applied Mathematics, Vol. 14, No. 3, May, 1966, pp. 496-501. 22. Bryson, A.E., Jr., and Hall, W.E., Jr., "Optimal Control and Filter Synthesis by Eigenvector Decomposition," Stanford University Report SUDAAR No. 436, Nov., 1971. 23. Dorf, R.C., Modern Control Systems, Addison-Wesley Publishing Co., Reading, MA., 1974, pp. 111-117. 24. Chen, C.-T., Analysis and Synthesis of Linear Control Systems, Holt, Reinhart, and Winston, New York, 1975, pp. 155-157. 25. Parsons, M.G., "Optimal Control of Linear Systems," University of Michigan, Department of Naval Architecture and Marine Engineering, unpublished class notes, Sept., 1976. 26. private communication Capt. McSweeney, M/V MASABI MINER, April 1979. 27. Desoer, C.A., Notes for a Second Course on Linear Systems, van Nostrand Reinhold Co., New York, 1970, pp. 48-50.

UWWW OF~ MKXVIOA I1111511 I 0 lllll lll lllllll ll 3 9015 03095 0680