Model Structures and Identiﬁcation for Fully Embedded Thrusters: 360-Degrees-Steerable Steering-Grid and Four-Channel Thrusters

: The European Watertruck + project introduced a new ﬂeet of self-propelled inland cargo barges to the European waters, in order to induce more sustainable freight transport in the European hinterland. An augmentation of the automation level of this ﬂeet could further advance their competitiveness and potentially pave the way for unmanned inland cargo vessels. The motion control of such a vessel forms a key component in this envisaged automation chain and beneﬁts from the knowledge of the capabilities of the propulsion system, which here envelops a 360-degrees-steerable steering-grid thruster in conjunction with a 360-degrees-steerable four-channel thruster. Therefore, this study details the mechanical design of both thrusters and lists their experimental towing-tank data. Furthermore, two different modelling methods are offered, one theoretically based and one using a multilayer neural network. A model structure comparison, based on a bias-variance trade-off, veriﬁes the adequacy of the theoretical model which got expended with an angle-dependent thrust deduction coefﬁcient. In addition, several multilayer feedforward neural network architectures exemplify their inherent capability to model the complex, nonlinear, ﬂow phenomena inside the thrusters. These identiﬁed model structures can additionally improve thrust allocation algorithms and offer better plant models to study more advanced control strategies.


Introduction
On average, current inland cargo vessels generate lower external cargo transportation costs compared to other modes of freight transport, where the authors in [1][2][3] define these costs as accidents, air pollution, climate, noise and congestion. Nevertheless, road-based cargo transport presently dominates the European hinterland freight transport [4]. Therefore, the European Commission targets inducing a cargo transport modal shift from road to rail and waterborne transport [5]. Subsequently, the European Watertruck + [6] project aims to facilitate this modal shift by constructing modular push vessels and barges which are either passive or self-propelled, in order to decouple sailing and transshipment time. The authors believe that an increase in the automation levels of these Watertruck + vessels, and inland cargo vessels in general, could further help the envisaged cargo shift. Furthermore, they believe that this automation augmentation could potentially pave the way for future unmanned or autonomous inland cargo vessels.
Therefore, in [7,8], they constructed a scale model of such a European class type I [9] self-propelled inland cargo barge. These self-propelled barges navigate with a non-conventional, over-actuated, fully-embedded actuation system configuration consisting of a 360-degrees-steerable lsteering-grid thruster in the bow in conjunction with a 360-degrees-steerable four-channel thruster in the stern. Conventionally, inland vessels tend to have one or more propellers, regularly ducted to protect them and to increase their performance, which are often placed in conjunction with multiple rudders to boost their manoeuvrability. Occasionally, the addition of an azimuth or bow thruster further improves their manoeuvrability [10,11]. In comparison, conventional marine vessels have propeller(s)-rudder(s) configurations positioned at the stern [12], and some marine vessels do carry more exotic propulsion systems such as transversal thrusters, azimuth thrusters, podded propellers, contra-rotating propellers, or even water jets for high speed vessels [13]. In addition, be aware that most small unmanned surface vehicles (USVs) also make use of propeller-rudder configurations [14]. However, some unmanned underwater vessels additionally utilise tunnel thrusters to increase their steering behaviour [15].
Within this aforementioned accumulation of more conventional actuation systems, the tunnel and azimuth thrusters seem to show the largest similarities with the thrusters of this study. The general design and performance of the tunnel thrusters are discussed in [16,17]. Whereas the authors in [18] review the different azimuth thruster designs, architectures, and working principles, the authors in [19] additionally review complementary thruster phenomena such as thruster-thruster or thruster-hull interactions, and the author in [20] published experimental data for an azimuth thruster during both static and dynamic operations. Unfortunately, the modelling literature for both systems is quite scarce. The existing research mostly uses the open-water propeller characteristics, as derived by [21], for its starting point on which additional flow phenomena can be added. For example, the authors in [22] built a mathematical manoeuvring model containing azimuth thrusters including thruster-thruster interactions. Similarly, the authors in [15] suggested a tunnel thruster model for an underwater vehicle exhibiting small yaw or pitch angles, and noted the scarcity of available models, probably due to the complexity of the flow phenomena at hand. Albeit certain parallels can be drawn between the actuation system in this study and these more conventional tunnel and azimuth systems, the usability and adequacy of these conventional modelling approaches need to be investigated. Therefore, this study aims to: (i) Detail the mechanical design of both the 360-degrees-steerable steering-grid thruster and the 360-degrees-steerable four-channel thruster nested inside the new fleet of Watertruck + barges. (ii) Discuss, list, and show experimental towing tank data from both thrusters at different azimuth angles and propeller speeds, at zero advance speed. (iii) Investigate the adequacy of the open-water propeller characteristics as the kernel of the theoretical thrust modelling approach. (iv) Incorporate azimuth-angle-dependent internal and external hull thrust deduction losses by means of an extension to the theoretical model from (iii) (v) Offer an additional artificial neural network modelling approach that might be useful to grip the current, and perhaps future, inherent complex flow phenomena occurring within both thrusters.
The successful completion of these aims offers the future researcher three different thrust characteristic descriptions (i.e., a data table (ii), an extended physically-based model (iii)+(iv), and an artificial neural network (v)) for both actuation systems (i). In addition to the insights that these data descriptions provide, their information can also be leveraged by advanced motion controllers, thrust allocation algorithms, and plant models. This paper continues as follows: first, Section 2 shows the constructed vessel and details its non-conventional over-actuated propulsion system. Afterwards, Section 3 expands the theoretical thrust model with angle-dependent thrust deductions and explains the multilayer neural network model structures. Subsequently, Section 4 lists the experimental towing-tank data and feeds these data to both modelling methods. Finally, Section 5 provides a comparison and discussion of the modelling results, and Section 6 concludes this study.

