1. Dynamic positioning The dynamic positioning module is a control module with the following input and output: Input: Position measurement Wind measurement Anchor line force measurement Thrust measurement Output: Desired resultant forces and moment from the thrusters A PID controller and a Kalman filter-based controller exist. In addition it is possible to use an external dynamic positioning controller that can be coded by SITEF Ocean as a shared library (DLL), or use the internal SINTEF Ocean DP-system (MTDP) which is used for model tests in the SINTEF Ocean Basin. 1.1. PID controller The controller converts position and velocity errors into a demand for thrust forces to correct the errors. A decoupling approach allows PID control parameters to be specified separately for surge, sway and yaw. PID characteristics The control law of the PID (Proportional + Integral + Derivative) controller is \[\begin{array}{ll}F_{T0}(t)&=K_D\dot \varepsilon (t)+K_p\varepsilon (t)+K_I\displaystyle \int_{0}^{t}\!{\varepsilon (\tau)}\textrm{d}{\tau}\\\\\varepsilon (t)&=x_0-x(t)\\\\\dot \varepsilon (t)&=\dot x_0(t)-\dot x(t)\end{array}\] where \(x(t)\) filtered position \(x_0(t)\) desired (target) position \(\varepsilon (t)\) position error \(F_{T0}\) desired control force from thrusters \(K_D\) velocity feedback gain \(K_P\) position feedback gain \(K_I\) integral feedback gain Figure 1 shows an example of the PID controller’s frequency response from the position error \(\mathrm {\varepsilon }\) to control force \(\mathrm {F_{T0}}\) . It is seen that for low frequencies (below \(\mathrm {\omega _{I}}\) ) the controller acts like an integrator. In the medium frequency range it has proportional action, and for high frequencies (above \(\mathrm {\omega _{D}}\) ) the controller behaves like a differentiator. With increasing frequency, the phase angle of the controller goes from \(-90\mathrm {deg}\) to \(+90\mathrm {deg}\) . The relations between the frequencies \(\mathrm {\omega _{I}}\) and \(\mathrm {\omega _{V}}\) and the feedback gains are \[\begin{array}{ll}K_p&=K_P'(1+\displaystyle \frac{\omega _I}{\omega _D})\\\\K_I&=K_P'\omega _I\\\\K_D&=\displaystyle \frac{K_P'}{\omega _D}\end{array}\] where \(\mathrm {K_P'}\) is the level of the horizontal asymptote in the figure. Usually \(\mathrm {\omega _I\ll\omega _D}\) , so that \(\mathrm {K_P'\approx K_P}\) As \(\mathrm {\omega }\) approaches zero the controller gain approaches infinity due to the integral term. As a result, there will be no static (mean) position error - which is the purpose of integral action. Without the integral term (\(\mathrm {K_{I}=0}\) or equivalently, \(\mathrm {\omega _{I}=0}\) ) the static error would be \(x_{\mathrm {stat}}=F_{\mathrm {ext}}/K_{P}\) , where \(F_{\mathrm {ext}}\) is the static component of the external disturbing force. The purpose of the derivative action (i.e. the velocity feedback) is to provide phase lead. As is seen from Figure 1, the derivative term causes the phase to approach \(+90\,\mathrm {deg}\) . In general, a certain amount of phase lead will be required in order to obtain stability. Filtering It appears from Figure 1 that the PID’s gain also approaches infinity for high frequencies. This will make the system very sensitive to high frequency noise in the position and velocity signals as this noise would be amplified and put through to the thrusters. For this reason, high frequency components are removed from the signals by filtering before they are fed to the controller. There are two filters. One serves the double purpose of 1) low-pass filtering the position measurement to reduce the amount of measurement noise, and 2) differentiating the filtered position. The result is the input \(\mathrm {x}\) and and its derivative. Frequency response of low-pass filter (solid lines) and combined low-pass filter/differentiator (dashed) shows the frequency response of the low-pass-filter/differentiator, which has the cut-off frequency \(\mathrm {\omega _{C}}\) as the only parameter. The cut-off frequency is the frequency at which the attenuation starts. Above this frequency the attenuation is \(20\,\mathrm {dB}\) (10 times) per frequency decade. Figure 1. Frequency response of a PID controller (Dashed lines are asymptotes) Frequency response of low-pass filter (solid lines) and combined low-pass filter/differentiator (dashed) Figure 2. Frequency response of wave filter. In Frequency response of low-pass filter (solid lines) and combined low-pass filter/differentiator (dashed), the frequency axis is scaled as \(\mathrm {\omega /\omega _{c}}\) where \(\mathrm {\omega _{c}}\) is the cut-off frequency. The other filter is a wave filter. Its task is to remove the wave frequency components of the position measurement signal to prevent them from entering the thrusters. This is because thrusters usually have no chance of counteracting first-order wave forces. [tm-fig-freq-wave-filter] shows the frequency response of a wave filter for different values of the strength parameter \(\mathrm {\eta _{WF}}\). The wave filter is of band stop type. It is characterised by a frequency parameter and a strength parameter. The frequency response of the filter is given in [tm-fig-freq-wave-filter] for a number of values of the strength parameter \(\mathrm {\eta _{WF}}\) , to be chosen in the range from \(\mathrm {0}\) to \(\mathrm {1}\) . The maximum attenuation of the wave filter is \(16.5\mathrm {dB}\) or about \(\mathrm {7}\) times. It is intended that optimum filter performance should be achieved simply by equating the filter frequency with the peak frequency of the wave spectrum. As a consequence of the typical asymmetry of wave spectra, the frequency parameter is chosen to the left of the filter’s centre frequency. It is perfectly possible to filter off the wave frequency motion using the low-pass filter only. It is just a question of choosing the cut-off frequency low enough. The wave filter, however, is better suited for this task, since it gives more attenuation for less phase lag than the low-pass filter. The wave filter does not reduce any high frequency measurement noise. If such noise exists and is a problem, the cut-off frequency of the low-pass filter must be chosen low enough to give a sufficient suppression of it. In any circumstance the cut-off frequency should not be chosen any lower than necessary, in order to introduce as little negative phase as possible into control loop. On the other hand, the cut-off frequency should be chosen well below one half of the sampling rate of the measurements. Mechanical analogy In a DP system, the position and velocity feedback gains can be interpreted as stiffness and damping coefficients, respectively, and this mechanical analogy can be useful in determining these gains. However, due to the influence of the integral term and the filters, the analogy should not be pulled too far. Both the integral term and the filters may affect the system’s frequency response substantially, and the controlled system’s behaviour should therefore not be judged on the basis of \(\mathrm {K_P}\) and \(\mathrm {K_D}\) alone. A procedure for selecting the control parameters is included below. Cross-coupling effects The above description applies to each of the three controllers for surge, sway and yaw. When the vessel’s centre of mass is chosen as the point to be positioned and the mooring system (if any) is fairly symmetric, the effects of cross-coupling between the motions will in general be small. In this case it may be good enough to control the motions independently of each other. If, on the other hand, a point near the vessel’s bow is to be positioned, large coupling effects must be expected. In this case effective control may be difficult to achieve with independent controllers. There are three kinds of cross-coupling effects: through mass forces, damping forces and mooring forces. SIMO uses the vessel’s mass matrix (dry mass + hydrodynamic mass) to "couple together" the independently calculated force outputs from the three PIDs. This happens automatically without participation from the user. In this manner compensation for the effect of mass coupling, which is well defined, is achieved. It is expected that also the damping crosscoupling is partially compensated for in this manner, so that approximate decoupled control of any point on a non-moored ship can be obtained. For a moored ship however, the characteristics of the mooring system must be taken into account when the PID parameters are selected. Procedure for tuning the PID controller Tuning a control system is no exact science. The following procedure should serve as a guideline. For each of the three motions surge, sway and yaw, determine position and velocity feedback gains \(\mathrm {K_{P}}\) and \(\mathrm {K_{D}}\) that give the wanted stiffness and damping. Set the integral gains \(\mathrm {K_{I}}\) to zero. Turn off the wave frequency filter (\(\mathrm {\eta _{WF}=0}\) ) and set the cutoff frequency of the lowpass filter to \(\mathrm {1/5}\) of the sampling rate. From an initial position of non-equilibrium and with no dynamic environmental forces acting, simulate the transient motion towards equilibrium (To avoid complications caused by possible cross-coupling effects, it is safest to do this for one motion at a time, simulating pure transients in surge, sway and yaw). The transients should be non-oscillating. If they are not, the damping coefficients must be reduced. Make a simulation with dynamic environmental excitation. Inspect the thruster outputs to see if the amount of wave frequency components is admissible. If it is, proceed to 5. If not, select the appropriate value for the strength parameter of the wave filter according to [tm-fig-freq-wave-filter]. With the selected wave filter, simulate transients again. The wave filter introduces negative phase in the control loop, which reduces stability. If any transient is unstable or shows oscillatory behaviour that does not die out rapidly (i.e. during one or two half cycles), stability must be improved. If it is unacceptable to reduce the strength of the wave filter, the remedy is to: 1) make the controller phase more positive at the frequency of oscillation, and/or 2) reduce the controller gain at the same frequency. Remedy (1) is carried out by reducing \(\mathrm {K_{P}}\) while keeping \(\mathrm {K_{D}}\) constant. Try this choice first, for example by halving the value of \(\mathrm {K_{P}}\) . Remedy (2) is effectuated by reducing \(\mathrm {K_{P}}\) and \(\mathrm {K_{D}}\) by the same amount (e.g. halving them). Try this remedy if (1) fails or is insufficient. If integral action is desired, determine the upper frequency limit \(\mathrm {\omega _{I}}\) for this action, i.e. the frequency above which the integral action becomes insignificant (cf. [tm-fig-freq-wave-filter]) Then set \(\mathrm {K_{I}=w_{I}\cdot K_{P}}\) and simulate transients. The integral action also deteriorates stability. If stability is unsatisfactory, \(\mathrm {K_{I}}\) must be reduced. If it is known that the amount of measurement noise (the present version of SIMO does not simulate such noise) necessitates harder low-pass filtering of the position measurements, the cut-off frequency of the filter must be reduced from the value set according to paragraph 1 above. This may then lead to an adjustment of the controller gains as outlined in paragraph 4, until a satisfactory transient motion is obtained for all three motion variables. The above procedure reflects the standard knowledge from control theory that for a feedback system to be stable, the open-loop gain must be less than \(1(0\,\mathrm {dB})\) at the crossover frequency (i.e. the frequency at which the open-loop phase is \(-180\,\mathrm {deg}\) .) 1.2. Kalman filter The Kalman filter module was originally developed at SINTEF, division of Automatic Control during the period 1975-77, and has later been maintained and slightly modified by MARINTEK. The system is based on an optimal state estimator which is running in parallel with the physical process (the SIMO vessel). With feedback from estimated positions and velocities, the estimator acts as an adaptive controller. The mathematical model of the vessel is divided in three main sections: Model of wave-induced motions, high frequency (HF) model. Vessel motions in calm water, low frequency (LF) model. Model for slowly varying forces from waves and current. Slowly varying forces due to current and second order wave forces can be modelled by two methods: Estimation of global current velocity components Estimation of local bias forces In both cases, other unknown forces will also be estimated. The bias esimation represents the integrator effect of the controller. Input data to the Kalman filter routines are: \(x_\mathrm{ref}^G\) global x-reference \(y_\mathrm{ref}^G\) global y reference \(\psi_\mathrm{ref}\) heading reference \(x_\mathrm{loc}^L\) x-coordinate of point on vessel to be positioned \(y_\mathrm{loc}^L\) y-coordinate of point on vessel to be positioned \(d_x\) longitudinal quadratic drag coefficient \(d_y\) transverse quadratic drag coefficient \(d_\psi\) rotational quadratic drag coefficient \(v_\mathrm{typ}\) typical velocity \(m_x\) body mass and added mass in surge \(m_y\) body mass and added mass in sway \(m_\psi\) body inertia and added inertia in yaw \(R(3, \, 4)\) Kalman filter gain matrix \(G(3, \, 2)\) linear controller matrix Each time step, the following measurements are input to the Kalman filter: \(F_{Mx}^L\) measured longitudinal forces \(F_{My}^L\) measured transverse forces \(F_{M\psi}\) measured moment \(x_M^G\) measured global x-position of body origin \(y_M^G\) measured global y-position of body origin \(\psi_M\) measured heading The output from the routines are: \(F_{Tx}^L\) longitudinal thrust demand \(F_{My}^L\) transverse thrust demand \(F_{T\psi}\) thrust moment demand The Kalman filter has the following state variables: \(x_{LF}^G\) global LF x-position \(y_{LF}^G\) global LF y-position \(\psi_{LF}\) LF heading \(\dot x_{LF}^G\) global LF x-velocity \(\dot y_{LF}^G\) global LF y-velocity \(\dot \psi_{LF}\) LF heading rate \(v_{cx}^G\) North current velocity if current is estimated \(v_{cy}^G\) West current velocity if current is estimated \(F_{Bx}^L\) local bias x-force if bias is estimated \(F_{By}^L\) local bias y-force if bias is estimated \(F_{B\psi}\) bias moment \(x_{HF}^L\) local HF x-position (surge) \(y_{HF}^L\) local HF y-position (sway) \(\psi_{HF}\) HF heading (yaw) \(\dot x_{HF}^L\) local HF x-velocity \(\dot y_{HF}^L\) local HF y-velocity \(\dot \psi_{HF}\) HF heading rate \(\omega_x\) surge oscillation frequency \(\omega_y\) sway oscillation frequency \(\omega_\psi\) yaw oscillation frequency In addition, the following important help variables are used: \(\delta x^L\) position error surge \(\delta y^L\) position error sway \(\delta \psi\) heading error \(\varepsilon_x^L\) innovation in surge \(\varepsilon_y^L\) innovation in sway \(\varepsilon_\psi\) innovation in yaw The differences between the position measurements from the vessel and estimated measurements in the Kalman filter are called innovation processes. State variables in the estimator are updated by multiplications of innovations with Kalman filter gains. The first step is calculation of the innovations. \[\begin{array}{ll} \varepsilon _x^G&=x_M^G-x_{LF}^G\\\\ \varepsilon _y^G&=y_M^G-y_{LF}^G\\\\ \varepsilon _x^G, \varepsilon _y^G &\rightarrow \varepsilon _x^L,\varepsilon _y^L\\\\ \varepsilon _x^L&=\varepsilon _x^L-x_{HF}^L\\\\ \varepsilon _y^L&=\varepsilon _y^L-y_{HF}^L\\\\ \varepsilon _\psi&=\psi_M-\psi_{LF}-\psi_{HF} \end{array}\] Updating LF state variables: \[\begin{array}{ll}\delta x_{LF}^L&=R_{11}\varepsilon _x^L\\\\ \delta \dot x_{LF}^L&=R_{12}\varepsilon _x^L\\\\ \delta y_{LF}^L&=R_{21}\varepsilon _y^L\\\\ \delta \dot y_{LF}^L&=R_{22}\varepsilon _y^L\\\\ \delta x_{LF}^L,\delta y_{LF}^L&\rightarrow\delta x_{LF}^G,\delta y_{LF}^G\\\\ \delta \dot x_{LF}^L,\delta \dot y_{LF}^L&\rightarrow\delta \dot x_{LF}^G,\delta \dot y_{LF}^G\\\\ x_{LF}^G&=x_{LF}^G+\delta x_{LF}^G\\\\ \dot x_{LF}^G&=\dot x_{LF}^G+\delta \dot x_{LF}^G\\\\ y_{LF}^G&=y_{LF}^G+\delta y_{LF}^G\\\\ \dot y_{LF}^G&=\dot y_{LF}^G+\delta \dot y_{LF}^G\\\\ \psi_{LF}&=\psi_{LF}+R_{31}\varepsilon _\psi\\\\ \dot \psi_{LF}&=\dot \psi_{LF}+R_{32}\varepsilon _\psi\\\\ \end{array}\] Updating HF state variables: \[\begin{array}{ll}x_{HF}^L&=x_{HF}^L+R_{14}\varepsilon _x^L\\\\ y_{HF}^L&=y_{HF}^L+R_{24}\varepsilon _y^L\\\\ \psi_{HF}&=\psi_{HF}+R_{34}\varepsilon _\psi\\\\ \end{array}\] Updating bias estimates. If current is estimated: \[\begin{array}{ll}\delta v_{cx}^L&=R_{13}\varepsilon _x^L\\\\ \delta v_{cy}^L&=R_{23}\varepsilon _y^L\\\\ \delta v_{cx}^L,\delta v_{cy}^L&\rightarrow\delta v_{cx}^G,\delta v_{cy}^G\\\\ v_{cx}^G&=v_{cx}^G+\displaystyle \frac{\delta v_{cx}^G}{\mathrm {max}(|v_{cx}^G|,\,v_\mathrm {typ})}\\\\ v_{cy}^G&=v_{cy}^G+\displaystyle \frac{\delta v_{cy}^G}{\mathrm {max}(|v_{cy}^G|,\,v_\mathrm {typ})}\end{array}\] If bias force is estimated: \[\begin{array}{ll}F_{Bx}^L&=F_{Bx}^L+R_{13}\varepsilon _x^L\\\\ F_{By}^L&=F_{By}^L+R_{23}\varepsilon _y^L\\\\ \end{array}\] In both cases: \[F_{B\psi}=F_{B\psi}+R_{33}\varepsilon _\psi\] Then the frequencies of high frequency motion are updated according to (Balchen et al., 1980) Thrust demand is computed as feedback from the LF state variables and feed-forward from measured forces (wind) and estimated forces (wave drift and current). Feedback from LF state variables: \[\begin{array}{ll}x_\mathrm {loc}^L,y_\mathrm {loc}^L&\rightarrow x_\mathrm {loc}^G,y_\mathrm {loc}^G\\\\ \delta x^G&=x_{LF}^G+x_\mathrm {loc}^G-x_\mathrm {ref}^G\\\\ \delta y^G&=y_{LF}^G+y_\mathrm {loc}^G-y_\mathrm {ref}^G\\\\ \delta \psi&=\psi_{LF}-\psi_\mathrm {ref}\\\\ \delta x^G,\delta y^G&\rightarrow\delta x^L,\delta y^L\\\\ \dot x_{LF}^G,\\dot y_{LF}^G&\rightarrow\delta \dot x^L,\delta \dot y^L\\\\ F_{Tx}^L&=G_{11}\delta x^L+G_{12}\delta \dot x^L\\\\ F_{Ty}^L&=G_{21}\delta y^L+G_{22}\delta \dot y^L\\\\ F_{T\psi}&=G_{31}\delta \psi+G_{32}\psi_{LF}\\\\ \end{array}\] If current is estimated: \[\begin{array}{ll}v_{cx}^G,v_{cy}^G&\rightarrow v_{cx}^L,v_{cy}^L\\\\ F_{Dx}^L&=d_x\|v_{cx}^L|v_{cx}^L\\\\ F_{Dy}^L&=d_y\|v_{cy}^L|v_{cy}^L\end{array}\] If bias force is estimated: \[\begin{array}{ll}F_{Dx}^L=F_{Bx}^L\\\\ F_{Dy}^L=F_{By}^L\end{array}\] The calculated thrust demands are compensated for estimated and measured bias: \[\begin{array}{ll}F_{Tx}^L&=F_{Tx}^L-F_{Dx}^L-F_{Mx}^L\\\\ F_{Ty}^L&=F_{Ty}^L-F_{Dy}^L-F_{My}^L\\\\ F_{T\psi}&=F_{T\psi}-F_{B\psi}-F_{M\psi}\\\\ \end{array}\] The LF model consists of three discrete differential equations defined by the vessel’s mass and damping coefficients. If current is estimated: \[\begin{array}{ll}F_{Dx}^L&=-d_x|v_{cx}^L|v_{cx}^L\\\\ F_{Dy}^L&=-d_y|v_{cy}^L|v_{cy}^L\end{array}\] If bias force is estimated: \[F_{D\psi}=-d_\psi|\psi_{LF}|\psi_{LF}\] In both cases: \[\begin{array}{ll} F_{Dx}&=-2d_xv_\mathrm {typ}\dot x_{LF}^L+F_{Bx}^L\\\\ F_{Dy}&=-2d_yv_\mathrm {typ}\dot y_{LF}^L+F_{By}^L \end{array}\] The accelerations are expressed by: \[\begin{array}{ll}\ddot x_{LF}^L&=\displaystyle \frac{1}{m_x}(F_{Dx}^L+F_{Tx}^L+F_{Mx}^L)\\\\ \ddot y_{LF}^L&=\displaystyle \frac{1}{m_y}(F_{Dy}^L+F_{Ty}^L+F_{My}^L)\\\\ \ddot \psi_{LF}&=\displaystyle \frac{1}{m_\psi}(F_{D\psi}+F_{T\psi}+F_{M\psi})\end{array}\] The HF model for each degree of freedom is an adaptive harmonic oscillator. The outputs from the oscillators are estimates of wave-induced motion in the vessel coordinate system. \[\begin{array}{ll}\ddot x_{HF}^L&=-\omega _x^2x_{HF}^L\\\\ \ddot y_{HF}^L&=-\omega _y^2y_{HF}^L\\\\ \ddot \psi_{HF}&=-\omega _\psi^2\psi_{HF}^L\end{array}\] The accelerations are integrated to get a priori estimates of global LF-states and local HF-states. Initial tuning of the Kalman filter-based controller may be performed as described below. The controller gains: \[\begin{array}{ll}\displaystyle G(1,\,1)=-\frac{2\pi }{T_x}m_x\qquad\qquad&G(1,\,2)=-2\xi _x\sqrt{G(1,\,1)m_x}\\\\ \displaystyle G(2,\,1)=-\frac{2\pi }{T_y}m_y\qquad\qquad&G(2,\,2)=-2\xi _y\sqrt{G(2,\,1)m_y}\\\\ \displaystyle G(3,\,1)=-\frac{2\pi }{T_\psi}m_\psi\qquad\qquad&G(3,\,2)=-2\xi _\psi\sqrt{G(3,\,1)m_\psi}\end{array}\] where \(\mathrm {T_x}\), \(\mathrm {T_y}\) and \(\mathrm {T_\psi}\) are the chosen/wanted natural periods for surge, sway and yaw, respectively. The damping factors are given as \(\mathrm {\xi _x}\) , \(\mathrm {\xi _y}\) and \(\mathrm {\xi _\psi}\), and they are typically set to a value of 0.7 which indicates a damping level of 70% of the critical damping. Depending on the natural periods of the motions (low frequency) and on the wave frequency motions, cut-off frequencies \(\mathrm {\omega _c}\) for filtering the wave frequency motions are normally chosen. Typically, the cut-off periods for surge, sway and yaw of \(\mathrm {T_c=2\pi /\omega _c}\) are given values between \(25\mathrm {s}\) and \(50\mathrm {s}\) . Then Kalman filter coefficients (for initial tuning) may be chosen as: \[\begin{array}{ll}R(1,\,1)=2\xi _x\omega _{c,x}\qquad&R(1,\,2)=\omega _{c,x}^2\\\\ R(2,\,1)=2\xi _y\omega _{c,y}\qquad&R(2,\,2)=\omega _{c,y}^2\\\\ R(3,\,1)=2\xi _\psi\omega _{c,\psi}\qquad&R(3,\,2)=\omega _{c,\psi}^2\end{array}\] 1.3. Velocity PI controller When the simulator is run in manual mode, the user may use an input device, such as a joystick, to give a reference in velocity. The controller is similar to the controllers with position reference, but the velocity controller uses errors in velocity instead of errors in position. The controller converts velocity errors into a demand for thruster forces to correct the errors. A decoupling approach allows the PI control parameters to be specified separately for surge, sway and yaw, since these are dominant when the reference point is not given in the centre of gravity. PI characteristics The control law of the PI (Proportional + Integral) controller is \[\begin{array}{ll}F_{T0}&=K_P\varepsilon (t)+K_I\displaystyle \int_{0}^{t}\!{\varepsilon (\tau)}\textrm{d}{\tau}\\\\ \varepsilon (t)&=\dot x_0(t)-\dot x(t)\end{array}\] where \(\dot x(t)\) filtered velocity \(\dot x_0(t)\) desired (target) velocity \(\varepsilon(t)\) velocity error \(F_{T0}\) wanted control force from thrusters \(K_P\) velocity feedback gain \(K_I\) integral feedback gain The PI controller may be tuned similarly to the tuning of the PID controller in Section PID controller. The damping term \(\mathrm {K_D}\) is omitted because feedback in acceleration is very sensitive to measurement noice. The wave filter is the same as the one described in Section PID controller. Coupling forces Guidance