Empirical Modeling of 2-Degree-of-Freedom Azimuth Underwater Thruster Using a Signal Compression Method

This paper presents an empirical modeling of a 2-Degree-of-Freedom (DoF) azimuth thruster using the signal compression method. The thruster has a gimbal mechanism with two servo motors and generates thrust in arbitrary directions. This mechanism can reduce the number of thrusters in an underwater robot and contribute to compact design. When an underwater robot is controlled with azimuth thrusters, the influence from the rotational motion of the thruster has to be considered, and a dynamic model of the azimuth thruster is needed. It is difficult to derive an analytical model because the system model depends on complicated fluid dynamics. In this study, empirical models of force and moment for rotational motion were derived for practical use through frequency analysis. A signal compression method can effectively extract the system model in the frequency domain from just the mechanically constrained frequency response. Experiments were carried out using a force/torque sensor that was connected to a cantilever in a water tank. The system model was analyzed with Bode plots, and the model coefficients were derived through curve fitting. The derived model was verified by a validation experiment.


Introduction
A typical mechanism for the propulsion of underwater robots is a thruster with a propeller. The dynamic model of a tilting thruster is essential for designing a controller for underwater robots. Due to the complex effects of nonlinear fluid dynamics, studies have derived the dynamic model using experiments. Yerger suggested that the thrust force is proportional to the signed square of the propeller's rotational speed [1]. Healy proposed a nonlinear model considering the fluid dynamics of the motor model and propeller blade [2]. Bachmayer presented a combination model consisting of an axial flow model and a rotational flow model. [3]. Blanke developed a three-state model based on dimensionless propeller parameters, thrust coefficients, and an advance ratio [4]. Kim proposed an advanced three-state model considering the surrounding flow rate and inflow angle [5][6][7]. Yu [8] and Spearenberg [9] numerically analyzed the effect of thrust ducts on different boundary conditions. Song conducted Bollard-pull experiments to find the relationship between axial thrust and input while considering the effect of thrust-thrust and thrust-hull interactions [10]. Boem used a system identification process to derive an accurate dynamic model for commercially available fixed-pitch variable speed thrusters [11]. Odetti devel-oped a novel thruster based on a pump-jet [12] and applied it to ASV (Autonomous Surface Vehicle) in shallow water [13].
Most underwater robots are driven by fixed thrusters, but tilting thrusters have recently been used in underwater robots [14]. A tilting thruster can reduce the number of thrusters in an underwater robot and contributes to compact design [15]. When an underwater robot is controlled with tilting thrusters, the influence of the rotational motion of the thruster has to be considered, unlike a fixed thruster. The rotating thrust force can be calculated by multiplying the steady-state thrust force by a trigonometric function of the rotation angle [16].
When a thruster is operated with rotation motion, additional forces and moments are generated, and they can have a complicated relationship. Studies have numerically analyzed dynamic models according to the rotation angle of the azimuth thrusters [17,18], but a model that takes into account the effects that occur while the thruster is rotating is needed. In our previous research, dynamic models for rotation motion were derived by oscillating a 1-Degree-of-Freedom (DoF) tilting thruster using sine waves [19]. Modeling errors occurred as a result of relying on discrete data in the frequency domain, and it was necessary to design a model-free robust controller to overcome them [20].
This paper presents a 2-DoF azimuth thruster, as shown in Figure 1. Like fish fins, the thruster has a 2-DoF direction of propulsion. There are two rotational joints with the gimbal mechanism. It can produce effects similar to those of bioinspired robots that produce twoway propulsion through tail motion [21]. In a steady state, the direction of the thrust is the same as the direction of the thruster. However, the thrust is bent and dispersed during the rotational motion of the azimuth thruster, and a reaction moment also occurs. It is difficult to analytically describe this phenomenon, because the system model highly depends on complicated fluid dynamics. In this study, a signal compression method was applied to derive an empirical model of the azimuth thruster. Signal compression methods are often used for the empirical modeling of mechanical systems and can extract the frequency response of a system in the whole range of mechanically constrained frequencies using a pseudo-impulse signal [22]. Signal compression methods have been used for the cylinder modeling of an automatic excavator [23] and friction modeling of a robot manipulator [24]. In repeated experiments, the velocity profile of the pseudo-impulse input was applied to the two rotational axes, and the frequency response of the longitudinal thrust force, lateral thrust force, and reaction moment were analyzed with a Bode plot. Model expressions were derived through a curve-fitting technique. The validity of the derived model was verified by comparing the experimental results with model-based simulation values. This work provides references on how to derive empirical models of vectoring thrusters that can be applied to a variety of underwater robots.
The remainder of the paper is organized as follows. The prototype of the 2-DoF azimuth thruster and the experimental apparatus are described in Section 2. The signal compression method for the empirical modeling is introduced in Section 3. In Section 4, the empirical models of thrust force and reaction moment are presented with an analysis of the Bode plot of the pseudo-impulse response. Section 5 shows a comparison between the model and experimental results. The validity of the empirical model is also discussed. Our conclusions are given in Section 6.

