**2. Ship Drive Train and Its Mathematical Model**

The ship propulsion simulation model is based on a model scale ship called "Tito Neri", which is shown in Figure 1. A detailed picture of its azimuthing thrusters is shown in Figure 2, and its main particulars are given in Table 1. A schematic representation of one of its two drive trains is given in Figure 3. It consists of a DC motor that drives an azimuthing thruster with a ducted fixed pitch propeller. Although not shown in the figure, the upper shaft is supported by a shaft bearing.

**Figure 1.** Tito Neri overview.

**Figure 2.** Tito Neri azimuthing thrusters from astern.

**Figure 3.** Schematic representation of Tito Neri drive train.

**Table 1.** Main particulars of Tito Neri.


The system behavior is governed by two differential equations that interact with each other. One is related to the (faster) electrical circuit, and the other related to the (slower) mechanical part of the drive train. The differential equation commonly used to model an electric DC motor circuit is given by

$$L\_a \frac{di\_a}{dt} = \mathcal{U}\_a - K\_c \omega\_{cm} - R\_a i\_a \tag{1}$$

in which *La* is the inductance, *ia* the current, *Ua* the supply voltage, *Ke* the motor coefficient, *ωem* the motor speed (in rad/s), and *Ra* the resistance.

The reduction ratio between motor shaft and intermediate vertical shaft *igb*,12 between intermediate vertical shaft and propeller shaft *igb*,23 and the resulting total reduction ratio *igb*,13 are defined by

$$i\_{\mathbb{S}^b,12} = \frac{\omega\_{\text{em}}}{\omega\_{\text{int}}}, \quad i\_{\mathbb{S}^b,23} = \frac{\omega\_{\text{int}}}{\omega\_{\text{p}}}, \quad i\_{\mathbb{S}^b,13} = i\_{\mathbb{S}^b,12} i\_{\mathbb{S}^b,23}$$

The differential equation for electric motor speed, assuming constant friction torque on all three shafts, is given by

$$I\_p \frac{d\omega\_{em}}{dt} = M\_{b,cm} - M\_f - \frac{M\_p}{i\_{\mathcal{S}^{b,13}}} \tag{2}$$

in which brake motor torque *Mb*,*em* is given by

$$M\_{b,em} = K\_c i\_a \tag{3}$$

and in which the total polar moment of inertia is given by

$$I\_p = I\_{p,1} + \frac{I\_{p,2}}{i\_{\mathcal{g}^b,12}^2} + \frac{I\_{p,3}}{i\_{\mathcal{g}^b,13}^2}$$

and in which the total friction is given by:

$$M\_f = M\_{f,1} + \frac{M\_{f,2}}{i\_{\mathcal{S}^{b,12}}} + \frac{M\_{f,3}}{i\_{\mathcal{S}^{b,13}}}$$

The propeller torque *Mp* and thrust *T* are modeled following Carlton [28], making use of the torque and thrust coefficients *KQ* and *KT* at advance ratio *J* = 0:

$$M\_p = \frac{\mathcal{Q}}{\eta\_R} = \frac{K\_{\mathcal{Q},l=0} \rho \,\omega\_{cm}^2 D^5}{\eta\_R 4\pi^2 \,i\_{\mathcal{S}^b,13}^2} \tag{4}$$

in which *Q* is the open water propeller torque, *η<sup>R</sup>* is the relative rotative efficiency, *ρ* is the water density, and *D* is the propeller diameter. Although not further used in this paper, propeller thrust *T* and bollard pull force *FBP* are modeled by

$$T = \frac{K\_{T,l=0} \text{ } \rho \text{ } \omega\_{\text{cm}}^2 D^4}{4\pi^2 \hat{t}\_{\text{g}b,13}^2} \tag{5}$$

and

$$F\_{BP} = k\_p T(1 - t) \tag{6}$$

in which *kp* is the number of operating propellers, and (1 − *t*) corrects for thrust deduction.

To summarize, the following system of differential equations describes the nonlinear system dynamics:

$$\begin{cases} L\_d \frac{di\_d}{dt} = \mathcal{U}\_d - K\_\varepsilon \omega\_{\varepsilon m} - R\_d i\_d\\ I\_p \frac{d \omega\_{\varepsilon m}}{dt} = K\_\varepsilon i\_d - M\_f - \frac{M\_p}{i\_{\mathcal{G} \cdot 13}} \end{cases} \tag{7}$$

in which *Mp* is given by Equation (4). The unknown parameters of this model are inductance *La*, resistance *Ra*, motor coefficient *Ke*, polar moment of inertia *Ip*, friction torque *Mf* , propeller torque coefficient *KQ*,*J*<sup>=</sup>0, and relative rotative efficiency *ηR*. However, *KQ*,*J*=<sup>0</sup> and *η<sup>R</sup>* are observationally equivalent, meaning that (with the sensors in this experimental setup) they cannot be distinguished from each other. Therefore, propeller torque coefficient and relative rotative efficiency are combined into a single combined unknown parameter *KQ*,*J*=<sup>0</sup> *<sup>η</sup><sup>R</sup>* , leaving a total of six unknown parameters.

Note that the unknown parameters *KT*,*J*=<sup>0</sup> and (1 − *t*) are not further considered in this paper due to difficulties in measuring the small thrust force during the experiment.

#### *Linearized Propulsion System Model and Step-Responses*

In this section the ship propulsion system model (7) is linearized, and its analytical step responses are given. Later these will be shown to be useful tools for the identification of the unknown parameters.

The linearization process of the ship propulsion plant in free sailing mode is described in detail in [29,30], although in both papers no electric circuit including DC-motor was included. Note that in the main text of this paper only the main results are given, and details of the notation and the full derivations are available in Appendixes A–C. The normalized and linearized versions of (7) are given by

$$
\pi\_{\rm em} \frac{di\_a^\*}{dt} = \frac{\mathcal{U}\_{a,0}}{R\_a i\_{a,0}} \delta \mathcal{U}\_a^\* - \frac{K\_c \omega\_{\rm em,0}}{R\_d i\_{a,0}} \delta \omega^\* - \delta i\_a^\* \tag{8}
$$

$$
\pi\_{\omega} \frac{d\omega^\*}{dt} = \delta i\_a^\* - 2\eta\_{rm,0} \delta \omega^\* \tag{9}
$$

in which the delta-asterisk notation indicates normalized difference as follows:

$$
\delta \dot{a}\_a^\* = \frac{\delta \dot{a}\_a}{\dot{a}\_{a,0}} = \frac{\dot{i}\_a - \dot{i}\_{a,0}}{\dot{i}\_{a,0}}
$$

such that for example a value of *δU*∗ *<sup>a</sup>* = 0.05 means a +5% perturbation from the nominal value *Ua*,0. The two integration constants are defined as

$$\pi\_{\varepsilon m} = \frac{L\_a}{R\_a}, \quad \pi\_{\omega} = \frac{I\_p \omega\_{\varepsilon m,0}}{M\_{b,\varepsilon m,0}} = \frac{I\_p \omega\_{\varepsilon m,0}}{K\_a i\_{a,0}} \tag{10}$$

The transmission efficiency *ηtrm*,0 is related to the friction torque *Mf* by

