$$\boldsymbol{\varepsilon} = \begin{bmatrix} \mathbf{s} \\ \boldsymbol{\varepsilon} \end{bmatrix} = \mathbf{R}\_a^T (\boldsymbol{\alpha}\_k) (\boldsymbol{p} - \boldsymbol{p}\_k) \tag{8}$$

where *RT α* is the transformation matrix from the inertial to the path-fixed frame given by

$$\mathbf{R}\_a^T(\boldsymbol{\alpha}\_k) := \begin{bmatrix} \cos \alpha\_k & \sin \alpha\_k \\ -\sin \alpha\_k & \cos \alpha\_k \end{bmatrix} \tag{9}$$

*s* is the along-track distance, and *e* is the cross-track error (see Figure 10). The cross-track error is defined as the orthogonal distance to the path tangential reference frame. From (8) and (9), the cross track error can be computed by

$$\boldsymbol{\pi} = - (\mathbf{x}\_n - \mathbf{x}\_k) \sin \boldsymbol{\alpha}\_k + (\boldsymbol{y}\_n - \boldsymbol{y}\_k) \cos \boldsymbol{\alpha}\_k. \tag{10}$$

Its time derivative, which will be used later, is consequently given by

$$
\dot{\boldsymbol{x}} = -\dot{\boldsymbol{x}}\_n \sin \alpha\_k + \dot{\boldsymbol{y}}\_n \cos \alpha\_k. \tag{11}
$$

From (1) and (2),

$$\begin{aligned} \dot{x}\_n &= -\mu \cos \psi - \upsilon \sin \psi \\ \dot{y}\_n &= -\mu \sin \psi + \upsilon \cos \psi \end{aligned} \tag{12}$$

so that

$$\dot{e} = -\left(\mu\cos\psi - \nu\sin\psi\right)\sin a\_k + \left(\mu\sin\psi + \nu\cos\psi\right)\cos a\_k. \tag{13}$$

The control objectives are formulated to drive the cross-track error to zero by steering the robot and controlling its forward speed. With the lookahead-based steering method, a fixed parameter, known as the lookahead distance Δ, which corresponds to a distance along the path ahead of point *s*, is used to define the LOS vector. The robot is steered so that its velocity vector is parallel to the LOS vector (Figure 10). The resulting velocity vector will have a component perpendicular to the path, driving the robot towards the path until the LOS vector is parallel to the path so that *e* → 0. Then, the *Guidance module* calculates the linear and angular velocities, which are sent to the Husky controller and afterward to wheel motors to let the robot move automatically between goal points. The *Guidance module*

sends the linear and angular velocity inputs to move the Husky robot automatically inside the working environment.

#### **3. Prototyping and Preliminary Experimental Tests**

The conceived and designed layout has been implemented, and first navigation tests have been performed. The ground robot platform is mechanically designed to handle suitable sensors. Thus, the mechanical and electronic systems used for our experiments were prototyped and integrated on top of the main platform. The final configuration of the mobile robot for the first experimental tests is shown in Figure 11a,b. The second control box arrangement is reported in Figure 11c. The main ambient awareness sensors installed on the robot are a 3D LiDAR Velodyne VPL16 and an Intel Realsense D455 camera. The robot is equipped with a mini-ITX computer (see [36]), while the second control box embedded computer is a Jetson TX2 platform. The preliminary test targeted an emulated orchard and a safe navigation along a collision-free path given multiple goal points.

**Figure 11.** Mechanical and electronic configuration of the experimental platform (UR5 manipulator not shown): (**a**) front view, (**b**) back view, and (**c**) second control box.

#### *3.1. Emulated Scenario and Test Inputs*

For the initial evaluation of the implemented framework, a lab-emulated orchard was created in which parallel rows of tree were considered fixed obstacles on the map. The experimental emulated scenario is shown in Figure 12a, and the Gazebo digital twin is shown in Figure 12b.

The length and width of the environment are fixed to 13.2 m and 6.6 m, respectively. For simulating an orchard, the spacing between rows is set to 2.2 m (see [49]). The length of each row (static obstacles) is fixed at 2.2 m long. The robot's working space is decomposed into rectangular grid cells, and each grid cell represents 10 cm in the robot's working environment. The goal tolerance was set to 0.1 m, and the heading tolerance was set to 0.1 rad. The UGV chosen localization system was based on a filtered odometry approach; the Intel RealSense D455 *RGB*-*D* camera or the Velodyne VLP-16 3D LiDAR were used (in different tests) for obstacle detection, and the extended *BNM* was adopted for path planning. The vision system was set to operate in front of the robot, setting a view size of 5 × 6 × 1.5 m (xyz coordinates), and the obstacles were boxes of different dimension. The bounding boxes were generated with an offset of 0.2 m to ensure an additional safety. For path-tracking, the lookahead distance was set to 1 m for the LOS-guidance method, and the heading error was set to 0.05 rad.

**Figure 12.** Robot working environment: (**a**) the real environment; (**b**) the simulated environment using Gazebo.

Given the input data to define the goal points and the working environment, the map generator node creates a map. Afterward, the map data are processed with the goal points through the path planner to compute a collision-free path between each pair of goal points.

The planned path between each two goal points formed by a discrete sequence of poses, i.e., waypoints, is fed into the *Guidance module*. At the same time, the output of the *Guidance module* is mapped to the wheel motors by a simulated low-level controller running onboard.

#### *3.2. Experimental Tests and Discussion*

As per the envisioned tasks, the mobile robot was sent to the proper orchard locations to perform predefined tasks. Then, we experimentally tested the system to evaluate the automatic navigation performance in a generated given-map environment. In the first test configuration, we evaluated the automatic navigation performance of the robot under static environmental conditions. Then, in the second test configuration, we evaluated the robot's navigation performance in a dynamic environment, i.e., with an unexpected obstacle along the path. During the experimental tests, we verified that the proposed framework generates an appropriate path, and the robot avoids static and dynamic obstacles without collisions.