Description of 2-DoF Azimuth Thruster
A prototype of a 2-DoF azimuth thruster was fabricated, as shown in Figure 2. The system consists of a thruster and two servo motors and has a gimbal mechanism for 2-DoF rotational motion. A counter mass is used so that the center of mass and the rotational axis intersect. The longitudinal thrust that is generated by the thruster and the tilting angles of the two servo motors are system inputs, and the actual thrust force and reaction torque are outputs of the system model. The thruster (RCD-MI60, 24VDC, RHINCODON, China) has a maximum thrust force of 5.5 kgf. The duct diameter is 95 mm, and the weight in air is 850 g. The gimbal mechanism is actuated by two BLDC (Brushless DC) motors (EC-max 30, 24VDC, Maxon, Switzerland). The total weight of azimuth thruster system is 8.56 kg. The generated thrust force and reaction moment are shown in Figure 3. The X-Y-Z frame is a fixed frame, and the x-y-z frame is a moving frame that is attached to the tilting thruster. The thrust force is bent and dispersed due to rotational motion, and it can be decomposed into f x , f y , and f z in the moving frame. These forces can be converted to f X , f Y , and f Z in the fixed frame using a rotation matrix. In this study, the frequency response of forces was analyzed in the moving frame. The reaction moment was analyzed within the framework of the gimbal mechanism. The reaction moment M Z generated by the rotational motion of the 1st axis is defined in the fixed frame, and the reaction moment M y for the rotational motion of the 2nd axis is defined in the moving frame. A diagram of the experimental apparatus is shown in Figure 4. In order to measure forces and moments, the 2-DoF azimuth thruster unit was installed at the end of a cantilever with a 6-DoF force/torque sensor (RFT80-6A01, ROBOTUS, Gyeonggi, Korea) attached at the clamped end. The thruster is driven by a pulse-width modulation (PWM) signal with a range of 55-95%, and the BLDC motors are driven by position controllers (EPOS2 25/2, Maxon, Switzerland). A real-time controller (NI CompactRIO, National Instruments Corp., Austin, TX, USA) controls the actuation and receives the values from the force/torque sensor. The measured forces and torques from the sensor are converted to the values of the system origin through the adjoint matrix of the frame transformation.

Signal Compression Method
Signal compression methods are useful tools for the frequency analysis of mechanical systems under physical limitations. The impulse response can show the dynamic characteristics of a system in the whole range of frequency. However, many mechanical systems cannot make an ideal impulse input. In the signal compression method, a pseudoimpulse input that is expanded with phase delay is applied to the system. The output that is obtained from the pseudo-impulse input is mathematically equivalent to the impulse response. After the equivalent output is compressed, the frequency response of the dynamic system can be observed within the desired frequency range. A diagram of the signal compression method is shown in Figure 5. In this study, the frequency response was analyzed for system modeling with the Bode plot before transformation into the time domain. The dynamic characteristics due to rotational motion are affected by the angular velocity, acceleration, and higher-order terms. In this study, the angular velocity was input to observe the system characteristics. The power spectrum of ideal impulse input is constant in the whole range of frequency. The power spectrum within the desired frequency was designed as shown in Figure 6a, considering the capability of servo motors. It can be derived by a function of the power spectral density as follows: where Q and a are parameters, and n is the number of data values. The range of frequency of the power spectrum is manipulated using the parameter a. In this study, a was set as 450, and the range of frequency was designed as up to 10 Hz, as shown in Figure 6a, considering the capability of the servo motors. The pseudo-impulse input was derived by signal expansion as follows: H(jn) is a function of the filter for signal expansion with a time delay. The time delay depends on the parameter b, which is 2100 in this case. Through the inverse fast Fourier transform (IFFT), the expanded pseudo-impulse signal can be obtained in the time domain as shown in Figure 6b, and it is applied to the system. The equivalent output is generated as shown in Figure 6c, and then the impulse response can be obtained via signal compression with an inverse filter. Finally, the Bode plot shows the frequency response of the system, as shown in Figure 6d.