$$
\eta\_{\rm rms,0} = \frac{M\_{b,cm,0} - M\_f}{M\_{b,cm,0}}
$$

When Equations (8) and (9) are put in state space notation, this results in the following linear system:

$$
\begin{bmatrix}
\frac{di\_d^\*}{dt} \\
\frac{d\omega^\*}{dt}
\end{bmatrix} = \begin{bmatrix}
\frac{1}{\tau\_{\omega}} & -\frac{2\eta\_{\ell m,0}}{\tau\_{\omega^\*}}
\end{bmatrix} \begin{bmatrix}
\delta i\_d^\* \\
\delta \omega^\*
\end{bmatrix} + \begin{bmatrix}
\frac{1}{\tau\_{\ell m}} \frac{lI\_{d,0}}{R\_d i\_{d,0}} \\
0
\end{bmatrix} \delta L I\_a^\*
\tag{11}
$$

The benefit of this notation is that it can easily be programmed and analyzed in software like MATLAB. Alternatively, the Laplace transfer function can be used. As derived in Appendix B, the two transfer functions from the supply voltage *δU*∗ *<sup>a</sup>* to the two state variables electric current *δi* ∗ *<sup>a</sup>* and rotation speed *δω*<sup>∗</sup> are

$$\frac{\delta i\_a^\*(s)}{\delta L l\_a^\*(s)} = \frac{(\pi\_{\omega,c}s + 1)\frac{\mathcal{U}\_{a,0}}{\mathcal{R}\_a l\_{a,0}}}{\pi\_{\rm cm}\tau\_{\omega,c}s^2 + (\pi\_{\rm cm} + \tau\_{\omega,c})s + 1 + \frac{1}{2\eta\_{\rm lrm}\eta\_{\rm l}}\frac{\mathcal{K}\_{\rm cl}\omega\_{\rm cen}\eta}{\mathcal{R}\_{\rm cl}\iota\_{\rm d}}}} \tag{12}$$

$$\frac{\delta\omega^{\*}\left(s\right)}{\delta L\_{a}^{\*}\left(s\right)} = \frac{\frac{1}{2\eta\_{\rm lrm}\eta\_{\rm l}}\frac{lI\_{a,0}}{R\_{a}\iota\_{a,0}}}{\tau\_{\rm cm}\tau\_{\omega,\rm c}s^{2} + \left(\tau\_{\rm cm} + \tau\_{\omega,\rm c}\right)s + 1 + \frac{1}{2\eta\_{\rm lrm}\eta\_{\rm c}}\frac{K\_{\rm c}\omega\_{\rm cm,0}}{R\_{a}\iota\_{a,0}}}\tag{13}$$

in which *τω*,*<sup>e</sup>* = *τω* <sup>2</sup>*ηtrm*,0 . The transfer function for current is recognizable as a summation of a bandpass system and a (lowpass) second-order system, while the transfer function for motor speed is (only) a second-order lowpass system.

Analytic expressions for the two poles *s*<sup>1</sup> and *s*2, the single zero *z*1, and the two DC-gains of the transfer functions are derived in Appendix B.

As derived in Appendix C the approximate response of motor speed to a unit step in supply voltage is given by

$$
\delta\omega^\*(t) \approx \mathcal{K}(1 - e^{s\_2 t})\tag{14}
$$

in which *K* = *Ua*,0 *Keω*0+2*ηtrm*,0*Raia*,0 . The response of current to a unit step in supply voltage is

$$\delta i\_a^\*(t) \approx K\_{LP} \left( 1 - e^{s\_2 t} \right) + K\_{BP} \left( \frac{1}{s\_2 - s\_1} e^{s\_2 t} - \frac{1}{s\_2 - s\_1} e^{s\_1 t} \right) \tag{15}$$

with *KLP* <sup>=</sup> *Ua*,02*ηtrm Keω*0+2*ηtrm*,0*Raia*,0 and *KBP* <sup>=</sup> *Ua*,0 *Laia*,0 .

#### **3. Applied Identification Techniques**

Once both the non-linear and the linearized plant models have been formulated, measurement data can be used to determine the unknown model parameters by making use of parameter identification techniques.

Many different identification techniques can be used, such as for instance the various possibilities that are included in the "system identification" toolbox of MATLAB. A possibility is to search for an optimal parameter vector *θ* by minimizing the (weighed) sum of squared errors given by the cost function *Vt*:

$$V\_t(\theta) = \frac{1}{N} \sum\_{t=1}^{N} \mathbf{e}(t, \theta)^T \mathbf{W}(\theta) \mathbf{e}(t, \theta) \tag{16}$$

where the error **e** is the difference between the vectors of measurement and simulation, containing all output signals that are to be taken into account:

$$\mathbf{e}(t) = \mathbf{y}(t)\_{measured} - \mathbf{y}(t)\_{model} \tag{17}$$

Another approach, which prevents usage of all time samples in the minimization algorithm and which ensures a balanced representation of system behavior throughout the frequency domain, is to define the cost function *Vω* directly in the frequency domain:

$$W\_{\omega}(\theta) = \frac{1}{N} \sum\_{\omega=1}^{N} \mathbf{e}(\omega, \theta)^{T} \mathbf{W}(\theta) \mathbf{e}(\omega, \theta) \tag{18}$$

in which the error is defined as the Euclidean norm of the error in the complex frequency response:

$$\mathbf{e} = ||\frac{Y(\omega)}{X(\omega)} - G(\theta, \omega)||\tag{19}$$

in which *<sup>Y</sup>*(*ω*) *<sup>X</sup>*(*ω*) indicates the measured frequency response data (FRD), model and *G*(*θ*, *ω*) indicates the modeled frequency response.

Within the two main groups "time domain approach" and "frequency domain approach", there are various possible refinements and alternatives. For an in-depth review, reference is made to Ljung [8].

The "goodness of fit" of a model with a given parameter set can be expressed in various ways. The quality metrics "FitPercent" and mean squared error "MSE" are used here:

$$FitPercent = 100 \left( 1 - \frac{||\mathbf{y}\_{measured} - \mathbf{y}\_{model}||}{||\mathbf{y}\_{measured} - \overline{\mathbf{y}\_{measured}}||} \right) \tag{20}$$

$$MSE = \frac{1}{N} \sum\_{t=1}^{N} \mathbf{e}^T(t)\mathbf{e}(t) \tag{21}$$

Equivalent versions of quality metrics can be defined for the goodness of fit in the frequency domain.

From the following non-exhaustive list of possible identification techniques, in this paper three different parameter identification procedures (1, 4, and 5) are applied to the "Tito Neri" drive train in the bollard pull condition:


Note that the frequency domain approaches 4 and 5 only differ in the way that they generate the experimental FRD. The subsequent parameter identification procedure is the same and aims to minimize cost function (18). 2 and 3 are not taken into account in the present work since they are investigated in open literature.

#### *3.1. Time Domain Identification: 1*

In this first method, for the sake of computational simplicity, the procedure to obtain parameters is split into two parts, assuming that the parameters do not change during the experimental time.

First, the stationary condition *<sup>d</sup>ω<sup>e</sup> dt* <sup>=</sup> 0 and *dia dt* = 0 is considered to reduce the number of unknowns, and a least-squares algorithm is applied. Starting from Equation (7) it is possible to obtain