Figures 13 and 14 show the planned path for a mobile robot using *BNM* in two testing scenarios as well as the performance of the navigation strategy in both testing scenarios. In the first test configuration (see Figure 13, the robot starts moving from initial position *g*0 (0, 0) to visit three fixed goal points (*g*1, *g*2, and *g*3) added to the pre-built map. The goal points were located at (−2.2, −2.2), (8.8, 2.2), and (4.5, −2.2). The red cells express the fixed obstacles, blue dashed lines represent the planned path of the robot, yellow circles in the figure define the goal points, and green cells represent the area related to the safety margin introduced to allow the robot to safely navigate and steer. In this study, we adopted a safety margin around obstacles to avoid the possibility of overlapping the paths traced by the robot overlapping with obstacles. As it can be appreciated in Figure 13a, the robot might move very close to the obstacle. Therefore, a certain safety margin for anti-collision has to be assured. Considering the robot dimensions and footprint, it has been defined with a constant-size of six grid cells. All fixed obstacles are given in parallel lines, and the goal points are given between obstacles. The robot localization sensors are utilized to determine where a mobile robot is located within the environment. In this test, a filtered-odometry, i.e., encoders and IMU, approach has been adopted. The current coordinates of the robot are compared with the predefined ranges of the fixed goal point, and if the robot is within the range, the coordinates of the next destination point are considered.

**Figure 13.** Simulated environment with fixed obstacles: (**a**) planned path using *BNM*; (**b**) path executed by the robot.

**Figure 14.** Simulated environment with fixed and dynamic obstacles: (**a**) planned and replanned path using *BNM*; (**b**) path executed by the robot.

In the second scenario, the experimental test was carried out with the addition of a dynamic obstacle, as shown in Figure 14. The robot starts moving from the initial position *g*0 (0, 0) to visit a fixed goal point (*g*1) located at (8.8, 0.0), and then the robot returns to the initial point. Since the robot sensors detect a new obstacle that occludes the pre-planned path, the *Planning module* generates a new path from the current position to the destination point (*g*1).

The plots show that the robot's executed path is coherent with the planned one. The mobile robot shows appropriate navigation accuracy and performs smooth motion along the environment with an average error ≤ 10 cm. The traveled path mostly depends on the robot's kinematics, the accuracy of the sensors, and the environment. High navigation speeds lead to slightly higher errors. The proposed *ROS* package always directs the robot to follow the planned path, and the robot moves close to the calculated path. The path starts from the robot's initial position, passes through the intermediate goal points, and then returns to the initial point. A final path is assembled by connecting the paths between goal points in an iterative way until the path is completed, and the length of the path is the sum of the lengths of the goal-to-goal paths. From the first experimental scenario and reported experimental results, the generated path length is 24.3684 m, close to the 24 m of the ideal *BNM*-planned path, thus showing good performance. Screenshots of the experimental tests at different locations are presented in Figure 15.

**Figure 15.** From (**a**) start to (**j**) end: sequence of screenshots of the simulated orchard environment with obstacles (*crop rows*) showing the navigation test of the robot prototype in the physical environment, with the fixed obstacles represented by yellow lines, and the dynamic obstacle, not present in the predefined map, represented by the box

Despite the overall good experimental results, the system may fail in certain circumstances. The perception module, for example, may fail detecting objects that have a size much smaller or much bigger than that of obstacles assumed during parameter tuning. This is due to the fact that when defining a cluster, the algorithm uses a *min* and *max* threshold value to constraint the number of points a cluster can contain. Additionally, because the planning module works based on the grid cells, increasing the cell size will reduce the system's accuracy. In the same way, the guidance model may fail to control the robot if we set a very small value for the acceptance radius. In this scenario, due to the noisy data from the localization sensors, the robot may never reach the goal area, and thus it would not switch to the next goal point.

#### **4. Conclusions**

Here, the conceptual design and experimental implementation of a modular mobile robotic system for agricultural field tasks is presented. Based on the functional design concept *HW* and *SW* modules are conceived, developed and integrated on a Husky unmanned ground vehicle. In particular, the software modules are implemented using *ROS*, following a high and low-level multi-layered approach, and allow the use of different sensors for localization and obstacle detection. Localization is implemented using a filtered-odometry method that can be integrated with indoor or outdoor positioning systems. Both an RGB-D camera and a 3D LiDAR sensor are configured for use on the platform, and either can be employed for obstacle avoidance. An extended Boundary Node Method, which has been specifically adapted for this work, is used for path planning. A modified Lookahead-based Line of Sight guidance algorithm is utilized for path following. Preliminary field tests in an emulated orchard scenario demonstrate that the system can perform path-following with a suitable accuracy, and obstacle detection, with path re-planning when needed. Planned future work will include intensive experimental testing, detailed data analyses, as well as performance assessments of the system's capabilities of executing the three *PA* scenarios (monitoring, harvesting, and spraying) in both emulated and real environments.

**Author Contributions:** Conceptualization and methodology, R.A.S., G.T., K.D.v.E. and R.V.; software, validation, and formal analysis, R.A.S. and G.T.; resources, K.D.v.E. and R.V.; writing—original draft preparation, R.A.S., G.T. and R.V.; writing—review and editing, all; supervision, G.C., K.D.v.E. and R.V.; project administration and funding acquisition, K.D.v.E. and R.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported in part by the "Reconfigurable Collaborative Agri-Robots (RE-COARO)" Südtirol/Alto Adige 4th Call (project #4122) and by the European Regional Development Fund (ERDF), FiRST Lab Project #FESR1084.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors are grateful to Matteo Malavasi, who designed and fabricated some of the electronic systems of the robot, and Josef Zelger, who designed and fabricated the mechanical superstructure.

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

#### **References**

## *Article* **An Intelligent Predictive Algorithm for the Anti-Rollover Prevention of Heavy Vehicles for Off-Road Applications**

**Antonio Tota 1,\*, Luca Dimauro 1, Filippo Velardocchia 2, Genny Paciullo 3 and Mauro Velardocchia 1**

**Abstract:** Rollover detection and prevention are among the most critical aspects affecting the stability and safety assessment of heavy vehicles, especially for off-road driving applications. This topic has been studied in the past and analyzed in depth in terms of vehicle modelling and control algorithms design able to prevent the rollover risk. However, it still represents a serious problem for automotive carmakers due to the huge counts among the main causes for traffic accidents. The risk also becomes more challenging to predict for off-road heavy vehicles, for which the incipient rollover might be triggered by external factors, i.e., road irregularities, bank angles as well as by aggressive input from the driver. The recent advances in road profile measurement and estimation systems make road-preview-based algorithms a viable solution for the rollover detection. This paper describes a model-based formulation to analytically evaluate the load transfer dynamics and its variation due to the presence of road perturbations, i.e., road bank angle and irregularities. An algorithm to detect and predict the rollover risk for heavy vehicles is also presented, even in presence of irregular road profiles, with the calculation of the ISO-LTR Predictive Time through the Phase-Plane analysis. Furthermore, the artificial intelligence techniques, based on the recurrent neural network approach, is also presented as a preliminary solution for a realistic implementation of the methodology. The paper finally assess the efficacy of the proposed rollover predictive algorithm by providing numerical results from the simulation of the most severe maneuvers in realistic off-road driving scenarios, also demonstrating its promising predictive capabilities.

**Keywords:** rollover detection; heavy vehicles; off-road applications; predictive algorithms; artificial intelligence; load transfer ratio

#### **1. Introduction**

Nowadays, vehicle rollover detection and prevention are two critical aspects that must be taken into account for the safety of car passengers and pedestrians to avoid dramatic fatal crashes [1] and accidents in the urban scenario. An analysis of the sequence critical events leading to loss of control, and hence to rollover of vehicles equipped with Electronic Stability Control (ESC) Systems is provided in [2], while typical scenarios and characteristics of rollover accidents are investigated in [3], referring to the test criteria of National Highway Traffic Safety Administration (NHTSA) [1,4]. Vehicles with a high position of the Centre of Gravity (CoG) are more prone to rollover, as it happens for buses [5,6], heavy commercial vehicles [7,8], and articulated heavy vehicles [9–14]. In a similar way, two-wheel vehicles are affected by stability issues [15,16]; hence, their components must be dynamically optimised [17], through a component to assembly dynamic analysis [18,19]. The rollover prevention is a typical challenge for the automotive sector, especially for fuel cell and hydrogen trucks, where the safety requirements demanded for their storage systems are stricter than for light-duty vehicle segments [20].

**Citation:** Tota, A.; Dimauro, L.; Velardocchia, F.; Paciullo, G.; Velardocchia, M. An Intelligent Predictive Algorithm for the Anti-Rollover Prevention of Heavy Vehicles for Off-Road Applications. *Machines* **2022**, *10*, 835. https:// doi.org/10.3390/machines10100835

Academic Editors: Marco Ceccarelli, Giuseppe Carbone and Alessandro Gasparetto

Received: 31 July 2022 Accepted: 8 September 2022 Published: 21 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The root cause of rollover can be addressed to external factors, i.e., road irregularities and its bank angle, which also influence the powertrain management [21], and to aggressive input from the driver on the steering wheel. Considering these factors, rollovers are sometimes categorized, as conducted in [22–24], into two main types: un-tripped rollovers and tripped rollovers. The former [25] is related to driver fast maneuvering on smooth roads, while the latter [26] occurs due to sudden impacts that may apply lateral or vertical forces to the vehicle tyres [27,28], hence possible root causes are impacts with guardrails or hits with road objects as curbs and bumps [29]. In the last decades many researchers have studied the vehicle rollover risk, whose knowledge is mandatory for developing practical rollover prevention systems, which may be implemented using different control logic, through a two-stages process. Indeed, the first step is the detection of rollover risk, whereas the second is its mitigation and control, e.g., through an active handling suspension controller [30] or by means of a Model Predictive Control (MPC) strategy through active front steering [31], or using a braking system for vehicle dynamics control [32]. Several rollover indices (RIs) have been proposed in literature, starting from static rollover index, as the Static Stability Factor (SSF) [33–35], which only considers the geometrical parameters of vehicles. The vehicle dynamics [36,37] plays a key role in the rollover risk phenomena, and hence also for the indicators adopted to prevent it, which consider the most important vehicle states for rollover, such as roll angle and roll rate [38–40], lateral acceleration [41–44], side slip angle [45], yaw angle [46] and yaw rate [47,48]. The load transfer ratio (*LTR*) [49–51], which is based on the computation of vertical tire forces, is one of the widely used RIs for dynamic simulations, and shows a direct measure of how close the vehicle is to rollover [49]. A predictive *LTR* [52] has been proposed to provide a better prediction of rollover propensity w.r.t traditional *LTR*, while in [53] the contour line of load transfer ratio (*CL*-*LTR*) is proposed for an accurate prediction of vehicle rollover threat. Moreover, time to rollover (*TTR*) [39,43] and the rollover index (RI) [38,54] are used widely for the rollover risk detection. At a glance, many RIs have been proposed in literature, with slight differences on the factors affecting the rollover, thus some RIs are more adapt for specific rollover conditions and not work well in other situations. For example, an extensive research investigation is focused on the untripped rollover risks and only a minimal attention is reserved to tripped rollover conditions which represent the most critical aspect for off-road and autonomous driving applications where the disturbance from the external environment play an important role [55]. Additionally, the artificial intelligence (AI) techniques are also proposed in literature [56,57] to obtain better anti-rollover prevention features. In particular, [58] has identified the Recurrent Neural Networks (RNN) as potential rollover risk-detector algorithm. However, by considering time and costs constraints [59], the big experimental data required to generate an efficient algorithm is not always possible.

Based on the widespread and in-depth bibliographic research found in the literature, to the knowledge of the authors, there are still many open questions on the rollover prevention assessment that the present activity attempts to solve. This paper aims to cover some relevant aspects in the the detection and prevention of the incipient rollover, with the following contributions:

• The numerical assessment of a real-time algorithm able to predict in advance the time to reach a specific load transfer ratio considered as the incipient rollover limit.

It is also important to note that one of the most important advantages related to the methodology described in the paper is the versatility of application to different powertrain and vehicle architectures, by providing an algorithm that requires the typical measurements available onboard the vehicle.

The manuscript is organised as follows: Section 2 describes the set of equations used for modelling the non-linear vehicle roll dynamics; Section 3 presents the two algorithms adopted for the load transfer ratio estimation; the the phase plane analysis of load transfer characteristics is carried out in Section 4; the ISO-LTR predictive time is introduced in Section 5 and its simulation assessment is shown in Section 6; finally, some conclusions are drawn in Section 7.

#### **2. Vehicle Model Description**

The evaluation of the rollover condition, together with the definition of a predictive index to prevent the incipient risk, is firstly carried out by investigating the influence of the most relevant quantities on the load transfer dynamics. A straightforward methodology is proposed with a three degrees-of-freedom (dofs) vehicle model to obtain an analytical and a straightforward approach for analysing the vehicle roll dynamics on flat and banked roads. The 3-dofs of the vehicle model, whose scheme is reported in Figure 1, are the vertical and the lateral vehicle motions (with respect to the road plane) and the relative roll motion between the sprung and the unsprung masses. The virtual roll axis is *R*, *Gs* is the CoG of the sprung mass and *Gui* is the CoG of the *ith* = *FL*, *FR*, *RL*, *RR* unsprung mass. The hypothesis behind the model are:

From Figure 1a, the following set of equilibrium equations for the whole vehicle is drawn :

$$\begin{cases} F\_{z\_R} + F\_{z\_L} - m\_s a\_{z\_s} - \sum\_i \left( m\_{\bar{u}\_i} a\_{z\_{\bar{u}\_i}} \right) - m\_\mathcal{S} \cos \phi\_\mathcal{R} = 0 \\ F\_{y\_R} + F\_{y\_L} - m\_s a\_{y\_s} - \sum\_i \left( m\_{\bar{u}\_i} a\_{y\_{\bar{u}\_i}} \right) - m\_\mathcal{S} \sin \phi\_\mathcal{R} = 0 \\ \left( F\_{z\_R} - F\_{z\_L} \right) \frac{T}{2} - m\_\mathcal{S} a\_{y\_s} (h\_R + h\_\mathcal{s} \cos \phi) - m\_\mathcal{S} \zeta (h\_R \sin \phi\_\mathcal{R} + h\_\mathcal{s} \sin \left( \phi\_\mathcal{R} + \phi \right)) + \\ -m\_\mathcal{S} a\_{z\_s} h\_\mathcal{s} \sin \phi - \sum\_i \left( m\_{\bar{u}\_i} a\_{y\_{\bar{u}\_i}} h\_\mathcal{u} \right) - m\_\mathcal{u} g h\_\mathcal{u} \sin \phi\_\mathcal{R} - G - I\_\mathcal{s} (\bar{\phi}\_\mathcal{R} + \bar{\phi}) = 0 \end{cases} (1)$$

where *FzR* = *FzFR* + *FzRR* and *FzL* = *FzFL* + *FzRL* are the total vertical forces on the right and left tyres, respectively. *FyR* = *FyFR* + *FyRR* and *FyL* = *FyFL* + *FyRL* are the right and left side total lateral forces, respectively. *ms* is the sprung mass, *mui* is the *ith* = *FL*, *FR*, *RL*, *RR* lumped unsprung mass (*mu* = ∑*i mui* ) and *m* = *ms* + *mu* is the vehicle total mass. *azs* and *azui* are the vertical components of the sprung and *ith* unsprung mass, respectively. *ays* and *ayui* are the lateral components of the sprung and *ith* unsprung mass, respectivelly. Moreover, the two angle *φR* and *φ* are the bank angle and the roll angle of the sprung mass w.r.t. the unsprung mass. The absolute sprung mass roll angle is *φA* = *φR* + *φ*. *hu* and *hR* are the heights of the unsprung mass and roll centre from the road, respectively, and *hs* is the distance between the sprung mass CoG from the roll centre. *g* is the gravity acceleration, *T* is the vehicle track width and *Is* is the roll sprung mass moment of inertia. *Rx* and *Ry* represent the horizontal and vertical components of the internal reaction in the virtual roll axis, respectively.

**Figure 1.** Free body diagrams of: whole vehicle (**a**), sprung and unsprung masses (**b**).

The term *G* groups together the relative vertical dynamics between the left and right unsprung masses and it is defined as reported in Equation (2):

$$G = \frac{T}{2} (m\_{u\_{\overline{R}\overline{R}}} a\_{z\_{\overline{R}}} + m\_{u\_{\overline{L}\overline{L}}} a\_{z\_{\overline{L}\overline{L}}} - m\_{u\_{\overline{R}\overline{R}}} a\_{z\_{\overline{R}\overline{R}}} - m\_{u\_{\overline{R}\overline{L}}} a\_{z\_{\overline{R}\overline{L}}}) \tag{2}$$

The roll equilibrium equation for the vehicle unspung mass from Figure 1b is:

$$\begin{aligned} \left( (F\_{z\_R} - F\_{z\_L}) \frac{T}{2} = \\ &= \left( F\_{y\_R} + F\_{y\_L} \right) h\_R - \left[ \sum\_i \left( m\_{u\_i} a\_{y\_{u\_i}} \right) + m\_{u} \mathfrak{g} \sin \phi\_R \right] \left( h\_R - h\_{u} \right) + K \phi + \mathcal{C} \phi - G \end{aligned} \tag{3}$$

where *K* and *C* are the roll stiffness and damping, respectively. Finally, roll equilibrium equation for the vehicle sprung mass from Figure 1b is:

$$\mathbf{I}\_s(\Phi\_\mathcal{R} + \Phi) + \mathbf{K}\phi + \mathbf{C}\phi = m\_s a\_{y\_s} h\_\mathfrak{s} \cos\phi + m\_s a\_{z\_\mathfrak{s}} h\_\mathfrak{s} \sin\phi + m\_\mathfrak{s} g h\_\mathfrak{s} \sin\phi \tag{4}$$

The heavy vehicle analysed within the present paper is defined by the parameters listed in Table 1.

**Table 1.** Vehicle parameters.

#### **3. Load Transfer Ratio Estimation**

The most accurate, and widely disseminated, index for detecting the incipient rollover risk is the Load Transfer Ratio (*LTR*) defined as the the relative vertical force on tires between the right and left sides of a vehicle:

$$LTR = \frac{F\_{z\_R} - F\_{z\_L}}{F\_{z\_R} + F\_{z\_L}} \tag{5}$$

When the *LTR* is equal to 0, the vertical load is equally distributed between the two sides, by means that the vehicle is far away from the rollover condition. The vehicle is very likely to rollover when the wheels on one side lift off the ground, condition verified when *LTR* = ±1.

However, the experimental acquisition of *LTR* is not feasible to realize since there are no sensors able to provide a direct measurement of the vertical tire load. For this reason, many authors have proposed multiple solutions to provide an accurate estimation of the *LTR* from the measurements commonly available onboard passenger cars, e.g., accelerations, linear and angular velocities, roll and pitch angles. Simple static [34,35], and dynamic indices have been used, in different works analysed in [24], to detect rollover risk during dynamic manoeuvres. The *LTR* proposed in [60] only considers the lateral

acceleration and the roll angle, while [61,62] neglects the unsprung masses, by placing the roll centre very close to the road plane at *hR* ≈ 0 and assuming the vehicle CoG height as the effective arm for lateral acceleration. Other proposed *LTR* indices are based on suspension parameters [50,63], by considering the effective torsional stiffness, torsional damping, roll angle and roll rate, and subsequently revised in [64] with the inclusion of the sprung mass lateral acceleration. RIs based on the tyre deflection is proposed in [65,66], also adapted for heavy vehicles. A more complex index, which considers the rolling motion, the lateral acceleration, and the time to wheel lift (TTWL) is proposed in [38], where the TTWL is achieved with a phase plane analysis of roll angle and roll rate. Although a RI for banked roads is already developed in [67], the road bank angle influence receives a lower attention from the literature since its effect is not as relevant as for off-road designed vehicles. Finally the index proposed in [22] tries to to consider both tripped and untripped rollovers, by developing a four dofs vehicle model.

However, most of the above mentioned formulations does not include the influence of road global, i.e., bank angle, and local, i.e., speed bumps or potholes, irregularities on the load transfer among the four wheels, which represents a fundamental aspect for the off-road driving scenario. The present section of the paper describes a model-based and a statistic formulation for the *LTR* estimation. The first approach aims at providing an analytical methodology to express the *LTR* as function of the vehicle states and road profile, able to predict their influences on the *LTR* dynamics. However, the accuracy of the model-based estimation is drastically affected by the model parameters uncertainties, a drawback that is improved by proposing an alternative approach based on the Recurrent Neural-Network (RNN) theory, thus allowing an improved estimation of the *LTR* directly from conventional measurements.

#### *3.1. Model-Based Estimation*

The *LTR* model-based estimation is obtained from the 3-dofs model presented in the previous section. By combining the roll equilibrium equations of Equations (1) and (4) and by considering the vertical equilibrium equation in Equation (1), the LTR is then calculated as follows:

$$LTR\_{\rm cell} = \frac{2}{T} \frac{K\phi + \mathcal{C}\dot{\phi} + m\_{\rm s}a\_{\rm y\_s}h\_{\rm R} + \sum\_{i} \left( m\_{\rm li\_i}a\_{\rm y\_{u\_i}}h\_{\rm u} \right) + \left( m\_{\rm s}gh\_{\rm R} + m\_{\rm u}gh\_{\rm u} \right) \sin \phi\_{\rm R} - G}{mg \cos \phi\_{\rm R} + m\_{\rm s}a\_{\rm z\_s} + \sum\_{i} \left( m\_{\rm u\_i}a\_{\rm z\_{u\_i}} \right)} \tag{6}$$

For a direct comparison against the *LTR* formulation commonly available in the literature, the particular cases of *a*) negligible unsprung dynamics and vertical accelerations (*mui* = 0 and *azs* = *azi* = 0) and *b*) negligible unsprung dynamics and vertical accelerations on a flat road (*mui* = 0, *azs* = *azi* = 0 and *φR* = 0) are also reported in Equation (7) and in Equation (8):

$$LTR\_{est,spring} = \frac{2}{T} \frac{K\phi + \mathcal{C}\dot{\phi} + m\_s a\_{y\_s} h\_R + m\_s g h\_R \sin\phi\_R}{mg \cos\phi\_R} \tag{7}$$

$$LTR\_{\rm cst,flat} = \frac{2}{T} \frac{K\phi + \mathcal{C}\dot{\phi} + m\_s a\_{y\_s} h\_R}{mg} \tag{8}$$

Three driving scenarios, shown in Figure 2, are built in IPG CarMaker® to simulate the realistic behavior of a heavy duty vehicle and to compare the efficacy and the reliability of the *LTR* formulations in Equations (6)–(8):

height amplitude of 10 cm) under the right vehicle side, called asymmetrical waves road in the rest of the paper.