Empirical Modeling
Variations of the thrust force and moment due to the rotational motion of the azimuth thruster were modeled based on the steady-state thrust force. Longitudinal thrust force is generated by a thruster, as shown in Figure 7a. The steady-state thrust force is proportional to the duty ratio of the PWM from the motor drive, as shown in Figure 7b. The empirical model expression of the steady-state thrust force for PWM is as follows: The steps of empirical modeling are as follows: (1) Repeat the rotational motion experiment, increasing the PWM by 5% from 55% to 95%. (2) Measure the force and moment by entering the velocity profile as shown in Figure 6b on the 1st and 2nd axes, respectively. (3) Calculate the force of the moving frame reference according to the angle of each axis using the adjoint matrix. (4) Draw the Bode plot of the pseudo-impulse response through the signal compression. (5) Derive the model through the curve fitting with the numerical method to maximize the correlation number.
As a result of the repeated experiments, waveforms for x-direction attenuation and yand z-direction generation of propulsion were similar, and the reaction moments showed different forms of waveforms. The models can be derived by the analysis of the Bode plot with the corner frequency and the slope of asymptotes [25], as shown in Figure 8. As the phase plot is noisy, we focused on the magnitude plot. The model expressions of the force and moment are defined as follows: The force and moment models have the same system order but slightly different equations. The coefficients of the model expression were obtained through numerical optimization that maximizes the correlation coefficients. The results of experiments with significantly different tendencies through repeated experiments were eliminated. The objective function of modeling is defined as the maximum correlation coefficient with repeated experimental results, and uncertainty in the data is handled.

Model for 1st-Axis Rotation Motion
When we placed the velocity profile of the pseudo-impulse input on the 1st rotational axis, the compressed response showed variation of the longitudinal thrust force in the x-direction, the lateral thrust force in the y-direction, and the reaction moment in the Z-direction. For each PWM, curve fitting was performed for three or more experimental values to extract the results with the highest correlation number, as shown in Tables 1-3. As an example, Figure 9 shows the Bode plot of the experimental values and the proposed model of F 1x (s) and F 1y (s) at 75% PWM.   The variations in each parameter for PWM allowed us to analyze the characteristics of the model, as shown in Figure 10. The average correlation number of F 1x (s) is 0.60. It is not a high value but is within a reasonable range. The gain K increases with increasing propulsion. Differences in parameters d and g indicate significant differences in the natural frequency of the longitudinal thrust force and the lateral force.

Model for 2nd-Axis Rotation Motion
When we placed the velocity profile of the pseudo-impulse input on the 2nd rotational axis, the compressed response showed the variation of the longitudinal thrust force in the x-direction, the lateral thrust force in the z-direction, and the reaction moment in the y-direction. The 2nd-axis rotational motion did not rotate the gimbal structure, unlike the 1st-axis rotational motion. When the moment was measured at the end of the cantilever, the moment in the y-direction was mixed with the force in the x-direction. It was too small to separate a meaningful value from the measurement. In this study, the reaction moment in the y-direction was neglected.
The model parameters for each PWM are shown in Tables 4 and 5. As an example, Figure 13 shows the Bode plot of the experimental values and the proposed model of F 2x (s) and F 2z (s) at 75% PWM. The tendency of the model parameters according to PWM was analyzed, as shown in Figure 14. Like in the 1st-axis motion, the gain K increases with increasing thrust. The model for the change in thrust forces differs somewhat from the 1st-axis motion. It can be seen that the movement of the gimbal structure has a considerable influence on the change in propulsion.