Material
This section starts with a brief description of the overall design philosophy of the new fleet of self-propelled inland barges in Section 2.1. Afterwards, it details the mechanical design of both the steering-grid and the four-channel thruster in Sections 2.2 and 2.3, respectively. Finally, it concludes with a detailed summary of the main differences between these thrusters and the more conventional actuation systems in Section 2.4.

Self-Propelled Inland Barges
During its first phase, the Watertruck + project will construct 12 unpropelled barges, 16 self-propelled barges, and 3 push vessels of European Class type I-II. Their final deliverable aims towards building 500 vessels in the European Class range I-IV [6]. In order to gain insight in the automation potential of these vessels, and in order to study the technological feasibility, and current or future limitations, of unmanned inland cargo vessels, the authors constructed a scale model of a self-propelled barge of European Class type I. The real-size counterpart has a length of 38.5 m, a beam of 5.05 m, and a maximal draught of 2.8 m, whereas the KU Leuven scale model has a similar geometry, scaled down with a factor of 8. Figure 1a shows four real-size Watertruck + barges of European Class type I, and one push vessel (slightly visible on the right), and (b) the KU Leuven scale model deployed as an unmanned, autonomous, inland cargo vessel. Figure 2a depicts a three-dimensional design drawing of this scale model where the embedded four-channel thruster (b), and steering-grid thruster (c) can be seen on the left-hand side and right-hand side of the vessel, respectively. This figure shows the vessel with open hatches in (a), which are the grey rectangular extensions on the top of the image, and sails with closed hatches in Figure 1b. As shown on Figure 2, the axes of the body-fixed reference frames which are used for all components throughout this paper point as follows: the x-axis towards the bow, y-axis towards starboard, and the z-axis towards the bottom. Both thrusters have a Kaplan Ka-series propeller without nozzle with a blade area ratio 0.65, and a pitch-diameter ratio of 0.95. The diameter of the steering-grid propeller measures 100 mm, whereas the diameter of the four-channel propeller measures 150 mm.

Mechanical Design Steering-Grid Thruster
The 360-degrees-steerable steering-grid bow thruster consists of two main parts: the propeller providing the thrust force, and the steering grid which orients the outflow of the water stream and thus provides the steering capabilities. Figure 3a shows a longitudinal cut of this thruster, i.e., along the xz-plane through the origin of Figure 2c. The propeller and its tilted shaft can be seen on the right side of the cut, where an anti-debris grid covers the inlet hole. The left-hand side shows the axis of the steering grid, but the grid itself is not shown to increase the image readability.  [7,8], and (c) tilted bottom view of the steering-grid thruster, adapted from [7,8]. Figure 3b shows an abstract top view of the bottom section of this thruster which illustrates the chosen angle convention of the internal control angle of the steering grid, θ g i . The drawn blue arrows show the orientation of the exiting water flow, which is opposite in direction compared to the thrust force, T g . It is assumed that internal control angle, θ g i , is equal to the output angle, θ g o , of T g , i.e., the flow exits the grid in alignment with the grid position in the xy-plane. Therefore, θ g i immediately gives the orientation of T g . On top of that, the steering grid has a static angle, α, of approximately 28°(curved surface), relative to the x-axis in the xz-plane, visible in Figure 2c. To illustrate this angle, the drawn blue arrows in Figure 3a show the exiting water flow in the xz-plane for θ can also be seen in in the xy-plane in Figure 3b on the left.