**Figure 2.** Driving scenarios created in IPG CarMaker® for the *LTR* estimation: double lane change manoeuvre on a flat road (**a**) , straight manoeuvre on a banked road (**b**) and straight manoeuvre on an asymmetrical waves road (**c**).

The first scenario is chosen to evaluate the influence of the maneuver severity, imposed by the driver behavior, on the *LTR* dynamics. The simulation results are reported in Figure 3, in terms of *LTR* estimation, driver commands, vehicle speeds, accelerations and roll angles.

Figure 3 shows that the manoeuvre is aggressive enough to push multiple times the *LTR* towards the rollover condition. It is clear that the lateral acceleration, imposed by the severity of the test and by the dynamic response of the driver, represents the main cause for the marked transfer loads from one vehicle side to the other one. Indeed, the road is flat and does not influence the roll dynamics meanwhile the unsprung mass only slightly affect the *LTR* dynamics, as shown by the comparison of *LTRest* against *LTRest*,*sprung* and *LTRest*, *flat*, due to the low frequency content of steering wheel angle imposed by the driver.

The second scenario aims at investigating the effect of the road bank angle on the vehicle transfer loads, as shown in the simulation results of Figure 4.

**Figure 3.** Double lane change manoeuvre at initial speed of 100 km/h on a flat road. From the top left side to the down right side: *LTR*, steering wheel angle and vehicle speed, vertical accelerations, lateral accelerations, roll angles and angular velocities.

The straight manoeuvre is selected to minimize the influence of lateral accelerations and to focus the attention on the bank angle effect on the load transfer dynamics. The road is banked up to 30 deg that represents a suitable threshold for approaching the rollover limit, identified by the condition *LTR* = 1. The general formulation of *LTRest* is able to provide a perfect estimation of the actual *LTR*, meanwhile the basic formulation *LTRest*, *flat* completely underestimate the load transfer behavior since it is valid only for flat roads. The estimation *LTRest*,*sprung* is closer to the actual *LTR* behavior than the *LTRest*, *flat*, but it shows a lower accuracy at high bank angles since the unsprung mass still play an important role in the roll dynamics, although *mu* represents only the 16% of the total vehicle mass. Finally, the last driving scenario is designed to evaluate the influence of the unsprung mass on the *LTR* estimation, by exciting the vehicle roll dynamics with a sinusoidal road profile applied only to the right side of the vehicle. The manoeuvre on the asymmetrical waves road, with a wavelength equal to the vehicle wheelbase, is run at the constant speed of 40 km/h, in order to excite the higher frequency content of the unsprung mass vertical dynamics. The simulation results are shown in Figure 5.

**Figure 4.** Straight manoeuvre at 30 km/h on a banked road. From the top left side to the down right side: *LTR*, steering wheel angle and vehicle speed, vertical accelerations, lateral accelerations, roll angles and angular velocities.

This manoeuvre can also evaluate the effect of a local road disturbances, i.e., road bump, pothole, etc. ..., which are typical conditions for the off-road applications and it has a different impact on the roll dynamics with respect to a global road geometry such us the presence of a bank angle. Indeed, the roll dynamics derived from the asymmetrical excitation induces a persistent condition of incipient rollover: during the small time frame of 1 s, the tires of the right vehicle side lift off three times. In this contest, the unsprung mass influence becomes much more important with respect the previous driving manoeuvres, and it is appreciable how the generic model-based estimation *LTRest* is still reliable in detecting the *LTR* dynamics, if compared against the other two formulations.

**Figure 5.** Straight manoeuvre at 40 km/h on an asymmetrical waves road. From the top left side to the down right side: *LTR*, steering wheel angle and vehicle speed, vertical accelerations, lateral accelerations, roll angles and angular velocities.

#### *3.2. Recurrent Neural Network Estimation*

As previously anticipated in Section 3, RNNs can be used to improve the estimation of rollover indicators, such as *LTR* and roll angle. The algorithm presented differentiates from the literature in terms of architecture, standardization method, inputs and outputs selection. In particular, the outputs are not defined as discretized quantities (higher or lower rollover risk factors), but they provide directly the estimation of the *LTR*. The RNNs method is chosen by considering the nature of the problem, since they belong to the category of artificial neural network designed to recognize patterns from a set of time histories.

The methodology adopted is characterized by two main phases. The first one consist in the generation of a large time domain dataset which is needed to successfully train the neural networks. Due to the big amount of required data, the time histories are generated by simulating the behaviour of the vehicle in the IPG CarMaker® environment. In particular, a set of 15 maneuvers (and related scenarios), are selected to induce the rollover risk. Each manoeuvre is 500 s long, to impose the same weight for the RNN training. In total, a starting dataset of 750, 000 points is obtained for each variable of interest, sampled at

constant frequency of 100 Hz. The variables obtained from IPG CarMaker® used for the neural networks inputs and outputs are shown in Table 2.

**Table 2.** RNN: inputs and outputs.

The approach developed suggests an integration between the neural network algorithm with the simulation environment (IPG CarMaker®), which represents its data provider. This leads to the following practical implications:

The second phase consisted in the training and testing of the AI algorithms, articulating the process in several steps with the Deep Learning Toolbox available in Matlab®. Firstly, a data pre-process and normalisation elaboration is carried out to obtain dimensionless quantities. The processed data is then adopted to elaborate the architecture of the RNN, by testing multiple solutions, leading to the final solution reported in Table 3.

**Table 3.** RNN: Architecture.

The features indicated in Table 3, corresponds to the inputs of the RNN, and so the responses for the outputs. In particular, regarding the RNN architecture selection, the number of hidden units (the neurons) is chosen through empirical rules. It has been observed that with more than 100 hidden units, the results do not improve in accuracy and the process becomes excessively time consuming, while a smaller number of neurons inevitably leads to a decline in terms of performance (measured by RMSE and loss function). The number and the typology of layers have been chosen following a similar criteria.

Given the RNN architecture, the algorithm training is ready to start by obtaining a neural network that, provided with the 8 inputs (Table 2), is able to estimate the 3 outputs with good precision and immediate response. This is also validated during the testing phase, where the features are known, and the vehicle responses not. The results related to three different maneuvers, generated in the IPG CarMaker® environment, are presented.

The first maneuver analyzed simulates the vehicle driving over the Bernina Pass. This is a completely unknown road to the RNN, and the results obtained, shown in Figure 6a,b are definitely promising in terms of accuracy (see the resultant Mean Absolute Error and Mean Square Error in Table 4). Moreover, the RNN estimates are collected in the lapse of 0.01–0.05 s. For the sake of brevity, results related to the roll angle have been omitted.

**Figure 6.** Estimation of the vehicle Load Transfer Ratio at front (**a**) and rear (**b**) axles driving over the Bernina Pass.

The second maneuver presented is obtained from one of the 15 database composing the training. The analysis aims at understanding if similar maneuvers to the ones already adopted for the RNN training, lead to an accuracy worsening, since the weights introduced in the RNN are influenced by the maneuvers set. Results in Figure 7a,b show how the RNN keep their high level of performance. There are only minimal differences that can be considered negligible for the required accuracy.

The Fishook manoeuvre, whose estimating results are shown in Figure 8a,b, produces the less promising results (as confirmed by Table 4). The reasons behind that is that the manoeuvre time length and speed profile are lower w.r.t. to the other manoeuvres introduced in the data set.

To provide a quantitative information related to the accuracy of the neural networks, the Mean Absolute Error (MAE) and the Mean Square Error (MSE) are reported for each manouvre in Table 4.

**Table 4.** Results: Mean Absolute Error and Mean Square Error related to the tested maneuvers.

**Table 4.** *Cont.*

**Figure 8.** Estimation of the vehicle Load Transfer Ratio at front (**a**) and rear (**b**) axles during a FishHook maneuver.

#### **4. ISO-LTR Phase Plane Portrait**

As mentioned in the previous section, the model-based *LTR* formulation, expressed in its more robust and generic form by the Equation (6), is not suitable for a realistic implementation of the LTR estimation but rather it represents a valid tool for the analytical correlation between the vehicle states and the *LTR* dynamics. One of the most adopted tool to evaluate the vehicle states dynamics is the phase plane analysis. By focusing the attention only to the relative roll dynamics between the sprung and the unsprung masses, the vehicle states selected for plotting the phase plane map are the relative roll angle *φ* and roll rate *φ*˙ between the sprung and the unsprung masses. The phase plane analysis has the fundamental advantage to describe the trajectory of the vehicle states in a unique plot that can also be obtained with complex and non-linear models or directly from experimental measurements. In this paper, the phase-plane plot is drawn by running a fast ramp steering manoeuvre at constant speed in the IPG CarMaker® simulation environment. According the NHTSA's standard (49 CFR Part 575), a steering wheel rate of 720 deg/s is imposed, and the test is repeated at different vehicle speeds with step of 10 km/h. An example of simulation results obtained during a fast ramp steering manoeuvre at 70 km/h on a 10 deg banked road is shown in Figure 9.

**Figure 9.** Fast ramp steering manoeuvre at 70 km/h on a 10 deg banked road. From the top left side to the down right side: *LTR*, steering wheel angle and vehicle speed, vertical accelerations, lateral accelerations, roll angles and angular velocities.

This aggressive manoeuvre is able to depict the transient behavior of the vehicle roll dynamics at a constant vehicle speed and road bank angle, through a trajectory line in the phase plane plot. The manoeuvre is repeated for different vehicle speeds, thus obtained the phase plane portrait for a predefined value of the road bank angle. An example of the phase plane portrait for a a 10 deg banked road is reported in Figure 10.

**Figure 10.** ISO-LTR characteristics obtained during a fast ramp steering manoeuvre on a 10 deg banked road at different speeds.

Each trajectory in the phase plane portrait represents the *φ*˙ − *φ* correlation during a fast ramp steering manoeuvre at constant vehicle speed and include the steady-state and the transient responses of the vehicle roll dynamics. Each trajectory shows a colour variation as function of the corresponding *LTR* value, passing from a lower rollover risk represented by a blue markers (*LTR* = 0) to a high level of rollover risk with the red markers (*LTR* = 1) indicating that the left wheels lifted off.

When the phase plane portrait (*φ*˙, *φ*) is built and defined for a specific boundary condition, i.e., the road bank angle, the ISO-LTR characteristics, i.e., the locus of constant vertical load transfers (*LTR* = *q*, with *q* = 0.4, 0.5, ... , 0.9), can be drawn, as shown in Figure 10 for *φR* = 10 deg. The ISO-LTR characteristics map is characterized by parallel lines, each corresponding to a constant *LTR* level.

By combining the vehicle roll equilibrium equation in Equation (1) with the roll equilibrium equation of the sprung mass in Equation (4), the following analytical expression of the linear ISO-LTR characteristics in the phase-pane is obtained:

$$\begin{split} \dot{\Phi} = -\frac{\mathcal{K}}{\mathbb{C}} \Phi - \left[ \frac{m\_{u} a\_{g\_{u}} h\_{R} + \sum\_{i} (m\_{u\_{i}} a\_{g\_{u\_{i}}} h\_{u})}{\mathbb{C}} \right] - \left[ \frac{q \frac{\mathcal{T}}{2} m\_{s} a\_{z\_{s}} + q \frac{\mathcal{T}}{2} \sum\_{i} (m\_{u\_{i}} a\_{u\_{i}}) - G}{\mathbb{C}} \right] + \\ + \left[ \frac{q \frac{\mathcal{T}}{2} m\_{s} \zeta \cos \phi\_{R} - (m\_{s} g\_{R} + m\_{u\_{i}} g\_{u}) \sin \phi\_{R}}{\mathbb{C}} \right] \end{split} \tag{9}$$

The slope of the ISO-LTR linear characteristics is always negative and it depends on the vehicle suspension system parameters through the total roll stiffness *K* and roll damping *C* coefficients. In particular, an increment of the total roll stiffness or, equivalently, a reduction of the total roll damping, would lead a lower ISO-LTR lines sensitivity to the relative roll rate *φ*˙.

It is possible to express the ISO-LTR characteristics in the phase plane through the following linear equation:

$$
\dot{\phi}\_{LTR=\eta} = k \phi\_{LTR=\eta} + n \tag{10}
$$

where *k* = −*K C* and *n* = *f ays*, *ayui* , *azs*, *azui* ,*φR* are the slope and the x intercept of the ISO-LTR corresponding to *LTR* = *q*, respectively.

To provide a proof of this conclusion, the influence of the the road bank angle on the phase-plane portrait and the distribution of the ISO-LTR lines are shown in Figure 11.

**Figure 11.** Comparison of the ISO-LTR characteristics for a flat and a 10 deg banked road.

The figure demonstrates that the ISO-LTR lines keep the same negative slope, not influenced by the road bank angle, whose main effect is a shift of the ISO-LTR characteristics towards smaller relative roll angles *φ*.