$$\begin{cases} 0 = K\_c i\_a - M\_f - \frac{M\_p}{I\_{gb,13}}\\ 0 = \mathcal{U}\_a - K\_c \omega\_{cm} - R\_a i\_a \end{cases} \tag{22}$$

These equations are rearranged in matrix notation as follows, separating known from unknown variables:

$$
\begin{pmatrix} 1 & -i\_d & 0 & c \\ 0 & \omega\_{\varepsilon m} & i\_d & 0 \end{pmatrix} \begin{pmatrix} M\_f \\ K\_\varepsilon \\ R\_a \\ \frac{K\_{Q,j=0}}{\eta\_R} \end{pmatrix} = \begin{pmatrix} 0 \\ \mathcal{U}\_a \end{pmatrix} \tag{23}
$$

where *c* is obtained from Equation (4):

$$\mathcal{L} = \frac{\rho \omega\_{cm}^2 D^5}{4\pi^2 i\_{\mathcal{S}^{b,13}}^3} \tag{24}$$

In this way the system consists of two equations and four unknown variables (*Mf* , *Ke*, *Ra*, and *KQ*,*j*=<sup>0</sup> *<sup>η</sup><sup>R</sup>* ), such that <sup>∞</sup><sup>2</sup> solutions exist. However, if measurements at *<sup>n</sup>* different steady state operating points are available, *n* quadruplets have to satisfy the system of Equation (23), resulting in the following over-determined system:

$$
\begin{pmatrix}
1 & -i\_{a,1} & 0 & c\_1 \\
\cdots & \cdots & \cdots & \cdots \\
1 & -i\_{a,n} & 0 & c\_n \\
0 & \omega\_{cm,1} & i\_{a,1} & 0 \\
\cdots & \cdots & \cdots & \cdots \\
0 & \omega\_{cm,n} & i\_{a,n} & 0
\end{pmatrix}
\begin{pmatrix} M\_f \\ K\_\varepsilon \\ R\_\pi \\ \frac{K\_{Q,i=0}}{\eta\_R} \\ \eta\_R \end{pmatrix} = \begin{pmatrix} 0 \\ \mathcal{U}\_a \end{pmatrix} \tag{25}
$$

The last can be written in general form, as follows:

$$A\mathbf{x} = \mathbf{b} \tag{26}$$

The system (25) cannot be solved in principle since it is overdetermined. Although an exact solution does not exist, an approximate solution to (25) can be determined by means of, for instance, a (weighed) least-squares approach; in our case we used unweighted least squares. The final goal, according to notation reported in (26), is to evaluate the vector **x** that minimizes the squared *l* <sup>2</sup>*norm* of the residual, naming *A*, **x**, **b**, the coefficient matrix, the unknown vector, and the constant terms vector, respectively. The quantity *S*(**x**) to be minimized is written as follows, in matrix notation:

$$S(\mathbf{x}) = ||\mathbf{b} - A\mathbf{x}||^2\tag{27}$$

Differentiating the above equation, and setting to zero the solution, it is possible to obtain the Normal Equation:

$$A^T A \mathbf{x} = A^T(b) \tag{28}$$

If *ATA* is non-singular, the solution is given by solving the linear algebraic system (28).

Once *Mf* , *Ke*, *Ra*, and *KQ*,*j*=<sup>0</sup> *<sup>η</sup><sup>R</sup>* are known, the second part of the procedure is deterministic. By using Equation (10), it is possible to evaluate the inertia *Ip* and the motor inductance *La* :

$$\begin{aligned} L\_d &= \pi\_{\varepsilon m} \mathbb{R}\_d\\ I\_p &= \pi\_{\omega} \frac{K\_d i\_d}{\omega\_{\varepsilon m, 0}} \end{aligned} \tag{29}$$

To obtain the two parameters, knowledge of the time constants *τem* and *τω* is necessary, and the step response of both current and motor speed reported in Appendix C is used. From Equation (14) and fixing whichever time, *t* ∗ (authors suggest to use the *t* ∗ when the response is at 63.2%), since parameters are time independent, it is possible to obtain *s*2, as follows:

$$s\_2 \approx \frac{1}{t^\*} \ln\left(1 - \frac{\delta \omega^\*(t^\*)}{K}\right) \tag{30}$$

Substituting the value of *s*<sup>2</sup> into Equation (A34) and remembering the difference between *τω*,*<sup>e</sup>* and *τω* gives

$$
\pi\_{\omega} \approx \frac{-\mathbb{C}}{s\_2} \eta\_{trm,0} \tag{31}
$$

The evaluation of *τω*, asit can beintuitive from thelast equations, itis anapproximate solution.

The electric time-constant *τem* is more challenging to estimate. As reported in Appendix C the step response of current could be obtained as a summation of two terms. The first term is represented by a low-pass filter in its simplified form and the second by a bandpass filter as reported in Equation (15). The total response is known from the experiment, and all terms describing the low-pass filter are known at this stage; so, numerically, it is possible to obtain the shape of the bandpass filter response over time. A specific time called *t* ∗ should be fixed, and at that time the value of *δi* ∗ *<sup>a</sup>*,*BP*(*t* ∗) can be obtained. After some adjustment the following relation is obtained:

$$\delta i\_{d,BP}^\*(t^\*) - \frac{lI\_{d,0}}{\tau\_{\rm em}R\_di\_{d,0}} \left(\frac{1}{-\frac{c}{\tau\_{\omega\_d}} + \frac{1}{\tau\_{\rm em}}}\right) \left(e^{-\frac{\xi}{\hbar \omega\_d t^\*}} - e^{-\frac{1}{\xi\_{\rm em}}t^\*}\right) = 0\tag{32}$$

From the previous equation, it is not possible to obtain a solution in closed form for *τem*, and numerical methods must be used (i.e., bisection methods or Newton–Raphson method). Eventually, using Equation (29) *La* can be obtained.

#### *3.2. Frequency Domain Approach Using Sinusoidal Input Voltage Signals 4*

The idea of the this method is to generate a sinusoidal voltage of a specific frequency and amplitude, to superimpose it on a constant voltage value *Ua*,0, and to apply the resulting signal as a voltage input to the system, while recording the response of current *ia* and electric motor speed *ωem*. Based on the input and response at each frequency, the gain and phase of the transfer functions of the system are estimated with a correlation-based single-frequency approach [8,32], in line with Figure 4. Since this method does not deliver the unknown parameters of the model directly, it is called a non-parametric identification method. The non-parametric frequency response data (FRD) model can however be used as a basis to determine the parameters of the model.

**Figure 4.** Block diagram of experimental setup.

In more detail the basis of the method is to generate two harmonic signals:

$$x(t) = X \sin(\omega t)$$

and the out-of-phase signal

$$x'(t) = X\sin(\omega t + \pi/2) = X\cos(\omega t)$$

of which the first signal is used to excite the system. The response of the system is

$$y(t) = \mathcal{Y}\sin(\omega t + \phi) + n(t)$$

in which n(t) is a noise signal which is assumed uncorrelated with input and output signals. Both input signals *x*(*t*) and *x* (*t*), in combination with the output *y*(*t*), are used to

determine the cross-correlations and auto-correlation according to