Mechanical Design Four-Channel Thruster
The 360-degrees-steerable four-channel stern thruster consists of two main parts: its propeller to induce thrust forces, and a 360-degrees-rotatable steering mechanism consisting of half a sphere with an opening of approximately 85°to orient its exiting water flow. Figure 4a shows a longitudinal section of this four-channel thruster, in the xz-plane in the origin of its body fixed reference frame. This section shows the propeller and its shaft together with the steering mechanism (of which only the bottom section is shown) which points towards the stern. This orientation would guide the exiting water flow towards the outlet channel on the left which points to the stern. The other shown outlet channel points to the bow, and the other two transversal outlet channels cannot be seen in this cut. Notice that the two shown outlets have a downwards angle of approximately 15°relative to the x-axis in the xz-plane, measured from the outlet of the internal steering mechanism-whereas the two transversal outlets have no downwards angle but are positioned a bit higher, directly in line of the internal outlet of the steering mechanism. This height-position difference of the outlets can be seen in Figure 2b.  i which orients the theoretical outflow of the water stream, denoted by the blue arrows. Due to the geometric differences of the thruster channels, one cannot assume that θ c i will equal the output angle, θ c o , of the resulting thrust force, T c . The main geometric design effects that cause this discrepancy between θ c i and θ c o seem to be: (i) the fact that the water flow can only exit through one (e.g., Figure 4d at θ c i = 180°), or a superposition of two (e.g., Figure 4d at θ c i = −45°) channels simultaneously, (ii) the different channel lengths and shapes (e.g., downwards bends), (iii) internal deflections of the water stream in the xy-plane when θ c i = 0°, 90°, 180°, or −90°, and (iv) a potentially remaining angular velocity of the water flow as there is no grid to align the outflow of the water. Therefore, in this study, the resulting four-channel thrust, T c , will be orthogonally decomposed in its longitudinal, T c x , and transversal, T c y , components for different propeller speeds, n c , according to:

Differences with Conventional Actuation Systems
In summary, the main design differences of both thrusters in this study compared to the more conventional thrusters appear to be the simultaneous occurrence of: (i) having the rotational axes of the propellers perpendicular (or almost perpendicular in case of the steering-grid) to the calm water plane, hence the propellers rotate in a plane parallel to the xy-plane, (ii) being completely nested inside the hull and thus not positioned directly in the ambient flow field, (iii) drawing in water from underneath the vessel hull, (iv) having both an inlet and an outlet positioned in the hull, not unlike jet systems, (v) having an internal deflection or rotation of the water stream from inlet to outlet, and (vi) having a 360-degrees-steerable outflow of the accelerated water provided by a steering mechanism, where the propeller itself does not change its position or orientation.

Modelling Structures and Identification Procedures
This study offers two different modelling approaches: one based on the open-water propeller characteristics, handled by Section 3.1, and one using artificial neural networks, explained by Section 3.2. Afterwards, Section 3.3 details the identification procedures for both modelling approaches.

Propeller Characteristics Thrust Model
This section uses the open-water propeller characteristics model of [21] as a starting point to construct new model structures. This open-water model takes a first order linearised lift force approximation of the propeller blades to calculate their theoretical thrust force. Given the fact that this lift force is a physical characteristic of the blade geometry, it is deemed suitable as a starting point to model both conventional and non-conventional actuation systems that use a propeller (for example, see its use for the tunnel and azimuth thrusters discussed in Section 1). Afterwards, this model will be expanded by adding the more complex flow phenomena that occur in the propulsion systems of this study. The model of [21] found the theoretical thrust force, T th , to be quadratically dependent on the revolution speed of the propeller shaft, n, and bi-linearly dependent on the shaft speed and axial inflow velocity, V A : It is common to introduce the non-dimensional thrust coefficient, K T , and non-dimensional advance ratio, J, to express the ship propeller performance, by introducing: where ρ is the water density and D p the propeller diameter. Moreover, for the conventional open-water propeller configurations, K T is often adequately approximated by a linear expression in J: Combining Labels (2), (3), (4), and (5) gives: which expresses the same information as Label (2) with: Furthermore, V A can be written as a fraction of the total speed of the vessel, U, by introducing the wake factor, w, which accounts for the speed reduction of the flow field at the inlet of the propeller: In addition, the available thrust force for the vessel, T, is a fraction of the generated T th due to propeller-hull-interaction losses, accounted for by the thrust deduction number, t, which is assumed to be speed-independent: In conclusion, by accumulating Labels (6), (8), and (9), T becomes: Note that (10) does not model the propeller dynamics, which can have a noticeable effect when manoeuvring at lower speeds [23]. Consequently, for some applications, it might be useful to extend T with these dynamics, as demonstrated in [24]. However, the experimental data of this study do not include transient dynamic measurements, and therefore their effects are omitted henceforth. Furthermore, as described in Section 4.1, the experiments were performed with zero advance speed, hence, K T = K T 0 . Combining this information, (10) tailors further down to: for this study. The model complexity of the thrusters in this study, which is summarised by Section 2.4, lies mainly in the angle-dependent thrust deduction, i.e., making t = f (θ i ). Moreover, both internal and external thrust deductions can occur as the water flows over the external hull, into the thruster inlet, and then through the nested thruster, and thus internal hull, itself, resulting in: where t i and t e respectively indicate the internal and external hull-interaction thrust deductions. Hence, (11) refines to: Finally, in order to avoid a model structure bias in the identification results, the theoretical thrust force will be investigated and modelled as an unknown polynomial function in n, T th = f (n) = T m (n). This approach does not force T th to equal T nn n 2 when U = 0. Accordingly, this method can thus validate whether the originally derived quadratic propeller-speed dependent thrust model indeed provides the best data fit. Ultimately, the thrust forces T to be identified according to the theoretic model in this study turn into:

Feedforward Multilayer Network Thrust Model
This section describes the architecture of feedforward multilayer neural networks that offer an alternative modelling approach to describe the complex flow phenomena discussed in this study. In general, artificial neural networks have one input layer, a number of hidden layers, and one output layer. However, a multilayer feedforward network with a non-polynomial activation function and only one hidden layer can already serve as a universal approximation for any nonlinear function [25]. This statement means that, regardless of the complexity of the flow phenomena that occur in our system, there will always exist network structures that can approximate the data sets. The capacity of these neural networks to model nonlinear characteristics drove the decision to add this alternative modelling approach.
In order to describe the general layout of a multilayer feedforward network, imagine the following example. If one takes an input vector, x ∈ R m , and an output vector, y ∈ R l , separated by one hidden layer with n h hidden neurons, and connects the hidden layer to y via W ∈ R l×n h and to x via V ∈ R n h ×m with the addition of a bias vector, β ∈ R n h , to the latter and an enveloping activation function, σ(·), the following matrix-vector expression will describe the model structure: which can also be expressed in an element-wise form, for i = 1, ..., l, according to: Figure 5 illustrates such a network hierarchy with the following model structure: two inputs (m = 2), one hidden layer of three neurons (n h = 3) of which the bias terms are not explicitly shown, two outputs (l = 2), and interconnection weights w ij , v ij . Be aware that these descriptions, (15) and (16), implicitly use a linear activation function for the output layer, although this is not necessary. Furthermore, the σ(·) can have any shape, but typically a nonlinear function, such as sigmoid, tangent, hyperbolic tangent, etc., is used in order to exploit their characteristics when one trains the network with experimental data. Note that, similarly to (15), a second hidden layer could be inserted, resulting in: with n h 1 and n h 2 neurons in the hidden layers, connection-weight matrices, W ∈ R l×n h 2 , V 2 ∈ R n h 2 ×n h 1 , V 1 ∈ R n h 1 ×m , bias vectors, β 2 ∈ R n h 2 , β 1 ∈ R n h 1 , and activation functions, σ 1 (·) and σ 2 (·). This hidden-layer insertion can be repeated to add multiple layers in an equivalent manner. The amount of hidden layers, their number of neurons, and the activation function types can all be chosen by the user. These design aspects offer great flexibility on the one hand, but an infinite amount of possible network architectures on the other hand. Accordingly, a few network designs will be shown in the results section which offers the possibility to compare their performance. These developed multilayer feedforward networks were trained to generate an output vector y (based on their input vector x), which aims to represent the desired output vector, i.e., the measured towing-tank data sets. More details of this network training procedure can be found in Section 3.3.2.

Model Identification Procedures
The theoretically derived thrust model of Section 3.1 and the neural networks of Section 3.2 provide a θ i -dependent and n-dependent model structure for the modelled thrust forces, T(n, θ i ), to which their respective experimental data sets, D(n, θ i ), can be fitted. Equation (18) shows the chosen cost function which calculates the least squares error, E, between the models and their data, i.e., between T(n, θ i ) and D(n, θ i ), for P data points. Sections 3.3.1 and 3.3.2 explain the procedure to minimize (18), and thus maximize the model fit, for each afore-mentioned model structure of which the maximal model complexity was iteratively determined based on its impact on the final residuals:

Identification Propeller Characteristics Model
The introduction of (14) into (18) produces the generic cost function for the model structures of Section 3.1: Within this study, the maximum polynomial orders for t(θ i ) and T m (n) were respectively chosen to be quintic and cubic; hence, in this situation, E would equal: (20) Furthermore, Section 4 also identifies the model structures of T m (n) = T nn n 2 or T nnn n 3 , given that the former aligns with the theoretical model of [21] (see Section 3.1), whereas for t(θ i ) all coefficients of its polynomial model will be identified. A Python (version 3.6) script performed the minimization of (20) by using the least_squares() function from the open source SciPy packages [26]. This function determined its step size of the optimization problem, i.e., changing the values of the coefficients of t(θ i ) and T m (n) in order to minimize (20), based on a trust region reflective method, based on the algorithm from [27].