#### **5. ISO-LTR Predictive Time**

The previous sections present the mathematical background to detect the incipient rollover risk through the evaluation of the *LTR*, its representation on the phase plane plot through the definition of the ISO-LTR characteristics and their perturbations due to the driver behaviour, the local and global road profile geometry. The present section aims at providing a predictive indication of the time the vehicle requires to cross one of the ISO-LTR lines, defined as the ISO-LTR Predictive Time (*ILPT*) in the rest of the paper.

The principle underlying the procedure, based on the formulation proposed by [53], is described by Figure 12 in the phase plane plot.

**Figure 12.** Schematic principle of the *ILPT* definition.

The figure shows three trajectories corresponding to different manoeuvres represented in the phase plane. After a generic time instant, the three trajectories reach the points *P*, *PI* , *PI I* approaching to the same roll angle but with different roll rates. The condition of incipient rollover is identified by the blue ISO-LTR line. The *ILPT* is defined as the time required to reach the intersection point between the ISO-LTR line and the tangent line at the current point of the phase plane trajectories. For example, the tangent line in *P* is expressed by:

$$
\dot{\Phi} = k\_P(\phi - \phi\_P) + \dot{\phi}\_P \tag{11}
$$

where *φ*¨*P*, *φ*˙ *P* and *φP* represent the relative roll acceleration, roll rate and roll angle corresponding to the point *P* of the phase plane trajectory and *kP* = *φ*¨*P φ*˙ *P* is the slope of the tangent line. The intersection point between the tangent line and the ISO-LTR line *Q* is then calculated as follows:

$$\begin{cases} \phi\_Q = \frac{k\rho\dot{\phi}^p - \dot{\phi}^p + n}{k\_P - k} \\ \dot{\phi}\_Q = \frac{k k \rho\dot{\phi}^p - k \dot{\phi}^p + k \rho n}{k\_P - k} \end{cases} \tag{12}$$

Hence, the *ILPT* represents the time required to reach point *Q* from point *P*, under the hypothesis to follow the tangent line expressed by Equation (11), as also described by [53]:

$$ILPT = \frac{\left[1 - \text{sign}(z)\right]}{2} \sqrt{\frac{\left(\phi\_Q - \phi\_P\right)^2 + \left(\phi\_Q - \phi\_P\right)^2}{\dot{\phi}\_P^2 + \ddot{\phi}\_P^2}}\tag{13}$$

with:

$$z = (\dot{\phi}\_P - k\phi\_P - n)(-\dot{\phi}\_P + k\phi\_P - n) \tag{14}$$

where *z* is introduced to saturate the *ILPT* to 0 when the point *P* crosses the corresponding ISO-LTR line. The calculation of the *ILPT* requires the knowledge of the current relative roll acceleration, roll rate and roll angle and the parameters *k* and *n* of the selected ISO-LTR line. The *ILPT* provides a time information only if it exceeds a predefined upper thresholds, i.e., 0.5 s above which the vehicle is far away from the rollover risk.

#### **6. Simulation Results**

The predictive capabilities of the *ILPT* are verified in IPG CarMaker® environment through the simulation of the following driving manoeuvres:

The simulation results are analyzed both in the phase plane and in the time domains. The phase-plane analysis provides a graphical interpretation for the *ILPT* calculation. On the other hand, the time domain analysis is carried out to asses the methodology accuracy in predicting the critical condition. Indeed, the reference time *TRef* , i.e., the time required to reach a defined *LTR* threshold *LTRth*, is backward evaluated from the whole post-processed time history and it is compared against the *ILPT* estimation. *TRef* is represented in the time domain as −45 deg inclined lines:

$$T\_{Ref} = -t + T\_{LTR=LTR\_{th}} \tag{15}$$

where *TLTR*=*LTRth* is the time that *LTR* = *LTRth*, and it is detected by a backward analysis of the *LTR* time history. The *TRef* represents a key quantity to asses the performance of any predictive algorithm to detect a particular event, i.e., the incipient rollover, but it can be only calculated when the whole time history up to the event occurrence is available and it is not suitable for a real-time implementation. The *LTR* threshold, here considered as potentially critical to provoke the vehicle rollover, is assumed to be ±0.8. The smaller the deviation between *ILPT* and *TRef* , the better the predictive capabilities of the proposed methodology.

#### *6.1. Fast Ramp Steering on a Banked Road*

The first manoeuvre is one from the simulation set adopted in the previous section to draw the ISO-LTR characteristics of Figure 10.

The resultant phase plane trajectory of the manoeuvre is reported in the topside of the Figure 13, where the graphical solution for the *ILPT* calculation is provided for four time instants. The resultant *ILPT*, elaborated in real-time during the manoeuvre, is plotted in the downside of Figure 13 where the *LTR* time history is also compared against the threshold *LTRth*.

**Figure 13.** Fast ramp steering manoeuvre at 30 km/h on a 10 deg banked road: phase plane trajectory (top), *LTR*, reference (*TRef* ) and estimated (*ILTP*) time to reach *LTR* = 0.8 (down).

At the beginning of the manoeuvre, when the *LTR* level indicates a substantial distance from the rollover risk, the *ILPT* (red circles) deviates from the *TRef* (red dashed line) to reach the condition of *LTR* = 0.8, event that occurs after 0.24 s from the application of the ramp steering. The time deviations between the estimated *ILPT* and the *TRef*

suddenly drops down after 0.05 s, when the *ILPT* starts tracking the reference time until the occurrence of the critical threshold *LTR* = 0.8. This deviation can be graphically explained by the phase-plane portrait in Figure 13. At the beginning of the ramp steering (left side of the trajectory), the trajectory strongly deviates from the tangent line (dashed lines). As soon as the trajectory is approaching the ISO-LTR line (blue continuos line), the tangent line begins to better approximate the future phase-plane trajectory thus allowing the convergence of the *ILPT* towards the *TRef* . It is also interesting to note that the *ILPT* slightly underestimates the *TRef* when approaching the corresponding ISO-LTR line, thus providing a more conservative estimate.

#### *6.2. Double Lane Change on a Flat Road*

The methodology is also applied to a second manoeuvre that is not used to extrapolate the ISO-LTR characteristics or any other results shown in the previous sections. This is the double lane change test which represent a critical manoeuvre from the rollover point of view. The time domain plot of the *LTR*, the *ILPT* and reference time are shown in Figure 14.

**Figure 14.** Double lane change manoeuvre at initial speed of 100 km/h on a flat road: *LTR*, reference (*TRef* ) and estimated (*ILTP*) time to reach *LTR* = 0.8.

During the double change manoeuvre, the vehicle crosses four times the *LTR* threshold. In all cases, the methodology is able to detect in advance the critical point as highlighted by the *ILTP* that starts dropping to zero when the vehicle is approaching the *LTR* threshold. Even during this maneuver, the *ILPT* tends to underestimate the reference time when the vehicle is getting closer to the ISO-LTR boundary, thus providing a conservative and safer prediction. During the last part of the maneuver, the algorithm provides a false positive indication of dangerous situation since the *ILTP* drops to 0.3 s even if the *LTR* did not reach the threshold of 0.8 after this point. This behavior is explained by the sudden reduction of the *LTR* at *t* = 30 s, which is successfully recovered by the stabilizing intervention of the driver.

Even though the vehicle does not rollover during the maneuver, the *ILPT* can predict all the potential dangerous situations in a time frame of a few tenths of seconds, which is suitable for an active control logic to intervene.

#### **7. Conclusions**

The activity presented in this paper aims at providing some usefull tools and methodologies to detect and prevent the incipient risk of rollover, applied to heavy vehicles designed for off-road driving applications. The main conclusions drawn from the results and the methodologies presented can be summarized as follows:

**Author Contributions:** Conceptualization, A.T., L.D., G.P. and M.V.; methodology, A.T., L.D. and M.V.; software, A.T. and F.V.; validation, A.T., L.D. and F.V.; formal analysis, A.T., L.D. and M.V.; investigation, A.T., L.D., G.P. and M.V.; writing—original draft preparation, A.T., L.D. and F.V.; writing—review and editing, A.T., L.D., F.V., G.P. and M.V.; visualization, L.D.; supervision, A.T., L.D. and M.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**

## *Article* **Discriminant Analysis of the Vibrational Behavior of a Gas Micro-Turbine as a Function of Fuel †**

**Vincenzo Niola \*, Sergio Savino, Giuseppe Quaremba, Chiara Cosenza, Armando Nicolella and Mario Spirto**

Department of Industrial Engineering, University of Naples "Federico II", Via Claudio 21, 80125 Napoli, Italy **\*** Correspondence<EMAIL_ADDRESS>
† This paper is an extended version of "Spirto, M.; Savino, S.; Quaremba, G.; Cosenza, C.; Nicolella, A.; Niola, V. Discriminant analysis of the vibrational dynamics of a gas microturbine with different fuels. In Proceedings of the 5th Jc-IFToMM International Symposium, the 28th Jc-IFToMM Symposium on Theory of Mechanism and Machines, Kyoto, Japan, 16 July 2022".

**Abstract:** Several studies were conducted previously on fuel and biofuel performance of microturbines. The present paper combines experimental and statistical approaches to study the vibrational behavior of a gas micro-turbine supplied with different pure fuels and admixed with rapeseed oils. Experimental tests carried out at different operating conditions have allowed us to build a classification model through using discriminant analysis. The classification model can distinguish the vibrational behavior occurring when the turbine is fueled with kerosene, or pure and admixed diesel with rapeseed oil. Moreover, the methodology has even allowed us to highlight differences in vibrational behavior caused by small amounts of rapeseed oil admixed in the fuel. The model reliability, in terms of Cohen's kappa, results in optimal data classification.

**Keywords:** gas micro-turbine; discriminant analysis; vibrational analysis; diesel blends; kerosene; biofuel; statistical index; accelerometer; confusion matrix; Cohen's kappa

#### **1. Introduction**

In recent years, significant effort has been made to lessen the environmental impact of operating gas turbines [1–3].

Micro gas turbines have the potential source to be an alternative power unit (APU) in range-extended electric vehicles (REEVs) [4].

Due to their improved thermal efficiency and reduced fuel consumption, diesel engines make up a significant portion of the passenger car industry in Europe. However, due to increasingly strict exhaust emission laws, particularly for NOx and PM, a more complex and expensive gas exhaust treatment system is required for most diesel vehicles, which also results in increased fuel usage. By 2030, several EU regulations and incentives would call for up to 45% bio-components in fossil fuel. It has been demonstrated that biofuels in the form of pure plant oils cannot be used for automotive applications in diesel engines efficiently because of technical issues brought on by their considerably higher viscosity, corrosive character, and increased exhaust smoke emissions [5].

Biofuels are, broadly speaking, any fuel source made from organic material, such as firewood, charcoal, animal fats and oils, dung, and vegetable oils. To mimic the performance and physical properties of fossil fuels, biofuels must be modified. The development of new power and propulsion technologies will necessitate a better control of the harmful chemicals generated by combustion sources. As energy demand increases, concerns about global warming could result in restrictions on greenhouse gas emissions from fossil fuel facilities. This has prompted substantial study into carbon capture and storage, as well as a growth in renewable energy sources. Many investigations of fuel preparation and emission characterizations of systems powered entirely or partially by pure vegetable oils are available [6–13].

**Citation:** Niola, V.; Savino, S.; Quaremba, G.; Cosenza, C.; Nicolella, A.; Spirto, M. Discriminant Analysis of the Vibrational Behavior of a Gas Micro-Turbine as a Function of Fuel. *Machines* **2022**, *10*, 925. https:// doi.org/10.3390/machines10100925

Academic Editors: Marco Ceccarelli, Giuseppe Carbone and Alessandro Gasparetto

Received: 27 September 2022 Accepted: 8 October 2022 Published: 12 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

A crucial aspect is the relation between spray quality and combustion performance in micro gas turbines burned with biofuels and biomasses [14,15].

Pure vegetable oils might cause unwanted vibrations by damaging the injection system or combustion device. The ideal answer would be to have a tool that could foresee and/or track the vibrational state of the combustion device in real time. By using appropriate quantifiers generated from microphone and accelerometer inputs, investigations based on acoustic and vibrational measurements seem to offer an intriguing diagnostic and predictive answer [16–20]. Other authors have proposed a neural-network-based tool that allows them to ensure protection and the safety measures against the instability phenomena in a gas turbine based on the modeling of its dynamic behavior [21].

In this context, Allouis et al. [22] have evaluated the impact of biofuel properties on emissions and performances of a micro gas turbine using combustion vibrations detection. Other methodologies described in [23–25] are employed for detecting anomalies in mechanical systems such as gas micro-turbines.

The present paper aims to study the vibrational behavior of a low-emission gas microturbine for power generation, fueled with different liquid fuels, including commercial diesel oil and its blend with pure rapeseed oil. Kerosene is the design fuel for this turbine. The authors of this current paper have presented a preliminary work on this topic [26].

In the first part, the study describes the experimental phase in which vibrational signals have been acquired through accelerometers properly mounted on the injector. Afterward, signal processing has been performed through discriminant analysis. A classification model has been developed to distinguish the vibrational behavior associated to the turbine fueled with kerosene, pure diesel, and admixed diesel with rapeseed oil.

#### **2. Materials**

The study system (the same as of [22]) is a Capstone 30 model 18-blade gas microturbine (Figure 1) with a maximum power output of 30 kW. Since its output power is less than 100 kW, it can be included in the micro-turbine category.

**Figure 1.** Turbine side view.

The maximum fuel flow rate is 10 kg/h, and the exhaust gas temperature is approximately 590 ◦C, while at the system discharge it is 276 ◦C (as nominal temperature). The test bench of the system is reported in Figure 2.

**Figure 2.** System test bench.