Comparison between Model and Experimental Results
To validate the feasibility of the derived model, we compared the experimental result for an arbitrary rotation path with the model-based simulation value. It is recommended that the longitudinal and lateral thrust force be interpreted on a moving frame basis for dynamic interpretation, but as an application, the forces of the azimuth thruster are converted into forces in the fixed frame based on the location of the attachment. We derived the forces and moments in the moving frame by entering the inputs applied in experiments into a simulation model on two axes. The forces were converted into the fixed coordinates through a rotational matrix. The Z-axis moment was directly compared with the experimental results without coordinate transformation. The simulation was performed in MATLAB (Mathworks Inc., Natick, MA, USA).

Verification Test
A verification experiment was conducted by entering commands such as step inputs of different angles on the two rotational axes. The actual rotation angle for the step input is shown in Figure 15. The first axis rotated 45 • , and the second axis rotated 30 • . The thruster was operated in advance with 65% PWM to start the rotational motion in a steady state. First, the initial offset due to gravity and disturbance was calibrated, and then the forces and torque were measured.

Simulation Based on the Proposed Model
As in the experiment, we simulated the response to a step input. Since the derived models are for rotational speed, they were transformed into models for the angle of rotation as follows: M(s) θ(s) = Ks 3 s 4 + as 3 + bs 2 + cs + d (8) Assuming that the forces due to the rotation of each axis are independent of each other, the longitudinal thrust was calculated by superposing the change in axial force of each axis onto the steady-state force as follows: The sum of the models of the derived x-axis force variation allowed us to obtain the attenuated x-axis thrust force. The rotation of the coordinate system by two axes is expressed in a matrix as follows: The results of force simulation in moving coordinates can be transformed into forces in fixed coordinates as follows: The experimental results and simulated values are compared in Figure 16. The desired value is the force in the fixed coordinate system for the actual rotation angle of the two axes when there is only ideal longitudinal propulsion. In other words, the desired forces are the decomposition force to the rotation angle measured by the encoder values of the two servo motors, based on the propulsion model of the fixed thruster. These are ideal values that do not consider force attenuation or tangential force due to the rotation of the thruster. In the initial position, the propulsion occurs in only the X-axis direction, but as the two axes rotate simultaneously, the direction of the propulsion changes. For each force and moment, the steady-state error and root mean square error (RMSE) are defined as performance indices and are shown in Table 6.  The simulation results from the derived model followed the waveform changes in the experimental results well. The change in force on the rotation of the two axes was modeled, but the model did not include high-frequency oscillations from the rotation of the propeller. As the motor encoders measuring the rotation angle are very precise, the steady-state error observed in the Xand Z-axial forces can be considered as interference from forces from the water surface and water tanks. Overshoots in the experimental values are caused by the vibration of the gimbal structure. It was observed that more vibrations caused by step inputs occurred than in the modeling experiments.

Conclusions
This study showed an empirical model of the force and moment of a 2-DoF azimuth thruster. A signal compression method was used to extract frequency response characteristics within mechanical limits. The longitudinal thrust and lateral thrust and reaction moments generated by the rotation of each axis were measured through repeated experiments, and the model was derived by curve fitting to the Bode plot obtained from the experimental results. We validated the derived model by comparing the simulation values with the experimental results on the step inputs of the two axes. From the comparison results, the causes of modeling errors were analyzed. This study shows the trend of change in coefficients of the azimuth thrust model for thruster force, and normalization for the sensitivity analysis of model coefficients will be studied in the future.
The experimental results showed that the gimbal structure has a significant impact on the force and moment of the thruster. The difference between the first-axis rotation model and the second-axis rotation model was also due to the gimbal structure. A design that minimizes the gimbal is needed when fabricating an azimuth thruster. The modeling procedure presented in this study could be applied to develop vectoring thrusters that control the direction of the propellant. It could also help to design controllers for underwater robots with azimuth thrusters.