Identification Feedforward Multilayer Network Model
The term feedforward refers to the calculation of y based on x, as these data run through the network of Figure 5 from left to right. Afterwards, the cost function (21) can be determined by the insertion of y from, for example (15) or (17), into (18): where y p represents the network output of the p-th data point, which would be generated by the p-th input point, x p . Within this study, the deepest developed network structure embodied two hidden layers; subsequently, its cost function can be defined by: Note that the p inputs for this network, x p , are the same inputs as for the data points D(n p , θ i,p ); hence, (22) refines to: Analogously to Section 3.3.1, this cost function (23) needs to decrease in order for the model fitness to increase. A common way to alter the distance between the data and the network is to use the back-propagation method as first suggested by [28]. This method derives an analytical expression of the gradient of the cost function for each neuron in the network by means of an iterative calculation from the right to the left in Figure 5, hence the name backpropagation, and uses this gradient to update the weights. Furthermore, one can augment this method by using the Levenberg-Marquardt algorithm to determine the updates of the connection weights [29]. This augmented back-propagation method was used in this study. Finally, in order to construct a network that generalises well to new inputs, one should avoid the over-fitting of data. Several regularisation methods exist to achieve this generalisation of the network. This study shows two regularisation methods: the early-stopping approach [30], and the Bayesian regularisation procedure of [31,32]. The early-stopping regulation divides x in a training, x t , and a validation, x v , set. Afterwards, it trains the network with the training data and tests the updated network with the validation data. When the error on the validation set, E(x v ), starts to rise, although the error of the training data, E(x t ), is still decreasing, the iterative algorithm stops, whereas, like the name suggests, Bayesian regulation uses the principles of the Bayesian probability theory as its regularisation principle. This regularisation is achieved by adding a term that penalises the size of weights, e.g., for a total of j weights w: to the cost function (23). In addition, within this approach, the Bayesian probability theory automatically embeds Ocam's razor principle by punishing too flexible and too complex module structures based on the Bayesian evidence principle [31]. Hence, this approach offers an automatic regularisation method, which boosts the generalisation of the trained networks. A collection of Matlab scripts constructed all the networks within this study. Consequently, the two afore-mentioned regularisation methods were applied by respectively calling the 'trainlm' [29], and 'trainbr' [33] function calls within the Matlab neural network environment. For the former, a data division of 70% training and 30% validation data was used. Note that this data separation is not necessary for the latter as it has an automatic regularisation approach as mentioned above.

Results
This section starts with listing and illustrating the raw experimental towing tank data from both thrusters in Section 4.1. Afterwards, Section 4.2 takes these data sets to identify the models described by (14) from Section 3.3.1, and similarly Section 4.3 trains the networks given by (15) and (17) from Section 3.3.2. All these results will be briefly discussed in each section, whereas an overall discussion and comparison between the results can be found in Section 5.

Experimental Towing-Tank Data
The experimental data for both thrusters were measured in a towing tank facility in deep water at zero advance speed, i.e., as a bollard pull test, and can be found in Section 4.1.1 for the steering-grid, and in Section 4.1.2 for the four-channel thruster.

Steering-Grid Thruster Data
The steering-grid thruster performed its bollard pull tests as a stand-alone device without enveloping ship hull hence only the internal angle-dependent hull losses occurred, i.e., t = t i = t i (θ g i ).
In addition, given the assumption that θ g i = θ g o , as mentioned in Section 2.2, the resulting thrust force data from the steering-grid thruster, D g , were measured in alignment with the angle of the steering grid, using a planar beam load cell of the type LCPB series of OMEGA with a maximal error of 1%. Table A1 summarises these different D g measurements at different θ g i and different rotational speeds of the propeller, n g , i.e., D g = D g (n g , θ g i ). Subsequently, Figure 6 outlines these data points on a polar plot, assuming symmetry over the longitudinal x-axes. The decrease in thrust force for θ g i ∈ [150, 180]°might be explained by the potential occurrence of a recirculation zone of water flow, as in this orientation the exiting water stream is pointed towards the inlet of the steering-grid thruster.

Four-Channel Thruster Data
The bollard pull tests of the four-channel thruster were performed inside half a ship hull, which was split in the transversal direction, i.e., cut in the yz-plane at midship. Consequently, under the assumption that the most crucial external thrust deduction effects, t e , will be captured by using half a ship hull, t(θ c i ) becomes the sum of t i (θ c i ) and t e (θ c i ). Given the information that for the four-channel thruster, θ c o = θ c i , as indicated in Section 2.3, the thrust forces were measured on both the longitudinal x-axis, D c x (n c , θ c i ), and transversal y-axis, D c y (n c , θ c i ), of the four-channel thruster for different θ c i and propeller speeds, n c . Therefore, similar to T c , the total resulting measured thrust force of the four-channel thruster, D c , can be derived from: Table A2 lists all the four-channel thruster data points, measured by a UDW3 force/torque sensor from AMTI, with a maximal error of 3% per component (mainly due to the possibility of cross-talk). Thereupon, Figure 7 illustrates these measurements in two different perspectives. Firstly, plot Figure 7a,b show the decomposed measured thrust forces for θ c i ∈ [0, −180]°at different n c . Note that the small fluctuations in n c , reported in Table A2, are accumulated around their mean value in order to plot these data sets with one line, resulting in for example: n c = 340 ± 25 rpm. Secondly, plot Figure 7c,d shows exactly the same data as plot Figure 7a,b, but transformed into θ c o and D c , calculated according to (25).  Figure 7c confirms the nonlinear mapping of θ c i to θ c o hypothesis, which appeared to stay consistent over different n c -sets. This n c -independent consistency in the angle mapping seems to hint at a structural, geometrical origin of this nonlinearity. The relatively constant θ c o for θ c i ∈ [−60°, −120°] seems to originate from the incapability of the thruster to produce a significant longitudinal force, i.e., D c x , within this angular domain allowing the transversal force D c y to dominate the thrust generation, which results in an almost transversal θ c o . The inverse seems to be true for θ c o when θ c i ∈ [−150°, −180°]. Given the currently available data sets, it remains hard to judge the explicit cause of these nonlinearities. However, Section 2.3 ended with a list of possible contributors to this nonlinear mapping.