Two different blends of diesel with rapeseed oil (respectively at 1 and 3% by volume) and kerosene (the turbine's designed fuel) are used as turbine fuels. The choice of such low percentages of rapeseed oil in the two additive fuels is necessary to characterize the machine's behavior from a vibrational point of view and to grasp the minimal variations resulting from the use of fuels with a very similar composition.

As for the choice of rotation speeds, speeds lower than the maximum one (90,000 rpm) were chosen to study the machine's behavior under conditions that gradually deviate from the maximum. Tests were conducted for turbine rotation speeds of 75,000, 80,000, and 85,000 rpm. A piezoelectric uniaxial accelerometer (PCB 352C22) and a data acquisition system (LMS SCADAS Mobile SM01) were used to sense and acquire the turbine vibrations. The accelerometer was mounted on a rigid bracket in-built with the injector. The need to use the bracket is due to the high temperature reached by the turbine: in the injector, the temperature can reach a high value to damage the accelerometer.

The accelerometer signals were acquired for a duration of 10 s with a sampling frequency of 102,400 Hz. As a preliminary step, the accelerometer was calibrated using a Bruel & Kjaer Calibration Exciter Type 4294. Below is the procedure followed to carry out the experimental tests at 85,000 rpm (the procedure is similar for the other two speeds):

Table 1 shows a summary of the acquisition data.

**Table 1.** Tests summary.

#### **3. Methods**

The acquired vibrational signals were analyzed by means of discriminant analysis of statistical indices directly calculated on the raw signals.

The method used for the creation of the classificatory model obtained through discriminant analysis is stepwise, i.e., the variables are not inserted into the model at the same time but in steps, i.e., only those variables with the greater discriminating weight are inserted. At each step, the variable with the lowest Wilks lambda value and, in turn, the highest F coefficient value is entered, i.e., the variable that contributes to better differentiate the groups. More details on this type of analysis can be found in the bibliographical references [26,27].

The chosen statistical indices for the analysis are asymmetry, kurtosis, shape factor, quadratic oscillation index, root mean square value, crest factor, non-normalized Shannon entropy, logarithmic entropy, synchrony index, Pearson's correlation index, Kendall's correlation index, and Spearman's correlation index [27–30]. The statistical indices used to create the three discriminating models (Section 4) will be briefly described. In the description, *n* is the length of the signal and *xi* is the *i*-th component of the signal. The quadratic oscillation index (*O*2) that provides indications of the oscillating phenomenon is described by Equation (1):

$$O^2 = \sqrt{\frac{1}{n-1} \sum\_{i=1}^{n-1} \left(\mathbf{x}\_i - \mathbf{x}\_{i+1}\right)^2} \tag{1}$$

The non-normalized Shannon entropy (*SE*) that provides a measure of the degree of order or disorder of a signal is described by Equation (2):

$$SE = -\sum\_{i=1}^{n} \mathbf{x}\_i^2 \log \left(\mathbf{x}\_i^2\right) \tag{2}$$

The asymmetry (*γ*), which is described in Equation (3), gives an indication of the distance between the mean value and the mode of the signal:

$$
\gamma = \frac{m\_3}{\sigma^3} \tag{3}
$$

where *m*3 is the third-order central moment of the signal, while *σ* is the standard deviation. The root mean square (*RMS*) value, which is described in Equation (4), gives an indication of the loads acting on the system and the system speed:

$$RMS = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \mathbf{x}\_i^2} \tag{4}$$

Finally, to assess the degree of accuracy and reliability of the classification, Cohen's kappa (*κ*) was calculated [31,32]. It is a statistical coefficient representing the degree of accuracy and reliability in a statistical classification defined from a confusion matrix, as reported in Equation (5):

$$\kappa = \frac{P\_0 - P\_{\mathfrak{e}}}{1 - P\_{\mathfrak{e}}} \tag{5}$$

*P*0 is the sum of the probabilities along the main diagonal of the confusion matrix, as given in Equation (6):

$$P\_0 = \sum\_{i=1}^{n} P\_{i,i} \tag{6}$$

where *Pi*,*i* is the probability of the generic element along the main diagonal of the confusion matrix. The *Pe* formula is given in Equation (7):

$$P\_{\varepsilon} = \sum\_{i=1}^{n} P\_{i,TOT} \cdot P\_{TOT,i} \tag{7}$$

where *Pi*,*TOT* is the total probability along the *i*-th row, and *PTOT*,*i* is the total probability along the *i*-th column of the confusion matrix.

#### **4. Results and Discussion**

The first operation carried out on the sampled signals was the calculation of the statistical indices indicated in paragraph 3. From the single sampled signal at a given speed of the micro-turbine powered with a given fuel, 125 values of the generic statistical index were obtained, i.e., one value was calculated every 100 revolutions of the turbine: this means that a total of 500 values of the generic index were obtained for each speed considered. An example-blocks diagram of the algorithm adopted to calculate the generic index *I* for the four signals (for each speed) was reported in Figure 3.

**Figure 3.** Algorithms implemented for calculating indices.

Subsequently, each of the calculated indices was used as a potential variable for the creation of the classification model for each speed. The software used to create the classification model returned only three variables (indices) out of all those inserted in the model. Table 2 summarizes the variables entered into the model, i.e., those variables that were able to explain 100% of the cumulative variance for each speed analyzed.

**Table 2.** Variables chosen for the model.

Starting from these variables, it is possible to write three discriminant functions, called *Fi*, for each analyzed speed:

$$F\_1 = a\_1 X\_1 \tag{8}$$

$$F\_2 = b\_1 X\_1 + b\_2 X\_2 \tag{9}$$

$$F\_3 = c\_1 X\_1 + c\_2 X\_2 + c\_3 X\_3 \tag{10}$$

where *Xi* are the variables (indices) chosen by the model and shown in Table 2, while *ai*, *bi*, and *ci* are the standardized coefficients associated with each function returned by the software.

The three discriminant functions (Table 3) are different for each turbine speed. The coefficients a, b, and c are different for each case study.

Table 4 displays the variance percentages contributed by the individual functions and the cumulative variance. For all three turbine speeds, it is possible to exceed 95% of the cumulative variance with only two functions. The third added function contributes a value of less than 5% to the cumulative variance, which is necessary to attain 100% of the total cumulative variance.

**Table 3.** Discriminant function.

#### **Table 4.** Variance percentage.

#### *4.1. Variables Included in the Classification Model*

In the following, the obtained indices for the fulfillment of the classification model at each speed (Table 2) are reported. All the following diagrams are shown as a function of the number of periods where the single period corresponds to 100 turbine revolutions.

Figure 4 shows the trends of the three indices (variables) inserted in the discriminant model for the speed of 75,000 rpm.

**Figure 4.** Indices for 75,000 rpm: (**a**) quadratic oscillation index, (**b**) non-normalized Shannon entropy, and (**c**) asymmetry.

In Figure 4a, which shows the trend of the quadratic oscillation index, the clear distinction between the four fuels can be seen immediately. Indeed, while kerosene and pure diesel are in the upper part of the diagram, the two admixed diesels clearly present lower values, with the diesel with 1% rapeseed oil averagely higher than the diesel with 3% rapeseed oil: this phenomenon can be seen in the greater tendency of diesel with 1% rapeseed oil towards the behavior of pure diesel. Since high values of this index emphasize that the system dynamic energy is better utilized, it can be said that using pure rather than admixed fuels makes the system better in usage use of dynamic energy. This can be explained by the fact that admixed diesel has an explosive power lower than the pure, reducing its energy output. This makes it possible to see that, by using biofuels, CO2 emissions and other pollutants can certainly be reduced at the expense of efficiency.

The distinction between pure and admixed fuels can also be seen in Figure 4b, in which the trend of non-normalized Shannon entropy is shown, by observing that admixed fuels exhibit higher values than pure fuels. This trend highlights how admixed fuels show little stationary and periodic vibrational cycles in contrast to pure fuels: this phenomenon is due to the reduction in explosive power caused by the presence of the oil and the resulting non-homogeneity of the mixture.

The asymmetry index is the last one adopted for modeling the 75,000 rpm case, shown in Figure 4c. By simply plotting the trend, it is not possible to show a difference between the four fuels used because of the overlapping of the curves and their highly chaotic course. The only possible observation is the main trend to negative values for all four fuels.

Figure 5 shows the trends of the three indices (variables) inserted in the discriminant model for the speed of 80,000 rpm.

**Figure 5.** Indices for 80,000 rpm: (**a**) non-normalized Shannon entropy, (**b**) asymmetry, and (**c**) root mean square value.

In Figure 5a, in which the trend of non-normalized Shannon entropy is shown, a clear distinction between pure and additive fuels can be seen, as in the case for 75,000 rpm. Pure diesel and kerosene present the lowest values for this index; diesel with 1% rapeseed oil gives higher values than the latter and diesel with 3% rapeseed oil gives the highest values among the four curves. This indicates how the vibrational dynamics of additive fuels are more irregular and unsteady than those of pure fuels and how the most irregular dynamics are in the case of diesel with 3% rapeseed oil, as expected, the latter being characterized by the higher percentage of rapeseed oil.

In Figure 5b, the asymmetry trend is reported, which tends to form two distinct clusters: the diesel group, including both pure and additive fuels, and the kerosene. The kerosene assumes negative values, which is not the case for diesel, where values are generally higher. This makes it possible to say that from a dynamic point of view, the use of diesel, pure or with rapeseed oil admixed, because it is characterized by a positive asymmetry, results in stronger vibrations than the use of kerosene, which, on the other hand, results in a much lower vibration. This was desirable since kerosene is the design fuel of the turbine.

Finally, Figure 5c shows the trend of the last of the three indices for the 80,000 rpm case, the RMS. The trend of the four curves allows a good distinction between pure and admixed fuels (note the separation of the red and the green curve from the other two). Furthermore, this graph emphasizes how the energy involved is greater in the case of pure fuels and, therefore, the explosive power decreases as the percentage of rapeseed oil admixed to diesel increases.

Figure 6 shows the trends of the three indices (variables) entered into the discriminant model for the speed of 85,000 rpm.

**Figure 6.** Indices for 85,000 rpm: (**a**) root mean square value, (**b**) quadratic oscillation index, and (**c**) asymmetry.

In Figure 6a, the RMS trend is shown. In this graph, the kerosene values are generally higher than those of the other fuels: as previously stated, this phenomenon was desired because kerosene is the turbine's design fuel, and the turbine works at a speed of only 5000 rpm lower than the maximum speed (90,000 rpm). Kerosene is the fuel with the best possible working conditions in terms of energy compared to all the others; on the other hand, the admixed diesel with the highest percentage of rapeseed oil has the worst working conditions in terms of energy.

In Figure 6b, a good distinction between pure diesel and admixed diesel can be seen, emphasizing how the turbine rotates smoothly and periodically using a pure fuel, as expected, compared to additive ones.

In the last diagram, i.e., the one in Figure 6c, the curves relating to admixed diesels are much more chaotic and assume a much wider range of values; those relating to pure fuels are much more stable and assume a much narrower range of values.

#### *4.2. Cluster Analysis*

As follows, the diagrams of the discriminant functions (Table 3) have been displayed for each speed. All functions are normalized in the range [0, 1].

Figure 7 shows the clusters of the discriminant function for the speed of 75,000 rpm.

**Figure 7.** Clusters for 75,000 rpm: (**a**) functions 1 and 2, (**b**) functions 2 and 3, and (**c**) functions 1 and 3.

Functions 1 and 2 (Figure 7a) explain most of the total variance, showing a clear distinction between pure diesel and kerosene, which occupy the right-hand side of the diagram, while on the other side there are the two additive diesel fuels, thus creating two families: pure and admixed fuels. It is possible to observe different behavior depending on whether a pure or non-pure fuel is used. Function 1 is associated with the quadratic oscillation index: the higher this index takes on values, the more the system is characterized by a better use of dynamic power. Indeed, it can be seen from the diagram that dynamic power is better used when using pure rather than admixed fuels: kerosene and diesel assume a

very high value of function 1 and, therefore, of the quadratic oscillation index compared to the admixed fuels. This is in full agreement with what was previously stated for the index diagram in Figure 4a, namely that using biofuels reduces emissions of pollutants at the expense of system efficiency. Function 2, which is obtained by adding to the term associated with the quadratic oscillation index a term relating to non-normalized Shannon entropy, the values are spread with respect to the centroids (yellow in the diagrams): as previously stated for this index, it highlights the presence of vibrational cycles that are not stationary and not periodic. It follows that the discriminant analysis at a speed of 75,000 rpm does not allow any clear distinction, since pure diesel describes power conditions like those of kerosene (function 1), while the admixed fuels assume high entropy values, which could compromise the integrity of the system.

Functions 2 and 3 (Figure 7b) do not allow a good distinction between the different types of fuels. Indeed, they explain only 13.1% of the total variance (Table 4) so that the different fuels make a single cluster, preventing any differentiation between the fuels.

Functions 1 and 3 (Figure 7c) allow a distinction between admixed fuels like in Figure 7a but less sharply: this is due to the lower discriminating power of function 3 compared to function 2 (9.3% of function 2 against 3.8% of function 3). Indeed, the overlap between the kerosene and pure diesel clusters is greater than the case in Figure 7a. Once again, the non-admixed fuels occupy the right-hand of the diagram, i.e., the dynamic power of the turbine is fully exploited (function 1 being composed of only the quadratic oscillation index), while the admixed fuels occupy the left-hand of the diagram, emphasizing an opposite behavior to the other two fuels. Finally, the cluster is much more dispersed with respect to the *y*-axis, particularly for the case of admixed diesel with 1% rapeseed oil.

Figure 8 shows the clusters of the discriminant function for the speed of 80,000 rpm.

**Figure 8.** Clusters for 80,000 rpm: (**a**) functions 1 and 2, (**b**) functions 2 and 3, and (**c**) functions 1 and 3.