$$R\_{xy} = \frac{1}{T} \int\_0^T X \sin(\omega t) \,\,\, Y \sin(\omega t + \varphi) dt + R\_{x\eta} = \frac{XY}{2} \cos\varphi + R\_{x\eta} \tag{33}$$

$$R\_{\mathbf{x'}y} = \frac{1}{T} \int\_0^T X \cos(\omega t) Y \sin(\omega t + \varphi) dt + R\_{\mathbf{x'}n} = \frac{XY}{2} \sin \varphi + R\_{\mathbf{x'}n} \tag{34}$$

$$R\_{\rm xx} = \frac{1}{T} \int\_0^T X \sin(\omega t) X \sin(\omega t) dt = \frac{X^2}{2} \tag{35}$$

where X is the amplitude of the input signal (in this case the amplitude of voltage *δU*∗ *a* ), and Y is the amplitude of the output signal under consideration (in this case the amplitude of motor current *δi* ∗ *<sup>a</sup>* or motor speed *δω*<sup>∗</sup> *em*). *Rxn* is the cross-correlation between input and noise, which reduces to zero with increasing measurement time. Division of Equation (33) by Equation (35) delivers the in phase (real) component of the frequency response:

$$\frac{R\_{xy}}{R\_{xx}} = \frac{Y}{X} \cos \varphi \tag{36}$$

while division of Equation (34) by Equation (35) gives the out-of-phase (imaginary) part of the response:

$$\frac{R\_{\mathbf{x'}\underline{y}}}{R\_{xx}} = \frac{Y}{X} \sin \varphi \tag{37}$$

Based on the real and imaginary components the gain and phase of the transfer function are calculated by

$$\frac{\chi}{X} = \sqrt{\left(\frac{R\_{xy}}{R\_{xx}}\right)^2 + \left(\frac{R\_{x'y}}{R\_{xx}}\right)^2} \tag{38}$$

$$\varphi = \arctan\left(\frac{R\_{x'y}}{R\_{xy}}\right) \tag{39}$$

By using this approach, the gain and phase can be determined experimentally for multiple appropriately spaced frequencies, resulting in a non-parametric frequency-response data (FRD) model. The results of the procedure are presented in Table 2. Subsequently, the procedure to derive the unknown system parameters from the obtained FRD model is based on minimization of the cost function (18).

The main advantage of the correlation method to determine an FRD model is the inherent high noise immunity.


**Table 2.** Experimental FRD based on single sinusoidal testing, followed by processing with the correlation approach.

#### *3.3. Noise Input Testing: 5*

An FRD model of a system can also be determined from the measured system response to a random input signal. This approach is often practical for processes that cannot be taken off-line for dedicated testing, but due to their nature do contain measurable random input disturbances. In this paper, a sequence of random supply voltage will be superimposed on the nominal supply voltage. The method is based on the relation between the transfer function *H*(*jω*), power spectral density of the input *Sxx*(*jω*), and cross-spectral density *Sxy*(*jω*) [8,31]:

$$H(j\omega) = \frac{S\_{xy}(j\omega)}{S\_{xx}(j\omega)}\tag{40}$$

The estimation of both the input power spectral density *Sxx* and the cross-spectral density *Sxy* requires sufficient length of data, and can be improved by application of suitable "windowing" and "smoothing", which can be done by averaging the spectrum derived from multiple segments of the total time-trace. Secondly, it is possible to increase the number of portions of a given time-trace by allowing a specific percentage of overlap between the parts.

For this method to work well, it is essential to ensure that the input signal is persistently exciting, which indicates that the signal power is sufficiently large for all frequencies of interest.

When using this method, the coherency *γ* usually is presented side by side with the estimated transfer function. It expresses the correlation between the input and output signal of the system with a value between 0 and 1, where 0 means no correlation and 1 means full correlation, thereby giving an idea of the quality of the estimated transfer function at different frequencies. Note that operations such as windowing, smoothing, and quantization of signals due to A-D conversion in the measurement system and noise in the measurement influence the coherency negatively.

#### **4. Experimental Campaign**

#### *4.1. Setup and Experimental Matrix*

The schematic experimental setup used is shown in Figure 4. The signal generator is operated via a customized graphical user interface and delivers the required voltage signal *Ua*,*set* to the amplifier, which in turn feeds the electric motor of the model scale ship with the supply voltage *Ua*. Two sensors are installed: a current sensor just before the electric motor and a 15 pulse encoder mounted on the motor shaft. The two sensor signals *ia* and *ωem*, together with the voltages *Ua*,*set* and *Ua*, are recorded with a data acquisition system. Although not discussed in detail in this paper, the transfer function of the amplifier itself could be determined experimentally, showing that the amplifier only causes a small drop in voltage (<1%), and a small phase lag (<2°) over the frequency range of interest.

Several experiments with varying sequences of voltage *Ua*,*set* have been done. The sampling rate of the data acquisition system was established based on the goals and duration of the specific experiment.

Trials were performed with the following input voltage signals: one staircase, nine sinusoidal waves with the different amplitudes and frequencies, a band-limited white noise input signal, and at the end a mix of the previous signals. Each identification technique uses data from a specific (set of) experiments. The final "mixed" test is used for validation purposes, as reported in Table 3.

**Table 3.** Experimental test matrix.


#### *4.2. Inspection of Current and Motor Speed Signals*

Initial measurements of the current revealed some unexpected behavior. The current signal showed a considerable amount of noise, and the reason was investigated. In particular, specific higher-order frequencies appeared when inspecting the FFT of the current signal. It was hypothesized that these higher-order frequencies, which are not captured by the linear or non-linear system model, could be caused by unmodeled system behavior. Examples could be, for instance, the gear-meshing frequency, shaft misalignment, unbalance in the shafting system, propeller blade passing frequencies, or cogging of the electric motor due to a discrete number of permanent magnets and the gaps in between them.

To obtain insight into the cause of the higher-order frequencies, an order tracking of current in the motor speed range from 220 to 1995 rpm was carried out, as shown in Figures 5 and 6. The figures reveal that although many harmonic frequencies were present in the current signal, the 6th and 12th harmonics of motor speed were particularly dominant. A similarly strong 6th and 12th harmonic were found when carrying out the test with disconnected gearwheels. Manual rotation of the motor shaft revealed a strong cogging effect at 6 times the motor shaft rate. Based on this it is concluded that the root cause of the higher-order frequencies lies in the interaction between rotor and stator of the electric motor.

Filtering has been considered to reduce the visually disturbing effect of coggingrelated harmonics from the plotted current signal. However, by filtering additional phase lag would be introduced, which would result in less steep current increase following a step in voltage, and it could reduce the amplitude of the current signal following a sinusoidal voltage input. In the end, it was decided to show the unfiltered current measurements.

**Figure 5.** Order waterfall plot.

**Figure 6.** FFT waterfall plot.

In addition, the motor speed signal showed unexpected behavior, which appeared to be caused by the sensor. A sketch of the encoder disk used in the experiments is shown in Figure 7. It is a round disk with 15 holes, which causes 15 pulses per revolution, generated by a photosensitive sensor. The motor speed is derived from the time interval between two upcoming flanks of the pulses. The resulting motor speed signal as shown in Figure 8 shows a repeating sequence of 15 motor speed values, indicating that the angle ΔΨ*i*,*<sup>j</sup>* between the holes varied slightly around 360/15 = 24°. No further corrections have been made to the signal, which explains the relatively "noisy" motor speed signal presented in the following sections.

**Figure 7.** Encoder disk.

**Figure 8.** Magnification of motor speed time history.

#### **5. Results and Discussion**

In this section, the results obtained with the different identification techniques are reported. Both the identification and validation analyses are carried out in both time and frequency domains.

In Table 4 the steady-state operating points recorded during the staircase experiment are reported. The time domain approach (1) uses all the five operating points to determine four out of the total of six unknown parameters. To determine the remaining two parameters *La* and *Ip* the transient response from operating point C to D is used.

The other identification approaches focus on operating point C. The reason to choose this point is that it corresponds to around 75% of the maximum supply voltage, which is a reasonable value thinking about the design of a full-scale propulsion plant.


**Table 4.** Evaluated operating points.

First, the resulting parameter sets of the different approaches are reported in Table 5 to appreciate the difference in terms of numerical value. The parameters derived from the spectral approach (5) are not reported as will be explained later. The table shows that the parameters obtained with the methods were of the same order of magnitude, but differences up to ≈100% were present. The effect of the different sets of parameters on the simulated system behavior is shown in the validation graphs.


**Table 5.** Identified parameters.

#### *5.1. Results Time Domain Analysis (1)*

The time domain identification method was used to derive the parameters from the staircase experiment. The supply voltage *Ua* during this experiment is shown in Figure 9, while the measured motor speed and current are shown in Figures 10 and 11.

Following the procedure outlined earlier, the five steady-state operating points during the staircase experiment were determined, and the parameters *Mf* , *Ke*, *Ra*, and *KQ*,*J*=0/*η<sup>R</sup>* were derived by the least-squares method. Subsequently, the transient response of motor speed and current, following the voltage step from C to D, was used to determine the parameters *Ip* and *La*.

The resulting set of parameters was implemented in the non-linear simulation model, and by using the staircase voltage signal as input, the model and its parameters are verified. The result is shown as the dashed red line in Figures 10 and 11. The motor speed matched the experimental data well: the stationary value errors were within 3% at all voltage levels. Close inspection of the transient responses showed that these were also captured

well. The simulated current signal had to be compared with very noisy experimental data, as explained earlier. Nevertheless, the static values seemed to be predicted well. Close inspection of the transient response shows that the simulation model could catch the timing and the initial steep slope of the current, but it was not able to represent the peak values in the current. It is concluded that this is either due to the limitations of the mathematical model, which might be too simple to describe the real physical phenomena, or due to the quality of the measured current signal.

**Figure 9.** Voltage time history of the staircase test.

**Figure 10.** Motor speed time history of the staircase test.

**Figure 11.** Current time history of the staircase test.

Figures 10 and 11 also show the results of the correlation approach, in continuous black line. However, since the staircase experiment was not used to determine the parameters using the correlation approach, this can be seen as validation of that method.

Figure 10 shows that the correlation approach predicted the motor speed behavior nearby the linearization point well, although the error between simulated and experimental data increased moving further away from the nominal operating point that was used in the sine experiments. Figure 11 shows that, compared to the time domain approach, the correlation approach was better able to predict the transient, although this method was also not able to catch the maximum current value.

To have an independent validation for the time-domain method, an experiment based on a mix of different input voltages, as shown in Figure 12, was used. Figure 13 shows that the parameter sets found by both methods led to similar dynamic behavior as the experiment, although a constant bias of around 50 rpm between simulated and sampled time histories was present. Figure 14 shows that the two parameter sets, in general, gave good correlation with the experiment, but both were unable to capture the maximum amplitudes of the current, which was particularly evident in the "sine wave" part from t = 10–13 s.

**Figure 12.** Voltage time history.

**Figure 13.** Motor speed time history.

**Figure 14.** Current time history.

#### *5.2. Results of Frequency Domain Analysis (4 and 5)*

The single frequency testing method was applied in the nominal operating point C that is defined in Table 4. The results of the nine experiments are plotted as asterisks data points in the Bode plots shown in Figures 15 and 16. The data points at 1000 and 5000 rad/s were discarded, as closer inspection of the time signals showed that the signal-to-noise ratio was too low to lead to meaningful results.

Based on the data points, the procedure as outlined above was followed, leading to the estimated parameters as listed in Table 5. The following values were iteratively determined from the experimental data points: *s*<sup>1</sup> = −2500 rad/s, *s*<sup>2</sup> = −9 rad/s, *z*<sup>1</sup> = −2 rad/s, *δω*∗ *δU*∗ *<sup>a</sup>* (*<sup>s</sup>* <sup>→</sup> <sup>0</sup>) = 1.24, *<sup>δ</sup><sup>i</sup>* ∗ *a δU*∗ *<sup>a</sup>* (*<sup>s</sup>* → <sup>0</sup>) = 0.71. Note that the locations of the poles and zero were read from the dB versions of the Bode plots.

**Figure 15.** Bode plot of *δω*<sup>∗</sup> *δU*∗ *a* .

**Figure 16.** Bode plot of *<sup>δ</sup><sup>i</sup>* ∗ *δU*∗ *a* .

To verify whether the estimation procedure was followed correctly, the found parameters were implemented in the transfer functions (12) and (13), which are plotted as solid lines in Figures 15 and 16. The agreement in trend and absolute numbers indicates that the procedure was followed correctly and that the linearized model can capture the reality well in the operating point under consideration. Validation in the time domain of the parameter set obtained with this method is reported in the previous section.

The shape of the transfer function for motor speed shows that up to 2 rad/s the response remained flat but then quickly dropped off due to the inertia of the drive train. The transfer function for current showed a flat response up to 1 rad/s and then started to rise due to the zero in the transfer function. Around 20 rad/s, it flattened out due to the inertia of the drive train. Somewhere after 1000 rad/s it dropped off, due to the electric pole, indicating that the current cannot follow the voltage variations anymore.

Figure 15 also shows the transfer function based on the parameter set derived with the time-domain approach. The response to low frequencies was good, but the drop in gain started slightly too early, which aligns with the low estimate of *Ip* in Table 5. In Figure 16 a substantial deviation from the data points is visible at frequencies higher than 10 rad/s, although the shape is clearly recognizable.

Finally, Figures 15 and 16 also show the results from approach 5. Between 0.4 and 400 rad/s the method resulted in a non-parametric frequency response that aligned well with the asterisk data points. In hindsight, the duration of the experiment should have been extended up to 1–2 min or even longer, instead of 30 s. A more extended trial would allow the estimation of the transfer function up to lower frequencies and would allow for further averaging over multiple data blocks to smooth out irregularities in the results. At frequencies above 400 rad/s, the signal-to-noise ratio dropped leading to bumps in the estimated frequency response.

Although the parameter estimation based on noise injection could be used to assess the unknown parameters from the frequency response, this is not performed here but is left for a further study on the potential of spectral methods for ship drive train identification.

#### **6. Future Outlook**

Application of the identification procedures on-board a real ship is expected to make the analytic derivations more complex because the system will, in that case, include other/additional components such as, for instance, a diesel engine and engine speed governor. This introduces at least one extra state equation due to the integral term in the PI(D) governor. An additional state, due the longitudinal equation of motion, and additional parameters would be added if the approach would be extended to free sailing instead of bollard pull conditions. The effect of such additions on the ability to determine parameters needs to be investigated in the future. On the positive side, it has to be noted that in reality it is not likely that all ship drive-train parameters are unknown, which helps to determine estimates of other unknown or more uncertain parameters. More work is required to investigate what parameter estimation procedure would be required for a real ship and drive train.

In the authors' opinion, in the future the presented algorithms could potentially be part of a condition-based maintenance system. By monitoring parameter variations of a propulsion drive train in real time, it could be possible to detect the degradation (or malfunction) of the machinery, and perhaps even to identify the root cause. For instance, an increase in friction coefficient *Mf* could mean wear in the bearings, an increase in *KQ* could mean that the propeller needs to be cleaned, etc. Another possible use of the presented techniques is to assess the correspondence between the design values with the real one, in fact, during the shipbuilding progress some change, or unexpected modification, can modify the original design values.

#### **7. Conclusions and Recommendations**

In this paper different parameter identification techniques were discussed and applied to experimentally determine the unknown parameters of a model scale ship drive train in bollard pull conditions.

A set of dedicated experiments was conducted using different DC voltage signals. In all tests the current was affected by a strong noise due to motor cogging. It is therefore recommended to use an electric motor with less strong cogging effect for future experiments. Moreover, the 15 holes encoder was found to give a low-quality motor speed measurement and should be improved.

Three different approaches to determine the unknown DC-electric propulsion plant parameters are discussed including their merits and weaknesses. For now, all three approaches remain candidates to be part of a (real-time) full-scale parameter identification system, which is one of the primary goals.

Two obtained parameter sets have been implemented in a simulation model, and the results were validated against independent measurements, both in the frequency and in the time domains. The time domain results obtained by implementing both parameter sets in the model compared well against the measurements, although differences were present.

In order to move towards firm conclusions about the value of the applied parameter estimation methods for condition monitoring, it is recommended to consider the sensitivity and uncertainty related to the approaches. This recommendation is supported by the relatively large differences between the parameter sets as determined in this paper, and the relatively small differences in time domain response.

**Author Contributions:** Conceptualization, A.V. and M.M.; methodology, A.V. and M.M.; software, A.V. and M.M.; writing—original draft preparation, A.V. and M.M.; writing—review and editing, A.V. and M.M.; experimental testing A.V. and M.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partly supported by Maritime by Holland (in Dutch: NML) and the Ministry of Economic Affairs of the Netherlands. Furthermore, this research was partly funded by University of Genoa within the program "Incentive for research periods abroad".

**Acknowledgments:** The authors want to acknowledge Eng. Vittorio Garofano of Delft University of Technology for his essential support during the experimental campaign.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**



#### **Subscripts and Superscripts**



#### **Appendix A. Normalisation and Linearisation**

Assume a variable that is the product of powers of other variables:

$$Z = c \,\, Y^{\varepsilon} X \tag{A1}$$

where *c* is a constant multiplier and *e* is a constant exponent. In an equilibrium point the variable *Z* equals

$$Z\_0 = c \, Y\_0^\varepsilon X\_0 \tag{A2}$$

Normalisation of Equation (A1) by Equation (A2) results in

$$\frac{Z}{Z\_0} = \left(\frac{\chi}{Y\_0}\right)^c \frac{X}{X\_0} \tag{A3}$$

If, by definition,

$$X^\* = \frac{X}{X\_0}, Y^\* = \frac{Y}{Y\_0}, Z^\* = \frac{Z}{Z\_0} \tag{A4}$$

then

$$Z^\* = Y^{\*\varepsilon} X^\* \tag{A5}$$

Now that the constant value *c* has been removed by the normalization, the next step is to remove the non-linearity from Equation (A5). Differentiation of Equation (A3) by using the chain rule gives

$$\frac{dZ}{Z\_0} = \left(\frac{\chi}{Y\_0}\right)^c \frac{dX}{X\_0} + \varepsilon \left(\frac{\chi}{Y\_0}\right)^{c-1} \frac{X}{X\_0} \frac{dY}{Y\_0} \tag{A6}$$

Near equilibrium *dX*, *dY* and *dZ* become small increments *δX*, *δY* and *δZ*. Division of *X* = *X*<sup>0</sup> + *δX* by *X*<sup>0</sup> delivers *<sup>X</sup> <sup>X</sup>*<sup>0</sup> <sup>=</sup> <sup>1</sup> <sup>+</sup> *<sup>δ</sup><sup>X</sup> <sup>X</sup>*<sup>0</sup> and likewise *<sup>Y</sup> <sup>Y</sup>*<sup>0</sup> <sup>=</sup> <sup>1</sup> <sup>+</sup> *<sup>δ</sup><sup>Y</sup> <sup>Y</sup>*<sup>0</sup> . Substitution of this in Equation (A6) gives

$$\frac{\delta Z}{Z\_0} = \left(1 + \frac{\delta Y}{Y\_0}\right)^{\varepsilon} \frac{\delta X}{X\_0}$$

$$1 + e\left(1 + \frac{\delta Y}{Y\_0}\right)^{\varepsilon - 1} \left(1 + \frac{\delta X}{X\_0}\right) \frac{\delta Y}{Y\_0} \tag{A7}$$

Taylor series expansion of Equation (A7) and neglecting the second and higher order terms leaves

$$\frac{\delta Z}{Z\_0} = \frac{\delta X}{X\_0} + \varepsilon \frac{\delta Y}{Y\_0} \tag{A8}$$

which by introduction of the shorthand notation for differential increment:

$$
\delta Z^\* = \frac{\delta Z}{Z\_0} = \frac{Z}{Z\_0} - 1 \tag{A9}
$$

this can be written as

$$
\delta Z^\* = \delta X^\* + e \,\delta Y^\* \tag{A10}
$$

The latter equation relates the relative change in output *Z* to the relative change in inputs *X* and *Y*, where the constant *e*, which was present as an exponent in the original Equation (A2), has changed to a constant multiplication factor. Secondly the multiplication of *X* and *Y* has turned into a summation. For further background on the linearization process, reference is made to Dorf and Bishop [33] and Franklin et al. [34].

The demonstrated concepts of normalization and linearization are the basis for the following section where they will be applied in the linearization of the system model.

#### **Appendix B. Derivation of Linearized System Model**

The electrical circuit of the DC motor is modeled by

$$L\_a \frac{di\_a}{dt} = \mathcal{U}\_a - K\_c \omega\_{cm} - R\_a i\_a \tag{A11}$$

All three right hand side terms vary around equilibrium:

$$\mathcal{U}\_a = \mathcal{U}\_{a,0} + \delta \mathcal{U}\_{a\prime} \quad \omega\_{cm} = \omega\_{cm,0} + \delta \omega\_{cm\prime} \quad \dot{\imath}\_a = \dot{\imath}\_{a,0} + \delta \dot{\imath}\_a$$

In static conditions the right hand side of Equation (A11) equals zero:

$$0 = \mathcal{U}\_{a,0} - K\_c \omega\_{cm,0} - R\_a i\_{a,0} \tag{A12}$$

This means that only the small increments are of importance:

$$L\_a \frac{di\_a}{dt} = \delta \mathcal{U}\_a^\* - K\_c \delta \omega\_{em}^\* - R\_a \delta i\_a^\* \tag{A13}$$

Division of all terms of Equation (A13) by nominal supply voltage minus the nominal emf (*Ua*,0 − *Keωem*,0) or alternatively by its equivalent *Raia*,0 gives

$$\frac{L\_{\rm a}}{R\_{\rm a}i\_{a,0}} \frac{d\dot{i}\_{a}}{dt} = \frac{1}{R\_{\rm a}i\_{a,0}} \frac{lI\_{a,0}}{lI\_{a,0}} \delta lI\_{a} - \frac{K\_{\rm c}}{R\_{\rm a}i\_{a,0}} \frac{\omega\_{cm,0}}{\omega\_{cm,0}} \delta \omega\_{cm} - \frac{R\_{\rm a}}{R\_{\rm a}i\_{a,0}} \delta \dot{i}\_{a} \tag{A14}$$

This can be shortened to

$$
\pi\_{em}\frac{di\_a^\*}{dt} = \frac{\mathcal{U}\_{a,0}}{R\_a i\_{a,0}} \delta \mathcal{U}\_a^\* - \frac{K\_\epsilon \omega\_{\epsilon m,0}}{R\_a i\_{a,0}} \delta \omega^\* - \delta i\_a^\* \tag{A15}
$$

in which the subscript *em* is intentionally dropped from *δω*∗ *em* because *δω*<sup>∗</sup> *em* = *δω*<sup>∗</sup> *p* and where

$$
\pi\_{\rm em} = \frac{L\_a}{R\_a} \tag{A16}
$$

The shaft dynamics including constant friction term are described by

$$I\_p \frac{d\omega\_{\rm em}}{dt} = M\_{b\,\rm arm} - M\_f - \frac{M\_p}{\dot{t}\_{\rm \%} \mu\_{\rm 13}} \tag{A17}$$

in which shaft inertia is assumed constant implying that change of mass of water, entrained by the propeller, is neglected. The brake motor torque is related to current by

$$M\_{b, \varepsilon m} = K\_{\varepsilon} \, i\_a \tag{A18}$$

The non-constant torque terms of Equation (A17) vary around equilibrium:

$$M\_{b,em} = M\_{b,em,0} + \delta M\_{b,em} = K\_e(i\_{a,0} + \delta i\_a)$$

and

$$M\_p = M\_{p,0} + \delta M\_p$$

such that:

$$I\_p \frac{d\omega\_{\rm cm}}{dt} = K\_\varepsilon (i\_{\rm a,0} + \delta i\_\rm a) - M\_f - \frac{M\_{p,0}}{i\_{\rm gb,13}} - \frac{\delta M\_p}{i\_{\rm gb,13}} \tag{A19}$$

In steady nominal condition the driving torque and the load-torque are equal:

$$0 = K\_{\mathfrak{e}} i\_{\mathfrak{a},0} - M\_f - \frac{M\_{\mathfrak{p},0}}{i\_{\mathfrak{g}b,13}} \tag{A20}$$

Subtracting Equation (A20) from Equation (A19) shows that only the small increments are of importance:

$$I\_p \frac{d\omega\_{cm}}{dt} = K\_c \delta i\_a - \frac{\delta M\_p}{i\_{\S^{b,13}}}$$

Normalizing all terms with nominal motor torque gives

$$\frac{I\_p \omega\_{\varepsilon m,0}}{K\_\epsilon i\_{a,0}} \frac{d\omega^\*}{dt} = \delta i\_a^\* - \frac{M\_{p,0}}{M\_{b,\varepsilon m,0} i\_{\overline{\varepsilon} b, 13}} \delta M\_p^\* \tag{A21}$$

in which the subscript *em* is intentionally dropped. The integration constant is defined as

$$
\pi\_{\omega} = \frac{I\_p \omega\_{cm,0}}{M\_{b,cm,0}} = \frac{I\_p \omega\_{cm,0}}{K\_c \ i\_{a,0}} \tag{A22}
$$

After noting that the multiplier in the second term of the right hand side of Equation (A21) can be written as

$$\frac{M\_{p,0}}{i\_{\rm gb,13}M\_{b,cm,0}} = \eta\_{rm,0} \tag{A23}$$

and implementing

$$
\delta \mathcal{M}\_p^\* = 2 \delta \omega^\*
$$

the normalised linearised differential equation for shaft rotation is given by

$$
\pi\_{\omega} \frac{d\omega^\*}{dt} = \delta i\_a^\* - 2\eta\_{trm,0} \delta \omega^\* \tag{A24}
$$

Introduction of the Laplace operator into Equation (A24) and re-arranging gives

$$\left(\frac{\tau\_{\omega}}{2\eta\_{rm,0}}s+1\right)\delta\omega^\* = \frac{1}{2\eta\_{rm,0}}\delta i\_a^\* \tag{A25}$$

which can be shortened by introduction of the effective time-constant:

$$
\pi\_{\omega,\varepsilon} = \frac{\pi\_{\omega}}{2\eta\_{\text{trunc},0}} \tag{A26}
$$

such that

$$(\pi\_{\omega,a}\mathbf{s}+1)\delta\omega^\* = \frac{1}{2\eta\_{\rm trm,0}}\delta i\_a^\* \tag{A27}$$

Similarly, introduction of the Laplace operator in the differential equation for current Equation (A15) and reordering gives

$$
\delta i\_a^\* = \frac{\left(\frac{\mathcal{U}\_{a,0}}{R\_{a^\* \!\!I\_{a,0}}}\right)}{\left(\mathcal{I}\_{cm}s + 1\right)} \delta \mathcal{U}\_a^\* - \frac{\left(\frac{K\_{\varepsilon}\omega\_{cm,0}}{R\_{\varepsilon}\mathfrak{r}\_{\mathfrak{a},0}}\right)}{\left(\mathcal{I}\_{cm}s + 1\right)} \delta \omega^\* \tag{A.28}
$$

Substitution of Equation (A28) into Equation (A27) and reordering gives the transfer function from supply voltage to rotation speed:

$$\frac{\delta\omega^{\*}}{\delta\mathcal{U}\_{a}^{\*}} = \frac{\frac{1}{2\eta\_{\rm lrm}\eta}\frac{\mathcal{U}\_{a,0}}{R\_{a}\mathfrak{i}\_{a,0}}}{\tau\_{\rm cm}\tau\_{\omega,\rm c}\mathfrak{s}^{2} + (\tau\_{\rm cm} + \tau\_{\omega,\rm c})\mathfrak{s} + 1 + \frac{1}{2\eta\_{\rm lrm}\eta}\frac{\mathcal{K}\omega\_{\rm cm,0}}{R\_{a}\mathfrak{i}\_{a,0}}}} \tag{A29}$$

In a similar way substitution of Equation (A27) into Equation (A28) and reordering gives the transfer function from supply voltage to current:

$$\frac{\delta i\_a^\*}{\delta \mathcal{U}\_a^\*} = \frac{(\tau\_{\omega,\mathcal{E}}s + 1)\frac{\mathcal{U}\_{a,0}}{\mathcal{R}\_a \mathcal{U}\_{a,0}}}{\tau\_{\varepsilon m}\tau\_{\omega,\mathcal{E}}s^2 + (\tau\_{\varepsilon m} + \tau\_{\omega,\mathcal{E}})s + 1 + \frac{1}{2\eta\_{\varepsilon m\eta,0}}\frac{\mathcal{K}\_{\varepsilon}\omega\_{\varepsilon m,0}}{\mathcal{R}\_a \mathcal{I}\_{a,0}}} \tag{A.30}$$

The characteristic equation of the two transfer functions Equations (A29) and (A30) is given by

$$
\pi\_{\rm \varepsilon m} \pi\_{\omega, \varepsilon} \mathrm{s}^2 + (\pi\_{\varepsilon m} + \pi\_{\omega, \varepsilon}) \mathrm{s} + 1 + \frac{1}{2 \eta\_{\mathrm{trm},0}} \frac{K\_{\mathrm{e}} \omega\_{\mathrm{em},0}}{R\_{\mathrm{a}} i\_{a,0}} \tag{A31}
$$

If we define

$$C = 1 + \frac{1}{2\eta\_{rm,0}} \frac{K\_t \omega\_{cm,0}}{R\_a i\_{a,0}} \tag{A32}$$

and

$$
\zeta = \frac{\tau\_{c\pi}}{\tau\_{\omega\mu}}
$$

then the characteristic equation can be written as

$$
\zeta x\_{\omega,\varepsilon}s^2 + (1+\zeta)s + \frac{C}{x\_{\omega,\varepsilon}}
$$

The two exact roots of Equation (A31) can now be determined by the ABC formula:

$$s\_{12} = \frac{-\left(1+\zeta\right) \pm \sqrt{\left(1+\zeta\right)^2 - 4C\zeta}}{2\zeta x\_{\omega,c}}$$

which can be written as

$$s\_{12} = \frac{-\left(1+\zeta\right) \pm \left(1+\zeta\right)\sqrt{1-\frac{4C\zeta}{\left(1+\zeta\right)^2}}}{2\zeta x\_{\omega,\varepsilon}}$$

The electrical time constant is much smaller than the effective time constant for the shaft; therefore, *ζ* 1. Application of Taylor expansion for the square root operation and leaving out second order terms gives

$$s\_{12} \approx \frac{-(1+\zeta) \pm (1+\zeta)\left(1 - 2C\zeta \frac{1}{(1+\zeta)^2} \dots \right)}{2\zeta x\_{\omega,x}}$$

Another Taylor expansion for the inverse square operation gives

$$s\_{12} \approx \frac{-(1+\zeta) \pm (1+\zeta)(1-2C\zeta(1-2\zeta...)\dots)}{2\zeta x\_{\omega,\mathfrak{c}}}$$

Further simplification gives the two approximate poles as

$$s\_1 \approx \frac{-1}{\zeta x\_{\omega, \varepsilon}} = \frac{-1}{x\_{\varepsilon m}} \tag{A33}$$

and

$$s\_2 \approx \frac{-\mathcal{C}}{x\_{\omega,\varepsilon}}\tag{A34}$$

Besides the two system poles, transfer function (A30) has a single zero which lies at

$$z\_1 = \frac{-1}{x\_{\omega,c}}\tag{A35}$$

The DC-gain of transfer function (A29) is given by

$$\frac{\delta\omega^{\*}}{\delta l I\_{a}^{\*}}(s \to 0) = \frac{l I\_{a,0}}{2\eta\_{rm,0}R\_{a}i\_{a,0} + K\_{\varepsilon}\omega\_{cm,0}}\tag{A36}$$

The DC-gain of transfer function (A30) is given by

$$\frac{\delta i\_a^\*}{\delta l I\_a^\*} (s \to 0) = \frac{2\eta\_{rm,0} l I\_{a,0}}{2\eta\_{rm,0} \mathcal{R}\_a i\_{a,0} + \mathcal{K}\_c \omega\_{cm,0}} \tag{A37}$$

#### **Appendix C. Step Response of Motor Speed and Current**

The exact response of motor speed to a unit step in voltage is given by

$$
\delta\omega^\*(t) = K \left( 1 + \frac{s\_1}{s\_2 - s\_1} e^{s\_2 t} - \frac{s\_2}{s\_2 - s\_1} e^{s\_1 t} \right) \tag{A.38}
$$

in which *K* = *Ua*,0 *Keω*0+2*ηtrm*,0*Raia*,0 . However, because |*s*1| >> |*s*2|, the step-response can be approximated by a first order system response:

$$
\delta\omega^\*(t) \approx K(1 - e^{s\_2 t})\tag{A39}
$$

The derivation of the step response of current starts with Equation (A30), which can be written as the summation of an overdamped second order lowpass (LP) system and a second order bandpass (BP) system:

$$G(\mathbf{s}) = G\_{LP}(\mathbf{s}) + G\_{BP}(\mathbf{s})\tag{A40}$$

The step response of the lowpass system is given by

$$
\delta i\_{a,LP}^\*(t) = K\_{LP} \left( 1 + \frac{s\_1}{s\_2 - s\_1} e^{s\_2 t} - \frac{s\_2}{s\_2 - s\_1} e^{s\_1 t} \right) \tag{A41}
$$

with *KLP* <sup>=</sup> *Ua*,02*ηtrm Keω*0+2*ηtrm*,0*Raia*,0 . Again, because |*s*1| >> |*s*2|, this can be approximated by a first order system response:

$$\delta \dot{a}\_{a,LP}^\*(t) \approx K\_{LP} \left(1 - e^{s\_2 t}\right) \tag{A42}$$

The step response of the bandpass part of the system is given by

$$
\delta i\_{a,BP}^\*(t) = K\_{BP} \left( \frac{1}{s\_2 - s\_1} e^{s\_2 t} - \frac{1}{s\_2 - s\_1} e^{s\_1 t} \right) \tag{A43}
$$

where *KBP* <sup>=</sup> *Ua*,0 *Laia*,0 . The total response of current to a unit step in voltage is the sum of Equations (A41) and (A43):

$$\delta i\_a^\*(t) = K\_{LP} \left( 1 + \frac{s\_1}{s\_2 - s\_1} e^{s\_2 t} - \frac{s\_2}{s\_2 - s\_1} e^{s\_1 t} \right) + K\_{BP} \left( \frac{1}{s\_2 - s\_1} e^{s\_2 t} - \frac{1}{s\_2 - s\_1} e^{s\_1 t} \right) \tag{A44}$$

or including the simplification:

$$
\delta i\_d^\*(t) \approx K\_{LP} \left( 1 - e^{s\_2 t} \right) + K\_{BP} \left( \frac{1}{s\_2 - s\_1} e^{s\_2 t} - \frac{1}{s\_2 - s\_1} e^{s\_1 t} \right) \tag{A45}
$$

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland www.mdpi.com

*Journal of Marine Science and Engineering* Editorial Office E-mail: jmse@mdpi.com www.mdpi.com/journal/jmse

Disclaimer/Publisher's Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Academic Open Access Publishing

mdpi.com ISBN 978-3-0365-9357-9