Identification of the Propeller-Characteristics-Based Thrust Models
Plugging the data sets from Section 4.1 into the cost function of (20) allows for the identification of the modelled thrust forces according to (14). By varying the order of the polynomials of t(θ) and T m (n) up to their respective maximums of quintic and cubic, one can compare the resulting model complexities based on the residuals of their cost functions. This comparison offers a bias-variance trade-off during the eventual model selection [34]. The final model selection depends on the requirements of its end user. Section 4.2.1 calculates T g (n g , θ g i ) by using D g (n g , θ g i ), whereas Section 4.2.2 calculates the orthogonal decomposition of T c (n c , θ c i ) into T c x (n c , θ c i ) and T c y (n c , θ c i ) by respectively using D c x (n c , θ c i ) and D c y (n c , θ c i ).

Steering-Grid Thruster Identification Propeller-Characteristics-Based Thrust Models
Due to the absence of the hull in the towing-tank experiments of Section 4.1.1, no external hull thrust deductions occurred (except for the small external surface of the thruster itself which will be neglected here), i.e., t(θ g i ) = t i (θ g i ). This paves the way to identify t e (θ i ) by the future increase in thrust deduction when the actuation sits inside its hull by using : t(θ i ) = t i (θ i ) + t e (θ i ). Table 1 summarises the residuals of the cost functions after performing the afore-mentioned model fits. The absolute value of these residuals has no meaning; however, their relative size does as it indicates the relative model fitness of the underlying model structure. The columns denote the different model structures of T m (n g ), whereas the rows show the different polynomial orders of t(θ g i ). For example, a cubic-order model for t(θ g i ) in combination with T m (n g ) = T nn n g 2 can be identified by: (26) Here, as visible in Table 1, the final cost of E would equal 5.38. To further clarify this representation of the results, Table A3 lists all the identified model coefficients for all the model structures of Table 1. This information sheds more light on the cost residuals; for example, Table A3e shows negative values for T nnn which would mean a thrust deduction with increasing n g , hence here the data seem to be over-fitted and thus losing physical meaning. Table 1 also seems to validate that taking T m (n) = T nn n 2 generates a simple but robust model which also aligns with the theoretically derived propeller-characteristics model of Section 3.1, as adding a linear propeller-speed dependency, T n , to T m (n) results in marginally lower residuals. Nevertheless, this linear term could be added if desired. Given the adequacy of T m (n) = T nn n 2 , Figure 8a plots the identified quadratic, linear, and constant models for t(θ g i ) with T m (n) = T nn n 2 , together with the data points D g (n g , θ g i ) to illustrate the model fits. In addition, given the rather constant thrust decrease for θ g i ∈ [150, 180]°, one could also perform the above-mentioned modelling procedure of T g only for the θ g i ∈ [0, 150]°domain, and assume T g (n g , θ g i ) = T g (n g , 150°) to be a constant value for θ g i ∈ [150, 180]°. Figure 8b illustrates this second methodology for T g , using a quadratic fit for t(θ g i ). Figure 9 provides another way to visualise these models fits by plotting T g over its n g domain for θ g i = 0°, 90°, and 180°by using a quadratic model fit for t(θ g i ) and T m (n) = T nn n 2 for the aforementioned full θ

Four-Channel Thruster Identification Propeller-Characteristics-Based Thrust Models
The four-channel thruster generated its data sets embedded inside half a hull, making it impossible to diversify between t i (θ c i ) and t e (θ c i ) as they both depend on the same angle. Thus, this section identifies the thrust deduction coefficients as a superposition of both external and internal hull losses. Note that future towing-tank experiments without a hull, as have been done for the steering-grid thruster, could identify t i separately and thus also give t e . Table 2 lists the different residuals of the cost functions of the model fitting algorithm for both T c x in (a), and T c y in (b). This table also gives the impression that T m (n) = T nn n 2 forms a simple and robust model structure to represent the theoretical thrust force, which again aligns with the theoretical derivation of Section 3.1. Here, too, adding T n (or T n and T nnn ) only generates marginal improvements compared to the quadratic n c -dependency. Furthermore, this table indicates that the polynomial order of t(θ c i ) is better modelled by an odd function for T c x and by an even function for T c y .

Identification of the Feedforward Multilayer Network Thrust Models
Similarly to Section 4.2, this section plugged the data from Section 4.1 into different feedforward multilayer network designs. Within the subsequent figures and tables, the abbreviations 'br', and 'es' respectively denote a Bayesian regulation and early stopping regulation scheme. These abbreviations will be followed by two numbers which indicate the amount of neurons in each hidden layer (e.g., br5-0, will only have one hidden layer with five neurons, whereas br3-3 would have two hidden layers with three neurons each). Furthermore, all networks used a tan-sigmoid transfer function, tansig(·), for the activation functions of the hidden layers, i.e., σ(·) = tansig(·). Finally, note that all the training procedures of the artificial neural networks within this study were initialised with random weights for the network coefficients. This random initialisation can have an impact on the final results. Therefore, all networks were trained five times and the subsequent tables will show the average, minimum, and maximum final residual cost for each network structure over all five of these training cycles. Analogously to Section 4.2, the model structures can be selected based on a bias-variance trade-off which depends on the requirements of the end user.