The highest discriminating power is associated with functions 1 and 2 (Figure 8a), where it is possible to notice a pattern that clearly distinguishes the four fuels: particularly, the kerosene has values of function 1 different from the remaining fuels. Observing the non-normalized Shannon entropy (function 1), the points are really close to their own centroids: the system vibrates more contained than the 75,000 rpm case. In function 2, the entropy term is more prevalent than the asymmetry: from the diagram in Figure 5b, the kerosene values are clearly lower than the other three fuels, an event that is not repeated in the plot of the clusters (Figure 8a). Finally, it shows the tendency of admixed diesel with 1% rapeseed oil to pure diesel to be greater than diesel with 3% rapeseed oil: this is due to the different percentages of rapeseed oil added to diesel. In general, it can be seen in the diagram that the four clusters are distinct, separated from each other and closer to their own centroids: this allows us to say that it is the best classification obtained.

Since functions 2 and 3 (Figure 8b) have the least discriminating power, they do not allow a good distinction between the different types of fuels as in the case of Figure 8a. Although pure diesel and admixed diesel with 1% rapeseed oil are well-distinguished

and separated from the others, kerosene and admixed diesel with 3% rapeseed oil tend to overlap and form a single cluster.

Functions 1 and 3 (Figure 8c) allow a good distinction between the various fuels, although the cumulative variance percentage of 77.9% is lower than in functions 1 and 2 (96.7%). Even though they are slightly less sharply defined than in Figure 8a (this is due to the lower discriminating power of function 3 compared to function 2), the four clusters are still separated so that the diesel and the kerosene family stand out. Therefore, the considerations made previously on function 1 for Figure 8a are the same for the same function but in Figure 8c, i.e., the clusters are concentrated to their centroids and the system vibrates smoother than in the case of 75,000 rpm. Function 3 does not allow a good distinction between the fuels: the four clusters have an elongated shape along *y*-axis and a same range of values, except for admixed diesel with 1% rapeseed oil.

Figure 9 shows the clusters of the discriminant function for the speed of 85,000 rpm.

**Figure 9.** Clusters for 85,000 rpm: (**a**) functions 1 and 2, (**b**) functions 2 and 3, and (**c**) functions 1 and 3.

In the diagram of Figure 9a, knowing that function 1 is associated with the RMS, function 1 indicates how much energy is involved in terms of signal power where function 1, relatively to the kerosene, assumes greater values than the other fuels, as expected since kerosene is the design fuel, and the turbine is at a speed of 85,000 rpm (closer to the maximum speed). Function 2 allows a clear distinction between pure and admixed diesel: the first is clearly distinguished from the latter with a well-defined cluster, while the admixed diesel with 1% and 3% rapeseed oil overlap and there is not a big difference between the two. Kerosene has the highest power because it has the highest RMS. Since kerosene is the most expensive fuel among the proposed ones, a good alternative is using the pure diesel, which is far cheaper than the kerosene but less powerful (lower RMS). The pure and admixed diesel fuels have about the same function 1 range values: it means that it is possible to use admixed diesel keeping a similar RMS. Increasing the amount of rapeseed oil, the power output of the turbine and the pollutants decrease (this is also confirmed by the diagram in Figure 6a).

Functions 2 and 3 (Figure 9b) do not allow a good distinction between the various clusters: the four fuels tend to overlap and mix. As said for the previous speed rates, this event is due to the low percentage of cumulative variance provided by functions 2 and 3 (26.8%).

As in the case with functions 1 and 2, in the case of functions 1 and 3 (Figure 9c) it is possible to distinguish, on the left hand, the diesel family and, on the right hand, the kerosene, regarding function 1; furthermore, there is no distinction between pure diesel and admixed fuels. As regards function 3, it can be observed that the assumed values by pure fuels is much more concentrated around their centroid, unlike admixed fuels, which present much wider and dispersed values.

#### *4.3. Classification Summary*

Finally, the confusion matrices and the Cohen's kappa for each of the analyzed speeds are reported to highlight the goodness of the classification performed. In all confusion matrices, there are target classes along the first row and output classes along the first columns.

In the 75,000 rpm case (Table 5), it can be seen that 90.4% of the samples are correctly classified, which is an optimal result. In the kerosene case, 100% of the samples are correctly classified, while the lowest percentage is in the case of admixed diesel with 3% rapeseed oil where 19.2% of the samples are included in the case of admixed diesel with 1% rapeseed oil. Finally in the pure diesel case, 84% of the samples are correctly classified, while the remaining are confused with kerosene. It is interesting to note that it is never confused with admixed diesel, except for one sample confused with admixed diesel with 3% rapeseed oil: this underlines the clear separation between the two families of fuels shown in the clusters of Figure 7a.

In the case at 80,000 rpm (Table 6), 100% of the samples is correctly classified: this is the best possible case of classification since no value of the chosen variables will be associated with a different predicted class. This result was desirable from the simple observation of the clusters in Figure 8a, which are all distinct and separate from the others, and this confirmed that this is the best possible classification.

**Table 6.** Confusion matrix for 80,000 rpm.

In the case at 85,000 rpm (Table 7), 93.4% of samples are correctly classified. Both kerosene and pure diesel fuel contribute to this high percentage as both report a 100% correct classification. The lowest percentage of correct classification is for admixed diesel with 1% rapeseed oil: 80% is for the 1% rapeseed oil while the remaining 20% is assigned to 3% rapeseed oil. Moreover, for the admixed diesel with 3% rapeseed oil, 93.6% of samples are correctly classified, while the remaining are confused with admixed diesel with 1% rapeseed oil. This was also evident from the diagram in Figure 9a in which the clusters of the two admixed diesels tend to overlap and form a single cluster. It is possible to say that the case analyzed at 85,000 rpm constitutes a good distinction between pure fuels (the latter being correctly classified as 100%) and admixed fuels.

**Table 7.** Confusion matrix for 85,000 rpm.

To confirm the optimal results obtained with the proposed classification method, Table 8 shows the Cohen's kappa values for each turbine speed. The values of this index are always greater than 0.8: this implies an optimal classification, especially for the case at 80,000 rpm.

#### **5. Conclusions**

From the simple observation of the selected indices of the discriminant model, the following results could be deduced:

For all three analyzed speeds, the diagrams of discriminant functions 2 and 3 do not allow a good differentiation of the fuels as expected given the low percentage of cumulative variance (Table 4).

As regards the diagrams of functions 1 and 3, they certainly allow better differentiation than the case of functions 2 and 3 because of the higher percentage of cumulative variance, but they still do not allow optimal differentiation for each speed.

The best differentiation is provided by the diagrams obtained from functions 1 and 2, which provide the highest percentage of cumulative variance. Indeed:

• As for the previous speed, the case at 85,000 rpm led to the distinction between the diesel and the kerosene family but not an optimal distinction between the two admixed fuels because of overlapping clusters.

The results show that the best possible classification occurs in the case at 80,000 rpm where the clusters are well-defined and spaced out. This classification model can be used for quality check of the purchased fuel, especially diesel. It could happen that the purchased fuel has been admixed with crude oils, which can negatively affect its quality because of different molecular structures. Another application can be the quantity check of rapeseed oil admixed with diesel. Comparing the vibrational behavior of the turbine powered by the supplied fuel with the vibrational response of the same with the desired fuel, it is possible to say whether the fuel purchased from the supplier meets the required standards or not.

Finally, it is possible to make some final observations on the individual velocities with respect to the diagrams obtained by plotting functions 1 and 2:

From the confusion matrices and the calculation of Cohen's kappa, it was possible to establish the worth of the classification method studied and to see that, for all three velocities studied, optimal results were obtained.

The operating conditions of a gas micro-turbine powered with different types of fuel was discussed in this work. As discussed in the current literature, vibrational analysis has already been addressed to study the behavior of gas turbines through different approaches. The innovation of this work is the application of discriminant analysis to accelerometric signals, highlighting how it can provide indications regarding the differences in the system vibrations when powered with pure fuels or biofuels.

**Author Contributions:** Conceptualization, C.C., M.S. and A.N.; methodology, V.N., S.S. and G.Q.; software, C.C., M.S. and A.N.; validation, V.N., S.S. and G.Q.; investigation, C.C., M.S. and A.N.; data curation, C.C., M.S. and A.N.; writing—original draft preparation, C.C., M.S. and A.N.; writing review and editing, V.N., S.S. and G.Q.; visualization, C.C., M.S. and A.N.; supervision, V.N., S.S., G.Q. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper received no external funding.

**Data Availability Statement:** Not available.

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

#### **References**

## *Article* **Tools and Methods for Human Robot Collaboration: Case Studies at i-LABS**

**Massimo Callegari \*, Luca Carbonari, Daniele Costa, Giacomo Palmieri, Matteo-Claudio Palpacelli, Alessandra Papetti and Cecilia Scoccia**

> Department of Industrial Engineering and Mathematical Sciences, Universita Politecnica delle Marche, ` 60024 Ancona, Italy

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** The collaboration among humans and machines is one of the most relevant topics in the Industry 4.0 paradigm. Collaborative robotics owes part of the enormous impact it has had in small and medium size enterprises to its innate vocation for close cooperation between human operators and robots. The i-Labs laboratory, which is introduced in this paper, developed some case studies in this sense involving different technologies at different abstraction levels to analyse the feasibility of human-robot interaction in common, yet challenging, application scenarios. The ergonomics of the processes, safety of operators, as well as effectiveness of the cooperation are some of the aspects under investigation with the main objective of drawing to these issues the attention from industries who could benefit from them.

**Keywords:** human-robot collaboration; human-in-the-loop; virtual reality; obstacle avoidance

#### **Citation:** Callegari, M.; Carbonari, L.; Costa, D.; Palmieri, G.; Palpacelli, M.-C.; Papetti, A.; Scoccia, C. Tools and Methods for Human Robot Collaboration: Case Studies at i-LABS. *Machines* **2022**, *10*, 997. https://doi.org/10.3390/ machines10110997

Academic Editors: Huosheng Hu, Marco Ceccarelli, Giuseppe Carbone and Alessandro Gasparetto

Received: 22 July 2022 Accepted: 7 September 2022 Published: 30 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Until the introduction of Industry 4.0 principles, robotics in factories was mostly about machines replacing laborers who were tasked with non-ergonomic duties. The exploitation of robots was almost limited to manipulation of heavy loads or in uncomfortable positions, execution of dangerous tasks due to toxic payloads or environment, and execution of monotonous repetitive operations. Nowadays, thanks to the wide spread of *collaborative robots* (or *cobots*), the trend is shifting, especially in small and medium enterprises (SME's) [1]. In fact, small batch production and high level of product customization make these industrial entities still based on the versatility of human labor. Cobots, for their part, had a chance to easily insert themselves in this productive paradigm for they have been specifically developed for coexistence with people.

In this scenario, the collaboration of the researchers of Universita` Politecnica delle Marche together with five companies operating in the technology sector resulted in the creation of a laboratory which aims to disseminate the principles of Industry 4.0 in the local industrial fabric. Among them, collaborative robotics plays a crucial role, alongside with the development of tools and strategies for *Human Robot Collaboration* [2,3]. The aim is that of achieving the seamless team dynamics of an all-human team, and to this goal the last decade of research focused on several aspects [4,5] going from the machines themselves (exteroception [6,7], collision avoidance [8,9], intrinsic safety of the mechanics [10]), to their relation capabilities with people (interaction modalities [11,12], control oriented perception such as gesture recognition [13] and gaze tracking [14,15]).

Many relevant studies of the recent past can be mentioned trying to assess some of the HRC related aspects. Rusch et al. [16] quantified the beneficial impact, under both ergonomic and economic point of view, of preliminary simulated design of HRC scenarios. Papetti et al. [17] proposed a quantitative approach to evaluate a HRC simulated application. Similar topics were analysed in [18] where authors considered such opposing effects

on collaborative production lines. However, besides the economic impact, ergonomics analyses surely play a fundamental role for injury hazard management [19,20]. This becomes even more important considering the conclusions drawn by Dafflon et al. [21]: in light of their claims, completely unmanned factories are not reasonably possible mainly due to the un-feasibility of such complex control systems. Therefore, processes involving humans (*Human In The Loop*, HITL) will be implemented more and more frequently, justifying the effort of developing dedicated simulation environments and evaluation metrics for anthropocentric approaches [22–25].

However, the perception of operator commands in HITL tasks remains a widely investigated topic [26–28], since it still represents a critical point of the HRC process. Robustness, effectiveness, and above all safety are the keywords to be kept in mind. Thus, the interest aroused by strategies not involving contact among machines and workers is well understendable. Dinges et al. [29] considered the use of facial expressions to assess aggravated HRC scenarios; Yu et al. [30] used a multi-sensor approach to sense both postures and gestures to be interpreted as commands. In all of these examples, the key role played by artificial vision is evident, as also featured by authors of [31].

In the following, some of these topics are exemplified by case studies approached by the researchers at i-Labs. The aim is that of providing an impression of the ongoing research developed to pursue the principles of Industry 4.0 for what concerns the cooperation among humans and machines. In particular, a simulated environment is used to evaluate the ergonomics of an industrial work-cell in the first example. The second case considers the use of virtual reality for testing of obstacle avoidance control strategies applied to a redundant collaborative industrial manipulator. The third case study shows how a simple yet effective gesture recognition strategy can be developed to control tasks executed by a dual arm cobot. For the last case study, an obstacle avoidance control algorithm has been implemented on a collaborative robot to validate the effectiveness of the proposed law.

#### **2. Simulation Methods in HRC Design**