Steering-Grid Thruster Identification Multilayer Networks
The networks for the steering-grid thruster used n g and θ g i for the input vector, x, as these inputs also generated the experimental data, i.e., D g (n g , θ g i ). The output vector of these networks, y, had only one element which constitutes the modelled thrust force, i.e., y = T g . To exemplify the cost function calculation and consequently the identification of its neural network, (27) calculates the residual for a network which models the steering grid with one hidden layer of three neurons by plugging (16) into (21), and thus by minimising: with n h = 3, m = 2, and l = 1 and thus its index is omitted. Accordingly, in matrix notation W ∈ R 1×3 , V ∈ R 3×2 , and β ∈ R 3×1 . The two input vector elements, i.e., x 1,p and x 2,p are n g p and θ g i,p , note that the subscript i from the latter refers to the internal actuation angle and not to an element notation. Depending on the selected training method, (27) can calculate the residual for es-3-0 or br-3-0, shown in Table 3 which lists the residuals of the final cost functions for all the identified network structures. Figure 12a plots the identification results of three different network topologies which only used one hidden layer which had two, three, or five neurons. Moreover, this figure compares the early stopping and Bayesian regulation approaches. Based on this plot, the Bayesian regulation outperformed the early stopping regulation. This performance difference can also be seen in Table 3a. In addition, Figure 12b shows similar identification results for two different network architectures which nested two hidden layers. Here too, one can observe the Bayes regulation method outperforming the early stopping approach, on which Table 3b

Four-Channel Thruster Identification Multilayer Networks
This section constructs different networks with two hidden layers and a Bayesian regulation approach for the four-channel thruster data. These specific network architectures were chosen based on their performance in the results of the steering-grid thruster, handled in Section 4.3.1. Although all the developed networks used n c and θ c i as their input values, three different approaches for the output, y, were selected to illustrate the network-architecture flexibility: (i) Fitting the orthogonally decomposed thrust forces separately, i.e., constructing two networks with one output each: one for deriving T c x based on D c x , and analogously one for finding T c y based on D c y .
(ii) Fitting the orthogonally decomposed thrust forces simultaneously, i.e., constructing one network with two outputs, T c x and T c y . (iii) Fitting the resulting thrust force and its orientation simultaneously, i.e., construction one network with two outputs, T c and θ c o . Table 4 compares the final cost function residuals of different network topologies for these three different cases. Firstly, regarding case (i), (a) summarises the residuals for T c x and (b) for T c y . Secondly, concerning case (ii), (c) shows the residuals of the chosen different network architectures to model this network which has two outputs. Thirdly, for case (iii), (d) shows a variety of model structures that represented this dual-output network. In addition, Figure 13 illustrates the fitness of several identified models for cases (i) to (iii). To plot these identified network structures, one of the five runs of the respective network designs was chosen at random. When one compares case (i) to case (ii), hence comparing the plots (a), (b) with (c), (d), the former seem to have less variance compared to the latter. However, caution should be taken when one compares these plots as they represent different network topologies, have different amounts of neurons and they are a random selection of one of the five network training cycles performed on each design.

Discussion
The pure experimental data of Section 4.1 already expose interesting results for both thruster designs. On the one hand, they show the potential occurrence of a recirculation zone of flowing water for certain positions of the steering grid. On the other hand, they highlight the rather biased output angles of the resulting thrust forces of the four-channel thruster. Nevertheless, additional experiments (e.g., setups with different channel geometries and lengths, propeller diameters, steering-mechanism designs, etc.) could provide complementary insights in the complex flow phenomena arising in these thrusters. In addition, experiments with U = 0 could identify the wake factor which might depend on both the vessel velocity and the thruster location, i.e., w = f (U, X t , Y t , Z t ). Similarly, these experiments could further investigate the adequacy of (5) for these non-conventional thrusters. Note that the ambient flow, when U = 0, does not decrease the thruster performance, but the complex interactions between the water outlet and the surrounding flow can decrease the effective thrust [35]. Furthermore, experiments with dynamic propeller speeds could verify the adequacy of the existing dynamic thrust models, e.g., [24]. Here, one should note that the ambient flow field tends to not have a large impact on the propulsion system dynamics [36].
Given the complexity of the flow phenomena at hand, this study offered two model structure methodologies. Firstly, Section 4.2 validated the theoretical propeller-characteristics thrust model structures of Section 3.1 for both thruster designs by means of a bias-variance trade-off. This trade-off indicates that the quadratic propeller speed dependency indeed forms a simple but generic description of the resulting thrust forces for both embedded thruster designs at zero advance speed, although more complex model structures could be selected if desired by the end user. Moreover, Section 4.2 further illustrated the successful modelling approach to capture the internal control-angle dependent thrust deductions that exist in both propulsion systems. Evidently, the aforementioned additional experimental setups could also help to further refine the propeller-characteristics based thrust model with the above-discussed possible extensions to this model. Secondly, Section 4.3 demonstrated the possibility to use the multilayer feedforward networks from Section 3.2 to grasp the complex water flows. Bearing in mind that these networks can serve as universal approximations of nonlinear functions, they offer the capability to model more complex flow phenomena that might occur during future experiments with dynamic propeller speeds or when U = 0. Figure 14 compares both model structure methodologies for the thrust forces generated by the four-channel thruster. Figure 14 (a) depicts the longitudinal modelled forces, T s x , whereas (b) shows the transversal modelled forces, T s y . For T s x , a quintic model describes t(θ s i ) in conjunction with T m (n) = T nn n 2 for the propeller characteristics-based model, whereas a Bayesian regulated network of two hidden layers, with three neurons each, forms the neural network structure. A similar approach describes the structures of T s y , only here t(θ s i ) has a quartic order. In both cases, i.e., (a) and (b), the neural networks seem to outperform the theoretically based models. Particularly in (b), the theoretical model generates a non-physical fit for θ c i ∈ [-150°, -180°] whereas the neural network does not. In (a), for θ c i ∈ [−60°, −120°] at n c = 1530 rpm, the theoretical model also offers an incorrect higher thrust force compared to the neural network which is closer to the data points. These differences could origin from the tendency of the higher order polynomials of t(θ c i ) to over fit to the data set, whereas the Ocam's razor principle from the Bayesian regulation helps to avoid this situation.
Finally, the current and future identified model structures can provide inputs to several applications. For instance, they can illustrate the achieveable forces for a certain thruster over its whole control domain as exemplified by Figure 15. Furthermore, these models can be used to augment certain links in the automation chain such as thrust allocation methods and motion controllers. For example, an optimal thrust allocation algorithm calculates the optimal thrust magnitudes and orientations for all onboard propulsion systems in order to achieve the desired force and moment output which the motion control calculated [37]. Several implementations of these algorithms exist, such as a quadratic programming approach [38,39], or by using a linear pseudo-inverse model [40]. However, these approaches tend to model the thrust forces as angle-independent, whereas the models of this study could be used to provide more accurate, control-angle dependent, representations of these forces. In addition, motion controllers could take the aforementioned and modelled constant thrust zones into account, similarly to incorporating dead-zones [41,42], further improving their performance. Furthermore, model predictive control architectures could also benefit from the suggested physically-based models as they can intrinsically integrate propulsion system limitations and states into their motion control design [14,43]. Moreover, these predictive controllers could also embed the developed neural networks [44].  Figure 15. Visualisation of the possible magnitudes and orientations of T g modelled by the quadratic thrust deduction model, i.e., polynomial order of t(θ g i ) equals two, and a pure quadratic n g -depedency, i.e., T m (n) = T nn n 2 , for a continuous area of (θ g i , n g ) pairs. The outer circumference of this plot shows the same information as the quadratic polynomial at 1500 rpm of Figure 8a.

Conclusions
In conclusion, this study offers the mechanical design and experimental bollard-pull data of both non-conventional thrusters on the one hand, and it provides two different modelling methods to represent these data on the other hand. The design and data themselves deliver more insight into the operational capabilities of the thrusters, and help to better understand the occurring flow phenomena inside these thrusters. The suggested modelling methods can be used to augment the automation level of the novel fleet of self-propelled inland barges. The theoretical propeller-characteristics based thrust model provides an adequate data representation which got confirmed by means of a bias-variance trade-off and can be extended with an angle-dependent thrust deduction function for both propulsion systems. In addition, a multilayer feedforward neural network with Bayesian regulation exemplifies a second modelling approach which is inherently capable to model complex nonlinear functions and could offer more protection against the non-physical over-fitting of data.

Nomenclature
(·) c superscript for four-channel thruster variables (·) g superscript for steering-grid thruster variables (·) x subscript longitudinal decomposition body-fixed frame (·) y subscript transversal decomposition body-fixed frame α static downwards angle steering-grid β bias vector σ(·) activation function θ i internal angle steering mechanism θ o output angle resulting thrust force n propeller speed n h amount of neurons in a hidden layer t thrust deduction number t i , t e thrust deduction number due to the internal, external hull t(θ i ) angle-dependent thrust deduction modelled as polynomial in θ i t 0 , ..., t 5 coefficients quintic polynomial of t w wake factor x input vector artificial neural network y output vector artificial neural network D experimental data E result cost function J advance ratio K T dimensionless thrust coefficient K T 0 , K T J constant and linear term of linear approxmiation K T T modelled thrust force T m (n) theoretical thrust force modelled as a polynomial in n T n , T nn , T nnn coefficients cubic polynomial of T th T th theoretical thrust force U total velocity vessel W, V connection-weight matrices X t , Y t , Z t three-dimensional position of thruster Appendix A