The method proposed in this section aims to support the human-oriented HRC design by merging the enabling technologies of industry 4.0 with the concepts of human ergonomics, safety, and performance. As shown in Figure 1, it starts with the human work analysis according to objective ergonomic assessments. The XSens™MVN inertial motion capture system is used to measure joint angles and detect awkward postures. The analysis results drive the collaborative robotic cell design; for example, allocating the nonergonomic tasks to the robot. The preliminary concept is then virtually simulated to identify critical issues from different perspectives (technical constraints, ergonomics, safety, etc.) and improve the interaction modalities. In these two phases, the NX (https://www.plm. automation.siemens.com/global/uk/products/nx/ (accessed on 1 January 2022)) and Tecnomatix (https://www.plm.automation.siemens.com/global/en/products/tecnomatix/ (accessed on 1 January 2022)) Process Simulate by Siemens are respectively used. The simulation contributes to the design optimization so that all the requirements are satisfied. The realization of the physical prototype allows the HRC experimentation, optimization, and validation before being implemented in the real production line. The experimentation phase includes a new ergonomic assessment to quantitatively estimate the potential benefits for the operator.

**Figure 1.** Simulation based method for human-oriented HRC design.

The industrial case study refers to the drawers' assembly line of LUBE Industries, the major kitchen manufacturer in Italy. Currently, a traditional robot and a CNC machine make the first part of the production line automated. Then, three manual stations complete the drawer assembly. Figure 2a shows one of the involved operators, equipped with 18 Xsens MTw (Wireless Motion Tracker), performing the tasks. The analysis involved all the labourers qualified for that specific task for a total of 6 operators (3 males, 3 females) all having average anthropometric characteristics. All the participants were informed about the goal of the study and the procedure. Also they were asked to read and sign the consent form.

The human work analysis highlighted a medium-high ergonomic risk (asymmetrical posture and stereotypy) for the operator dedicated to screwing. Before the objective evaluation, this workstation was mistakenly considered the most ergonomic of the three. Accordingly, the HRC design was based on the need to assign the screwing task to the robot to preserve the operator's health. Figure 2b shows the preliminary concept of the workstation that includes the following elements: conveyor belt; Universal Robot UR10e; collaborative screwdriver with a flexible extension; an aluminum structure where the robot is fixed by four M9 holes; L-shaped squaring system; clamping system; dispenser for feeding the screwdriver, easily and safely accessible by the operator to reload it; two plexiglass fences to reduce collisions risk.

The simulation highlighted several critical issues to be solved. To overcome some safety and technical problems, a new type of screwdriver was designed and implemented, which provides a maximum torque of 12 Nm and weighs 2.3 kg. To improve the line balancing the first two manual stations were aggregated in the cooperative robotic cell. However, the ergonomics simulation (Figure 2c) highlighted potential ergonomics risks for the operator during the drawer rotation.

**Figure 2.** (**a**) human work analysis; (**b**) preliminary concept of cooperative robotic cell; (**c**) HRC simulation; (**d**) physical prototype at i-Labs.

Then, an idle rotating roller was inserted so that the operator can manually turn the conveyor belt comfortably and safely. The simulation also showed that the new cooperative cell increased the current takt time. Consequently, an autoloader (FM-503H) was introduced to avoid unnecessary cobot movements. It consists of a tilting blade (suction) that loads the screws from a collection basket. Finally, the human-robot interaction has been improved by introducing the Smart Robots (http://smartrobots.it/product/ (accessed on 1 January 2022) ) vision system, which monitors the area of cooperation and enables the robot program according to human gestural commands. It allows the management of a high number of product variants.

In the new HRC workstation (Figure 2d), installed at i-Labs, negligible ergonomics risks, according to RULA (Rapid Upper Limb Assessment) and OCRA methods, were observed. As summarized in Table 1, also the Pilz Hazard Rating (PHR), which is calculated by (1), was reduced to acceptable values.

**Table 1.** Performance comparison between the "as-is" workstation, the preliminary concept of cooperative robotic cell and the final prototype.

$$PRR = DPH \times PO \times PA \times FE \tag{1}$$

where:

As shown, the new production line respects the cycle time (57 s) imposed by the CNC machine and the introduction of a collaborative robot led to better line balancing (+10%). The lab test was aimed at recreating the actual working conditions, therefore both a male and a female volunteer (actually two researchers of the i-Labs group) have been asked to take part in the investigation. Also, in this case, their anthropometric characteristics can be considered average.

#### **3. HITL: Gesture Recognition Examples**

In the search for effective ways to increase collaboration between workers and robots, new types of communication can be considered. Voice commands, for example, are typical in many systems today, such as vehicles, home automation systems, assistive robotics, etc., but they are not suitable for the industrial environment, where noise and the coexistence of multiple workers in a shared space make their implementation impossible. On the other hand, image sensors, such as standard or RGBD cameras, can be exploited to capture workers' gestures, which can then be interpreted by artificial intelligence algorithms to generate a command to be sent to the robot.

A distinction should be made between applications that require the use of wearable devices and those that rely on direct hand and body recognition. Regarding the former type, the most common technology for hand gesture recognition are glove devices. As part of their work at i-Labs, the authors developed the Robotely software suite for the specific purpose of overcoming the need for wearing gloves by adopting a hand recognition system using a machine learning approach [32]. As shown in Figure 3, the software exploits MediaPipe libraries in order to recognize 20 landmarks of the hand and identify up to three different hand configurations that can be converted to commands which can be sent to whatever robot or system. A specific test case was carried out at i-Labs: a PC with a standard webcam was used to run Robotely and connected to a Universal Robots UR10e cobot in order to sequence some operations between the operator and the robot.

**Figure 3.** Example of gesture recognition using a standard webcam and Robotely software [32].

A second example is shown in Figure 3 where the commercial Smart Robots system is integrated with the ABB YuMi cobot. Smart Robots is a programmable system based on RGBD sensors directly integrated into the robot control software. Again, three different hand gestures can be recognized and converted into commands, such as stop in Figure 4a, calling a program variation such as a quality inspection (Figure 4b), and opening the gripper for part rejection (Figure 4c). In addition, it is possible to define areas of the workspace corresponding to the activation of different program sequences: when the operator's hand enters a predefined area of the scene (once, or typically twice to avoid inadvertent activation) a program state is changed. In Figure 4d), for example, the operator covers a letter drawn on the table with his hand to order the robot to compose that letter with LEGO bricks.

**Figure 4.** (**a**–**d**): example of gesture recognition phases using Smart Robots on YuMi Robot.

#### **4. Safe Autonomous Motions: Collision Avoidance Strategies**

This section recalls the obstacle avoidance algorithm introduced by authors in [8,9] and presents some experimental results. The research in this field is presented hereby as an advancement in the direction of safe human-robot coexistence. Also, the resulting control algorithm is at the base of the test case presented in the subsequent section, which makes use of virtual reality tools to test and calibrate the empirical parameters introduced by the algorithm.

For a start, it is possible to consider the velocity kinematics of a generic manipulator owning a number of actuators greater or equal to 6 (depending on its degree of redundancy). In matrix form, the velocity kinematics of such a system can be written as **x**˙ = **Jq**˙ , where **J** is the 6 × *n* arm Jacobian (with *n* number of the joints composing the kinematic chain of the manipulator, *n* ≥ 6) and **q**˙ is the vector of joint rates. In such terms, the vector of variables is **q** = *q*1 ... *qn T* where *q*1 ... *qn* are the *n* joint variables of the arm kinematic chain. With such notation, the Jacobian **J** depends on the structure of the manipulator. The inverse of the Jacobian **J** of the redundant system can be obtained as a damped inverse:

$$\mathbf{J}^\* = \mathbf{J}^T (\mathbf{J} \mathbf{J}^T + \lambda^2 \mathbf{I})^{-1} \tag{2}$$

where *λ* is the damping factor, modulated as a function of the smaller singular value of the Jacobian matrix (the interested reader is addressed to [8,9] for further details). The inverse **J**∗ can be used to compute the joint velocities needed to perform a given trajectory with a Closed-Loop Inverse Kinematic (CLIK) approach:

$$\dot{\mathbf{q}} = \mathbf{J}^\*(\dot{\mathbf{x}} + \mathbf{K}\mathbf{e}) \tag{3}$$

where **x**˙ is the vector of planned velocities, **K** is a gain matrix (usually diagonal) to be tuned on the application, and **e** is a vector of orientation and position errors (**e***r* and **e***p*), defined as:

$$\mathbf{e} = \begin{bmatrix} \mathbf{e}\_r \\ \mathbf{e}\_p \end{bmatrix} = \begin{bmatrix} \frac{1}{2}(\mathbf{i} \times \mathbf{i}\_d + \mathbf{j} \times \mathbf{j}\_d + \mathbf{k} \times \mathbf{k}\_d) \\ \mathbf{P} - \mathbf{P}\_d \end{bmatrix} \tag{4}$$

In (4) the subscript *d* stands for desired planned variable, **P** is the position of the end-effector, while **i**, **j** and **k** are the unit vectors of the end-effector reference frame.

The collision avoidance strategy is then implemented as a further velocity component (to be added to the trajectory joint velocities) capable of distancing the end-effector and the other parts of the robot from a given obstacle. Such contribution is a function of the distance among the obstacle and every body of the robotic system. To this purpose, the bodies have been represented by means of two control points, **A** and **B** referring to Figure 5. Called **C** the centre of a generic obstacle, the distance among it and the segment *AB* can be differently computed in three different scenarios:

$$d = \frac{| (\mathbf{C} - \mathbf{A}) \times (\mathbf{B} - \mathbf{A}) |}{|\overrightarrow{AB}|} \tag{5}$$

• (**C** − **A**) *T*(**B** − **A**) ≥ |*AB*| 2, Figure 5c: point **C** is closer to any other point, therefore *d* = |*BC*|.

**Figure 5.** (**a**–**c**): obstacle to robot body distance in the three considered cases.

It is worth remarking that for this paper only spherical obstacles were considered, although similar approaches can be developed for objects of any shape starting from the distance primitives here defined.

Now the set of repulsive velocities for the ith body of the system can be introduced as:

$$\dot{\mathbf{q}}\_{r,i} = \psi v\_r \left( k \mathbf{J}\_{\mathbf{A},i}^\* \frac{(\mathbf{A} - \mathbf{C})}{|\overline{A}\overline{\mathbf{C}}|} + (1 - k) \mathbf{J}\_{\mathbf{B},i}^\* \frac{(\mathbf{B} - \mathbf{C})}{|\overline{B}\overline{\mathbf{C}}|} \right) \tag{6}$$

where:

At this point, the CLIK control law (3) can be completed as:

$$\dot{\mathbf{q}} = \mathbf{J}^\*(\dot{\mathbf{x}} + \mathbf{K}\mathbf{e}) + \sum\_{i=1}^m \dot{\mathbf{q}}\_{r,i} \tag{7}$$

being *m* the number of segments used to describe the manipulator.

At last, some experimental results are shown as demonstration of the obstacle avoiding strategy consistency. Going a little into details, the CLIK control law of the KUKA KMR iiwa (a 7 axes redundant collaborative manipulator) was built using 6 different segments, while the two obstacles were modelled as spheres. As shown by the test rig image (Figure 6), the robot carries a collaborative gripper by Schunk which is also considered for the definition of the robot segments (cfr Figure 7). The results are presented in terms of joints positions and rates, and pose and velocity of the robot end-effector (EE). Each variable is referred, for comparison, to the values collected during the execution of the same trajectory without obstacles in the robot workspace. The robot was controlled via Matlab with a cycle frequency of ∼ 150 Hz. It is worth remarking that the visualization provided in Figure 7 is a plot built on experimental data. The trajectory under investigation is a simple linear motion with constant orientation of the EE. Two obstacles have been put in the robot workspace, one directly on the planned trajectory (*O*1 in Figure 7) and the other in the space occupied by the robot non-terminal bodies (*O*2 in Figure 7). The results show how the control law is able to follow the given trajectory even in presence of the two obstacles. As expected, obstacle *O*1 (met approximately in the time span 2–4.5 s, gray area in graphs of Figure 7) prevents the robot EE from maintaining the desired pose: the control law permits the dodging of *O*1 deviating the bare minimum from the planned motion. Regarding obstacle *O*2 (met in the span 5.5–7.5 s, purple areas), the control law exploits the robot redundancy to avoid the contact among the robot elbow and the obstacle while maintaining the EE on the right trajectory and orientation.

**Figure 6.** Experimental test rig: the KUKA iiwa LBR avoiding two different obstacles in its workspace.

**Figure 7.** Experimental results: the KUKA iiwa LBR avoiding two different obstacles in its workspace.

#### **5. Virtual Reality Based Design Methods**

In this section a virtual implementation of the obstacle avoidance control strategy is presented. Aside the experimentation previously shown, the algorithm has been tested in advance in a simulated interactive scenario to demonstrate the feasibility of the control law, and to tune the parameters used to optimize the robot response to a dynamic obstacle. Such an approach, schematically shown in Figure 8, represents an efficient paradigm for design and virtual testing of HITL applications. Starting in a simulated environment, the interaction workflow can be designed considering also the presence of labourer and their impact on the production. A further optimization step can be added to optimize the interaction among humans and robots exploiting virtual reality. The VR tools permit one to disregard the modeling of labourer behaviour since their presence is played by actual humans who bring in the design process not only the actual expertise of the production line final users, but also the randomness of human actions and movements with further advantages in terms of security features design.

For this purpose, the idea was to create a *virtual reality* (VR) application in which a real obstacle is inserted in the simulation loop, while the robot is still completely virtual. The advantages of this approach are many: firstly, the algorithms can be modified and fine-tuned ensuring safety for the operators, secondly these changes can be made easily and quickly, reducing the development time of the real application.

The control architecture used for the VR application is shown in Figure 9, while Figure 10 is the experimental rig. The core of the system is the VR engine, running on a standard PC, developed by the SteamVR development suite with Unreal Engine 4 used to model the kinematics of the KUKA robot. A classic HTC VR set is used to equip the operator: the HMD provides an immersive three-dimensional representation of the workspace, which in this case is simply the robot mounted on a stand; the controller, manually managed by the operator, is rigidly connected to a virtual obstacle, shaped like a sphere, which can be moved in the virtual space to interfere with the robot's motion. The control of the robot is executed by the same PC in a parallel thread developed in Matlab, with a frame rate set to a typical value for communication protocols with robots (e.g.,125 Hz or 250 Hz, depending on the manufacturer). The frame rate of the VR engine is set to 60 Hz, which is the standard display refresh rate in VR applications. The communication between the two executed threads is realized by TCP/IP protocol. A classic pick and place task is simulated, thus the robot controller sends joint rotations to the VR engine at each time step, similarly to what is done in real robotic systems. The trajectory of the robot is updated in real time by the collision avoidance algorithm if the obstacle enters in a safety region of the manipulator. The position of the obstacle is known once the VR engine reads the coordinates of the controller held by the operator and sends this information to the robot controller.

**Figure 9.** Architecture of the control loop for the VR application.

**Figure 10.** VR system components.

#### **6. Concluding Remarks**

The paper showed the results obtained by the i-Labs laboratory in the wide field of research on human-robot interaction. Many technologies have been investigated and experimented to horizontally approach the issue of close cooperation among labourers and machines, going from off-line design, to on-line safety oriented control strategies.

The first case study showed a possible approach to the design and optimization of an assembly work-cell. The optimization, oriented at enhancing the labourer ergonomics and safety, allowed the realization of the prototype work-cell for the cobot assisted assembly of furniture pieces. The outcome quality has been quantitatively estimated.

Similar aspects of Human In The Loop operations have also been investigated in the second case study, which was about the realization of a visual based interaction paradigm among humans and robot. In this case, a depth camera was used to recognize the labourer commands which allowed the execution of simple tasks by the robot. This kind of contactless interaction pushed a little bit further in direction of security, which have been the core focus of the third case study.

The third case introduced an innovative control algorithm for active obstacle avoidance, applicable to any cobot. The strategy, based on the method of repulsive fields of velocity, was tested on a redundant industrial cobot to evaluate the possibilities offered in terms of both safety and task execution ability. At last, a Virtual Reality approach for the detail design an exploitation of the obstacle avoidance strategy has been described in the fourth case study.

The technologies which have been separately presented in this manuscript, actually represent some of the enabling technologies in the field of Industry 4.0 and 5.0 which place the human operator at the very center of the production process. Such a knowledge base can be exploited for the enhancement of industrial processes, especially in SMEs, which represent the main matter of the local industrial tissue.

**Author Contributions:** Conceptualization, G.P. and A.P.; methodology, G.P. and A.P.; software, L.C.; validation, D.C. and C.S.; formal analysis, L.C.; investigation, M.-C.P.; resources, M.C.; data curation, C.S.; writing—original draft preparation, L.C.; writing—review and editing, M.C.; visualization, C.S.; supervision, G.P. and A.P.; project administration, M.-C.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partly funded by the project URRA, "Usability of robots and reconfigurability of processes: enabling technologies and use cases", on the topics of User-Centered Manufacturing and Industry 4.0, which is part of the project EU ERDF, POR MARCHE Region FESR 2014/2020–AXIS 1–Specific Objective 2–ACTION 2.1, "HD3Flab-Human Digital Flexible Factory of the Future Laboratory", coordinated by the Polytechnic University of Marche.

**Data Availability Statement:** No publicly archived dataset has been used or generated during the research.

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

#### **References**

## *Article* **Configuration-Dependent Substructuring as a Tool to Predict the Vibrational Response of Mechanisms**

**Jacopo Brunetti 1, Walter D'Ambrogio 1,\* and Annalisa Fregolent 2**

**Abstract:** Dynamic substructuring allows us to predict the dynamic behavior of mechanical systems built by linking together several subsystems, whose dynamic behavior is known. The classical formulation, originally conceived for invariant systems, was extended by the authors to include mechanical systems made by invariant subsystems that may be coupled in different configurations. A mechanism is a typical example of a mechanical system built by coupling together invariant subsystems; during its motion, it can take several configurations that significantly affect its vibrational behavior. Therefore, the configuration-dependent substructuring approach can provide meaningful insights into the dynamic behavior of the mechanism. In this paper, the proposed approach is exploited to evaluate the vibrational behavior of a three-point linkage, a widely used mechanism to connect agricultural tractors to operating machines, considering a significant range of operative configurations. The proposed substructuring approach is able to predict the frequency response functions, the natural frequencies and the mode shapes of the mechanism in a wide range of configurations.

**Keywords:** configuration-dependent substructuring; dynamics of linkages; three-points linkage; vibrations

#### **1. Introduction**

Dynamic substructuring allows us to predict the dynamic behavior of mechanical systems built by coupling together several subsystems, whose dynamic behavior is known. The classical formulation [1] deals with invariant systems. Applications on configurationdependent systems are presented in [2,3] without exploiting dynamic substructuring. Dynamic substructuring was extended in recent years to include mechanical systems made by invariant subsystems that may assume different configurations during motion. Specifically, the problem of frictional sliding contact and of rolling contact are respectively considered in [4,5], while friction-induced vibrations are tackled in [6,7]. More generally, configuration-dependent substructuring can deal with mechanical systems composed of subsystems in relative motion with respect to each other, such as mechanisms and linkages as outlined in [8]. In fact, a linkage is a typical example of a mechanical system built by coupling together invariant subsystems; during its motion it can take several configurations that significantly affect its vibrational behavior. A classical approach to deal with the dynamics of mechanisms with deformable links is provided by the framework of multibody systems dynamics [9,10]. The multibody approaches are typically able to perform the kinematic and dynamic analysis of each link in the time domain by considering at the same time rigid body motion and vibrations. This is particularly useful in simulations aimed to verify the behavior of the mechanism under specified loads. However, a frequency domain analysis of such transient responses might not provide significant results for other design purposes mainly because the input forces might not properly excite the mechanism.

**Citation:** Brunetti, J.; D'Ambrogio, W.; Fregolent, A. Configuration-Dependent Substructuring as a Tool to Predict the Vibrational Response of Mechanisms. *Machines* **2022**, *10*, 1146. https://doi.org/10.3390/ machines10121146

Academic Editors: Marco Ceccarelli, Giuseppe Carbone and Alessandro Gasparetto

Received: 2 August 2022 Accepted: 28 November 2022 Published: 1 Decmber 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

On the contrary, configuration-dependent substructuring can be directly formulated in the frequency domain. On one side, this allows the description of a complex link, that would be difficult to model properly, in terms of an experimentally determined Frequency Response Function. So far, multibody approaches do not allow the use of experimentally identified components. On the other side, the vibrational behavior of the mechanism can be analyzed at every position using the configuration-dependent frequency response function from which configuration-dependent natural frequencies and vibration modes can be identified. This quasi-static information set is very important for design purposes since it highlights most of the potential dynamic problems that can occur during the mechanism's operation.

Since the proposed approach is quasi-static and a linearization of the system is performed on each configuration, the effects of friction and of nonlinearities in general, are not accounted. Friction contact problems can be accounted in the subsructuring framework as shown in [6,7]. However, in this case, the problem must be solved in the time domain and no experimentally derived models can be used. For these reasons, the configurationdependent substructuring approach cannot be seen as an alternative to multibody analysis but rather as a technique providing different and complementary information.

The original contribution of this paper lies in envisaging the application of substructuring techniques, typically devised for structures, to mechanisms and linkages. Specifically, as a proof of concept, configuration-dependent substructuring is applied on a three-point linkage, a widely used mechanism to connect agricultural tractors to operating machines [11]. The vibrational behavior of this system is evaluated throughout a significant range of operative configurations. Frequency response functions, natural frequencies and mode shapes of the mechanism are predicted for each analyzed configuration.

#### **2. Substructure Coupling in the Frequency Domain**

A dynamic system made up of *n* coupled subsystems is considered. Each subsystem *r* can be described using the mass, stiffness and damping matrices **M**(*r*), **K**(*r*) and **C**(*r*), from which the dynamic stiffness matrix can be computed as **Z**(*r*)(*ω*) = **K**(*r*) − *ω*2**M**(*r*) + *iω***C**(*r*).

For a given linear and invariant subsystem *r*, the equation of motion can be expressed in the frequency domain as:

$$\mathbf{Z}^{(r)}(\omega)\boldsymbol{\mu}^{(r)}(\omega) = f^{(r)}(\omega) + \mathbf{g}^{(r)}(\omega) \tag{1}$$

where:

**Z**(*r*)**:** dynamic stiffness matrix of subsystem *r*;

*u*(*r*)**:** vector of displacements of subsystem *r*;

*f*(*r*)**:** vector of external forces acting on subsystem *r*;

*g*(*r*)**:** vector of connecting forces with other subsystems (internal constraints).

By writing the equation of motion of the *n* subsystems in a block diagonal format, it is obtained, after leaving out the frequency dependence:

$$\mathbf{Z}\mathbf{u} = f + \mathbf{g} \tag{2}$$

with

$$\mathbf{Z} = \begin{bmatrix} \mathbf{Z}^{(1)} & & & \\ & \ddots & & \\ & & \mathbf{Z}^{(n)} \end{bmatrix}, \quad \mathbf{u} = \begin{Bmatrix} \mathbf{u}^{(1)} \\ \vdots \\ \mathbf{u}^{(n)} \end{Bmatrix}, \quad \mathbf{f} = \begin{Bmatrix} \mathbf{f}^{(1)} \\ \vdots \\ \mathbf{f}^{(n)} \end{Bmatrix}, \quad \mathbf{g} = \begin{Bmatrix} \mathbf{g}^{(1)} \\ \vdots \\ \mathbf{g}^{(n)} \end{Bmatrix}$$

To couple together the *n* subsystems, compatibility and equilibrium conditions must be enforced. Compatibility at the interface DoFs means that any pair of corresponding DoFs *u*(*r*) *l* and *u*(*s*) *m* , i.e., DoF *l* on subsystem *r* and DoF *m* on subsystem *s* must share the same displacement, that is *u*(*r*) *l* − *u*(*s*) *m* = 0.

Generally, this condition can be written as:

$$
\mathbb{B}u = \mathbb{0} \tag{3}
$$

where each row of **B** refers to a pair of corresponding DoFs.

Equilibrium implies that internal constraint forces must be balanced. They arise when two subsystems are connected together at a pair of corresponding DoFs. Therefore, the sum of internal constraint forces must be zero for any pair of corresponding DoFs, i.e., *g* (*r*) *l* + *g* (*s*) *m* = 0.

Moreover, if DoF *k* on subsystem *q* is not a coupling DoF, *g* (*q*) *k* = 0. Generally, the equilibrium conditions can be written as:

$$\mathbf{L}^T \mathbf{g} = \mathbf{0} \tag{4}$$

where **L** is a localization matrix.

The system of Equations (2)–(4) provides the so-called three-field formulation, defining the coupling between any number of subsystems:

$$\begin{cases} \mathbf{Z}\boldsymbol{\mu} = \boldsymbol{f} + \mathbf{g} \\ \mathbf{B}\boldsymbol{\mu} = \mathbf{0} \\ \mathbf{L}^T \mathbf{g} = \mathbf{0} \end{cases} \tag{5}$$

In the dual assembly [1,12], all DoFs are retained, i.e., each coupling DoF among two substructures appears twice. The equilibrium condition *g* (*r*) *l* + *g* (*s*) *m* = 0 at a pair of coupling DoFs is ensured by selecting *g* (*r*) *l* = −*λ* and *g* (*s*) *m* = *λ*. Therefore, the connecting forces can be written in the form:

$$\mathbf{g} = -\mathbf{B}^T \boldsymbol{\lambda} \tag{6}$$

where the Lagrance multipliers *λ* represent the intensities of connecting forces.

The equilibrium condition (4) can be rewritten as:

$$\mathbf{L}^T \mathbf{g} = -\mathbf{L}^T \mathbf{B}^T \boldsymbol{\lambda} = \mathbf{0} \qquad \forall \boldsymbol{\lambda} \tag{7}$$

The Equation (7) is always satisfied, thus the three-field formulation (5) reduces to:

$$\begin{cases} \mathbf{Z}u + \mathbf{B}^T \lambda = f \\ \mathbf{B}u = \mathbf{0} \end{cases} \tag{8}$$

To eliminate *λ*, the following steps can be performed. The first of Equation (8) becomes:

$$
\mu = -\mathbf{Z}^{-1}\mathbf{B}^T\boldsymbol{\lambda} + \mathbf{Z}^{-1}\boldsymbol{f} \tag{9}
$$

Equation (9) can be substituted in the second of Equation (8) giving:

$$\mathbf{B}\mathbf{Z}^{-1}\mathbf{B}^{T}\boldsymbol{\lambda} = \mathbf{B}\mathbf{Z}^{-1}\boldsymbol{f} \quad \Rightarrow \quad \boldsymbol{\lambda} = \left(\mathbf{B}\mathbf{Z}^{-1}\mathbf{B}^{T}\right)^{-1}\mathbf{B}\mathbf{Z}^{-1}\boldsymbol{f} \tag{10}$$

By substituting *λ* in the first of Equation (8), one finally gets:

$$\mathbf{Z}u + \mathbf{B}^T \left(\mathbf{B}\mathbf{Z}^{-1}\mathbf{B}^T\right)^{-1} \mathbf{B}\mathbf{Z}^{-1}f = f$$

$$\Rightarrow \qquad u = \left(\mathbf{Z}^{-1} - \mathbf{Z}^{-1}\mathbf{B}^T \left(\mathbf{B}\mathbf{Z}^{-1}\mathbf{B}^T\right)^{-1}\mathbf{B}\mathbf{Z}^{-1}\right)f\tag{11}$$

Since **H**(*r*) = [**Z**(*r*) ] −1 is the Frequency Response Function (FRF) matrix of the *r*-th subsystem, it can be written:

> **Z**−1 = **H** = ⎡ ⎢ ⎢ ⎣ **H**(1) ... **H**(*n*) ⎤ ⎥ ⎥ ⎦ (12)

Therefore, Equation (11) becomes: