2.1.1. Model scale

Once the CAD model of the power plant to be inspected is obtained, it is appropriately scaled down to fit the software working environment. We define the model scale $S_L$ as the ratio of length of the facility $L_p$ to the length of the CAD model $L_m$ .

(1) \begin{equation}S_L=\frac{L_p}{L_m} \end{equation}

2.2.1. Camera field of view

For a chosen dimension d, and effective focal length f, the angular field of view $\alpha$ can becalculated as,

(2) \begin{equation}\alpha =2\arctan\left(\frac{d}{2f}\right) \end{equation}

The camera field of view can be measured horizontally, vertically, and diagonally (Fig. 3). For an Intel Real®Sense R200 camera, the horizontal $\alpha_h$ vertical $\alpha_v$ diagonal $\alpha_d$ angular field of view for the RGB camera, are 60°, 45° and 70°, respectively [Reference Keselman, Woodfill, Grunnet-Jepsen and Bhowmik24].

Figure 3. UAV inspection camera angular FOV (field of view).

2.2.2. The UAV system

The aerial platform used in this project consists of an off-the shelf quadrotor vehicle Intel® Aero Ready to Fly Drone (Fig. 4). However, this paper focuses on the offline trajectory generation and validation in a simulation environment to show the capabilities of the approach. All data and parameters are taken from the said UAV system for the corresponding calculations.

Figure 4. Drone platform and reference frames.

2.3.1. UAV kinematics

To analyze the motion of the drone, an inertial frame {A} (or world frame) and a body coordinate frame {B} are included. {A} is selected at a convenient location that acts as the global reference for the plant to be inspected. In the other hand, {B} has its origin coincident with the center of gravity of the quadrotor, the z axis downward, and the x axis in the direction of the flight, following the aerospace convention (Fig. 4). Besides, as depicted in Fig. 5, the position of the drone in inertial frame {A} is given by the vector $\boldsymbol{r}\in\mathbb{R}^{3}$ .

Figure 5. Kinematic Diagram: frames definition and inspection distances in a power plant boiler case study.

The inspection cameras on the UAV are installed in the chassis and are near the depth sensors; therefore, the camera location can be considered as the point of origin/zero point for the drone. Notably, the UAV propellers expand beyond this point by $a_d$ and from the extreme endpoint of the propeller to the closest point on the inspection surface can be designated as $d_s$ . As such, the total inspection offset from camera at any time is

(3) \begin{equation}\ d =d_{s}+a_{d} \end{equation}

For a safe flight path planning, a minimum drone offset distance $d_s$ is introduced. Also, due to the design of the drone, the rotor blades expand a certain distance from the camera that also needs to be considered along with the camera focal distance itself. Therefore, the offset distance should be $d \geq d_s+a_d$ and $d>f$ at any time (Fig. 5). As such, the position vector for the UAV camera is

(4) \begin{equation} {{\boldsymbol{r}} } = \begin{bmatrix} x+d & y & z \end{bmatrix}^{T}\end{equation}

For a quadrotor, the angular motion state about the three rotation axes is represented by six state variables: the angular velocity vector, $\boldsymbol{\omega\in\mathbb{R}^{3}}=[\omega_x, \omega_y, \omega_z]^T$ about the three axes of the body frame and the Euler angles $\boldsymbol{\beta\in\mathbb{R}^{3}}=[\phi,\theta,\psi]^T$ , representing the roll, pitch and yaw angles, respectively. The 3-2-1 (Z-Y-X) Euler angles are used to model the rotation of the quadrotor in the global frame. The rotation matrix for transforming the coordinates from {A} to {B} is given by the rotation matrices about each axis ( R ^′ $(\phi)$ , R ^′ $(\theta)$ , and R ^′ $(\psi)$ ). Therefore, the rotation matrix $\boldsymbol{R}\in\mathbb{R}^{3 \times 3}$ that expresses the orientation of the coordinate frame {B} with respect to global reference frame {A} is obtained.

(5) \begin{equation} \boldsymbol{R} = \boldsymbol{R}'{\rm (\phi)} \boldsymbol{R}'(\theta) \boldsymbol{R}'(\psi) = \begin{bmatrix} C_\theta C_\psi & C_\theta C_\psi & -S_\theta \\[4pt] S_\phi S_\theta C_\psi & C_\phi C_\psi + S_\phi S_\theta S_\psi & S_\phi C_\theta \\[4pt] C_\phi S_\theta C_\psi + S_\phi C_\psi & C_\phi S_\theta C_\psi - S_\phi C_\psi & C_\psi C_\theta \end{bmatrix}\end{equation}

where $ C_\phi $ and $ S_\phi $ are abbreviations of cos( $\phi$ ) and sin( $\phi$ ), respectively, similarly for $\theta$ and $\psi$ . The components of angular velocity of the quadrotor, $\omega_x$ , $\omega_y$ and $\omega_z$ are related to the derivatives of the roll, pitch and yaw angles as,

(6) \begin{equation} \begin{bmatrix} \omega_x \\[4pt] \omega_y \\[4pt] \omega_z \end{bmatrix} = \begin{bmatrix} 1 & 0 & -S_\theta \\[4pt] 0 & C_{\phi} & S_{\phi} C_\theta \\[4pt] 0 & - S_{\phi} & S_{\phi} \end{bmatrix} \begin{bmatrix} \dot{\phi} \\[4pt] \dot{\theta} \\[4pt] \dot{\psi} \end{bmatrix}\end{equation}

2.3.2. UAV dynamics

The four rotors on the quadrotor are driven by electric motors powered by electric speed controllers. The rotor speed is $\overline{\omega_{i,}}$ and the thrust is an upward vector,

(7) \begin{equation}T_i = b \overline{\omega_i}, \quad i = 1,2,3,4 \end{equation}

In the quadrotor’s negative z-direction, where b $>$ 0 is the lift constant, dependant on number of blades, air density, the cube of the rotor blade radius and the chord length of the blade. Assuming the quadrotor is a symmetric rigid body, its equations of motion can be written as [Reference Corke25],

(8) \begin{equation}\boldsymbol{f}= \boldsymbol{f}_g + T \end{equation}

(9) \begin{equation} \boldsymbol{\tau} = \boldsymbol{J}{\boldsymbol{\dot\omega}} + \boldsymbol{\omega} \times \boldsymbol{J}\boldsymbol{\omega}\end{equation}

where $\boldsymbol{f}\in\mathbb{R}^3$ and $\boldsymbol{\tau}\in\mathbb{R}^{3}$ are the force and torque applied to the vehicle, respectively; $\boldsymbol{J}\in\mathbb{R}^{3\times 3}$ is the inertia matrix of the quadrotor, and $\boldsymbol{\omega}\in\mathbb{R}^3$ is the angular velocity vector of the UAV, $f_g$ is the force due to the gravity, and T is the total thrust generated by the propellers. The resultant force can be then expressed as

(10) \begin{align} \boldsymbol{f} &= \begin{bmatrix} 0 \\[4pt] 0 \\[4pt] mg \end{bmatrix} + \boldsymbol{R} \times \begin{bmatrix} 0 \\[4pt] 0 \\[4pt] T_1 + T_2 + T_3 + T_4 \end{bmatrix}\end{align}

Pairwise differences in rotor thrusts cause the vehicle to rotate. The rolling torque about the vehicle’s x axis is generated by the moments

(11) \begin{equation} \tau_x = a_{CM} T_4 - a_{CM} T_2 \end{equation}

$a_{CM}$ is the distance from the quadrotor center of mass to the rotor axis (Fig. 5). Converting the equation in terms of rotor speeds using (7)

(12) \begin{equation}\tau_x = a_{CM} b (\overline{\omega_4}^2 - \overline{\omega_2}^2) \end{equation}

Similarly, for the torque about its y axis,

(13) \begin{equation}\tau_y = a_{CM} b (\overline{\omega_1}^2 - \overline{\omega_3}^2) \end{equation}

The torque applied to each propeller by the motor needs to overcome the aerodynamic drag,

(14) \begin{equation}{Q}_{i} = k \overline{\omega_i}^2 \end{equation}

where k depends on the same factors as b. From the torque, a reaction torque on the airframe is transmitted that is responsible for the “yaw” torque about the propeller shaft in the opposite direction to its rotation. The total torque about z axis is

(15) \begin{align} \begin{aligned}\tau_z &= Q_1 - Q_2 + Q_3 -Q_4 \\[4pt] &= k(\overline{\omega_1}^2 - \overline{\omega_2}^2 +\overline{\omega_3}^2 -\overline{\omega_4}^2)\end{aligned} \end{align}

where the signs are different due to the rotation directions (counter clockwise taken as positive for rotor 1 and 3) of the rotors. Therefore, yaw rotation is achieved by carefully controlling all of the 4 rotor speeds. The applied torque $\boldsymbol{\tau}$ is calculated from (12), (13), (15) and (7) obtaining

(16) \begin{align} \boldsymbol{\tau} &= \begin{bmatrix} \tau_x \\[4pt] \tau_y \\[4pt] \tau_z \end{bmatrix}=\begin{bmatrix} a_{CM}b(\overline{\omega_4}^2 - \overline{\omega_2}^2) \\[4pt] a_{CM}b(\overline{\omega_1}^2 - \overline{\omega_3}^2) \\[4pt] k(\overline{\omega_1}^2 - \overline{\omega_2}^2 +\overline{\omega_3}^2 -\overline{\omega_4}^2) \end{bmatrix} = \boldsymbol{J}{\boldsymbol{\dot \omega}} + \boldsymbol{\omega}\times \boldsymbol{J} \boldsymbol{\omega} \end{align}

By combining Eqs. (7) and (10), the following expression is obtained

(17) \begin{align} \begin{bmatrix} T \\[4pt] \boldsymbol{\tau} \end{bmatrix} &= \begin{bmatrix} -b & -b & -b & -b \\[4pt] 0 & -a_{CM}b & 0 & a_{CM}b\\[4pt] a_{CM}b & 0 & -a_{CM}b & 0\\[4pt] k & -k & k & -k \end{bmatrix} \begin{bmatrix} \overline{\omega_1}^2 \\[4pt] \overline{\omega_2}^2 \\[4pt] \overline{\omega_3}^2 \\[4pt] \overline{\omega_4}^2 \end{bmatrix} = \boldsymbol{A} \begin{bmatrix} \overline{\omega_1}^2 \\[4pt] \overline{\omega_2}^2 \\[4pt] \overline{\omega_3}^2 \\[4pt] \overline{\omega_4}^2 \end{bmatrix}\end{align}

The moments and the forces of the quadrotor are function of the rotor speeds. The matrix A is constant, and the equation can be rewritten using matrix inversion for obtaining the rotor speeds given the torques and forces [Reference Corke25],

(18) \begin{align} \begin{bmatrix} \overline{\omega_1}^2 \\[4pt] \overline{\omega_2}^2 \\[4pt] \overline{\omega_3}^2 \\[4pt] \overline{\omega_4}^2 \end{bmatrix} = \boldsymbol{A}^{-1} \begin{bmatrix} T \\[4pt] \tau_{x} \\[4pt] \tau_{y} \\[4pt] \tau_{z} \end{bmatrix}\end{align}

2.3.3. UAV control system

As depicted in Fig. 6, the UAV position and attitude are controlled by a series of controllers built-in in the system. A position controller receives the desired coordinates to track in the position vector set-point $\boldsymbol{r}_{sp}\in\mathbb{R}^{3}$ that is compared with the actual estimated drone position $\boldsymbol{\hat{r}}\in\mathbb{R}^{3}$ to obtain the position error $\boldsymbol{\tilde{r}}\in\mathbb{R}^{3}$ to compensate for. The desired coordinates are the final coordinates obtained from the “Post-processing/toolpath extraction phase.” The UAV’s X, Y and Z axis are set via the flight controller, inline with the AM/CAM software, and the origin of the inertial frame, {A}, (Fig. 5) is where the UAV is placed to start the aerial inspection flight. A proportional controller calculates then the desired velocity $\boldsymbol{v}_{sp}\in\mathbb{R}^{3}$ that is limited according to the hardware capabilities. Then, a PID controller in cascade will determine the required thrust T to achieve the desired altitude. While, through a conversion stage, the thrust and the desired yaw angle $\psi_{sp}$ determine the desired attitude $\boldsymbol{\beta}_{sp}\in\mathbb{R}^{3}$ that is compared with the estimated attitude $\boldsymbol{\hat{\beta}}\in\mathbb{R}^{3}$ to obtain the orientation error $\boldsymbol{\tilde{\beta}}\in\mathbb{R}^{3}$ that will be compensated by a proportional controller, leading to the attitude rate set point $\boldsymbol{\dot{\beta}}_{sp}\in\mathbb{R}^{3}$ that at the same time is related to the estimated attitude rate $\boldsymbol{\dot{\hat{\beta}}}\in\mathbb{R}^{3}$ to determine the attitude rate error $\boldsymbol{\dot{\tilde{\beta}}}\in\mathbb{R}^{3}$ that will be regulated by a PID controller that generates the control torques, that together with the thrust command the drone to the desired position and orientation. Finally, the required motor’s angular speed can be calculated by (18). References [finite] and [global] The estimated attitude and attitude rates are directly inferred from the IMU data, while the position estimation is more complex and cannot rely on a GPS signal. For solving this issue, we propose implementing a well-tested visual inertial odometry (VIO)-based sub-system [Reference Nikolic, Burri, Rehder, Leutenegger, Huerzeler and Siegwart7] using a stereo camera integrated in the payload. As for the simulation, the actual UAV’s position is corrupted by noise, which simulates the VIO system’s output working in a feature richenvironment.

Figure 6. UAV controllers for position and attitude control.

3.1.1. Model scaling

Boiler height commonly ranges between 12–30 m. For this study, a 17 m high boiler is designed. Due to the tool size limitation, overall power plant height is scaled down to 170 mm. Therefore, model scale $(S_L)$ is 100:1. The $S_L$ is good for all use cases, including actual case scenario of a real-life boiler height range (typically 12-30m). The AM tool ’‘Slicer’’ is also good for scaling and trajectory planning of Z axis from 120-300 mm range. Once trajectory planning is done in model phase, the coordinates or the UAV waypoints is converted to the real-life scaling using the scale of 100:1, such that a 100 mm model distance between coordinates converts to 10 m in actual use.

3.1.2. Camera field of view

For calculation, the stereo camera angle of view is taken. Average camera focal length is 30 mm. Here, a 500 mm inspection offset distance for UAV camera is chosen. From (2), with $\alpha_{h} = 60^{\circ}$ , $\alpha_{v} = 45^{\circ}$ , and $f = 5$ mm, the inspection grid horizontal length, h, and vertical length, v, are calculated:

(19) \begin{equation}h = 2f\tan\frac{\alpha_h}{2}\! \Rightarrow h = 5.8\end{equation}

(20) \begin{equation}v = 2f\tan\frac{\alpha_v}{2}\! \Rightarrow v = 4.1\end{equation}

Considering a 30% overlap of adjacent pictures for a sharp and high image quality inspection profile, the camera inspection field of view is taken as,

(21) \begin{equation}h_1 = 70\% of h = 4.06 \approx 4 \end{equation}

(22) \begin{equation}v_1 = 70\% of v = 2.87 \approx 2.5 \end{equation}

Therefore, the model scale yields to vertical length field of view ( $v_1$ ) of 2.5 mm; hence, this is the maximum step down/feed for CAM and layer to layer height in AM. From the scale ratio, the application field of view is calculated to be 400 mm $\times$ 250 mm approximately (Fig. 7) and the number of parallel toolpath (slices) is 68.

Figure 7. UAV camera FOV (field of view) grid.

3.1.3. Drone velocity range

Velocity of drones range from 10 mph to as high as 70 mph. For this project, an average speed needs to be calculated. Average drones used for photography usually clock between 30–50 mph, for which the drone camera must have a suitable fps (frames per second) value. For the maximum speed, $v_{max} =$ 50 mph $= 22.235 $ m/s, the frame per second (fps) for the drone camera is the ratio of total distance to the camera grid horizontal length $(h_1)$

(23) \begin{align} \begin{aligned} {\rm fps} & = \frac{22,352}{400} \\ &= 55.88 \end{aligned} \end{align}

Taking an average speed of 20 mph or 8.9408 m/s, that is, the drone would need to travel 8940.8 mm in 1 sec. Therefore,

(24) \begin{align}\begin{aligned} {\rm fps} &= \frac{8940.8}{400} \\ &= 22.352 \end{aligned} \end{align}

The proposed camera has 60 fps for a depth resolution of $480\times360$ , which is well within our calculated range. Also, as mentioned previously, a 30% overlap has been accounted for a higher image quality. As such, both photography and video for inspection can be used well with the inspection payload for an average drone speed of 20 mph.

3.2.1. Method 1: CAM toolpath generation

A. Toolpath Generation with 2D Contour Operation: In the “Manufacture” workspace of Fusion 360, the machining properties are setup, that is, setting up the appropriate operation type, defining work coordinate system, choosing apt stock offset, post-processing setup, etc.

At first, a relevant sized stock block is chosen. The stock offset is equal to the inspection offset, “d” (Eq. (3)). In order to optimize the toolpath and thereby obtaining a proper UAV flight trajectory, the toolpath needs to be uniformly distributed and well wrapped around the model (Fig. 8). As such, a 2D contour milling operation in the CAM environment of the Fusion 360 software is selected. All machining parameters are aptly chosen which gives a uniform toolpath around the powerplant (Table I). The endpoint of a 2.5 mm diameter ball end point cutting tool chosen for this operation, machines out the part from the stock and a “toolpath’’ is obtained.

Figure 8. Toolpath generated around the powerplant CAD model for 2D contour operation.

Table I. CAM Contour (2D) machining parameters.

Two very important parameters are maximum step down and stock offset. “Maximum step down” parameter represents the length that is traveled when the cutting tool goes to next layer, and this corresponds to the hatch distance value in AM (slice). “Stock offset” parameter corresponds to the vertical distance from the inspection surface which is correlated to the distance of the camera to capture inspection profile photos and videos. These parameters are recommended considering the minimum and smallest distance value that can result in dense inspection profile; however, different values can be set for these parameters as operation requires. For example, here a minimum safe distance is taken as500 mm, but it can be increased depending on the type of UAV and the capability of its camera for producing desired outcome.

Mentionable is that the ‘‘blue’’ path lines describe cutting the stock or machining while the ‘‘yellow’’ and ‘‘red’’ path lines describe tool retraction and ramping, respectively. Notably, due to the tool orientation, the overhang, gaps and inward slopes are avoided, as expected from a CAM operation. Although the parallel toolpath around the inward slope of the exhaust is acceptable due to the camera’s superior zoom capabilities, however, to get the best possible wrapped toolpath of the section, tool orientation is changed.

B. Toolpath Generation with Modified Tool Orientation: For the inner gap of the powerplant model, again 2D contour milling operation is chosen, other parameters are the same (Table II). The Z axis of the tool orientation along with Y axis are flipped 180° (Fig. 9(a)). Z axis directs the tool entry to the workpiece. For this case, it is the same as assuming that the workpiece has been flipped for operating from the bottom surface of the setup, whereas the X axis is setup according to the manufacture operation coordinate principle. The inner gap between the boiler and exhaust can also be addressed with “Slot” milling operation (Fig. 9(b)). This is a step toolpath that requires a different set of parameters (Table II).

Table II. CAM Slot machining parameters.

Figure 9. Toolpaths generated with modified tool orientation and operation for the component gap and inward slope of the powerplant CAD model.

Like the earlier operation, Z axis is pointing downwards (flipped). Although the toolpath is able the reach the shallow gap on the upper chamber of the exhaust, due to the aerodynamic restriction of a UAV, it is not considered and the operation ceases near the shallow zone of the “particulate recirculation” component of the power plant.

Overall for all operations, ‘‘smooth transition’’ and ‘‘plunge in’’ features were chosen to avoid sharp change in UAV flight heading. From the generated toolpath, additional noise is observed for the CNC tool safe retrievals/lead-ins/lead-outs. Mentionable is that the tool time information, for example, feed rate, revolution, surface speed, is not used in the final flightpath generation. Only the tool end point coordinates are extracted as point cloud; therefore, these are not essential to this project. Later, an UAV appropriate timestep is set up.

3.2.2. Method 2: AM extruder path

In order to obtainthe external coordinates with a pre-set offset, a CAM-based approach was used previously. The CAM base method has some challenges that are inherent to the CAM characteristics, for example:

1. 1. Machining of certain features is difficult, for example, inward slopes larger than 45, convex surfaces, etc.

2. 2. Difficult to reach shallow gaps between design parts due to tool size limitation. The tool tilt angle is also limited due to physical constraint and workpiece orientation etc.

To overcome these challenges, an “Additive Manufacture (AM)”-based method is introduced with significant advantages. Essentially, a path that is followed in a CAM milling operation and a 3D printer (in AM) is very similar, one being subtractive, tool traveling in -Z axis (top to bottom) and the other being additive, that is, tool travels in +Z axis (bottom to top) (Fig. 10). During the operation of both manufacturing methods, the tool moves one layer at a time.

Figure 10. 3D printing layers preview of the model in Cura 3.6.0.

A commercially available 3D printer software (slicer) “Ultimaker Cura 3.6.0” is used to generate a 3D shell of the power plant with a pre-set offset equal to 5 mm, similar to CAM method (Table III). The layer height corresponds to the vertical travel of the UAV during inspection. No structure support is chosen, and mostly default parameters have been opted for getting a general shape of the powerplant. This operation gives extruder/toolpath coordinates at the selected offset or the shell itself. This covers all surface areas, including the shallow gaps and slopes. The G-code generated a bottom to upward Z axis values, along with XY coordinates, layer by layer (Fig. 10). Considering all the advantages in the final trajectory generation, the AM method is finalized for use.

Table III. AM parameters.

3.3.1. Direct data extraction

Using a post-process add-on, the Cartesian coordinates of the tool endpoints are exported (Fig. 11). This is the first step of path planning. There are 109,575 number of points in the ‘‘exhaust’’ component of the powerplant setup alone. The point cloud is very dense; hence, a snapshot of the coordinates is shown (Fig. 11). In this step, all the points are taken and organized under one global coordinate. In the second and third operations, tool axis, that is, Z axis and Y axis, was flipped. In this step, the negative Z and Y coordinate point values are overturned, and the data are reset and arranged with Z values going low to high in a CSV file. Mentionable, this method can only be used for 3 axis milling operations.

Figure 11. A sample of the initial Cartesian coordinates (xyz) exported in excel sheet as. csv file.

3.3.2. Targeted G-code decoding with MATLAB

Both CAM and AM methods essentially send instructions in a coded format to the machine (CAM/CNC/3D printers) it needs to run, called the G-code files. As mentioned earlier, G-code files with detailed machine instruction not only contain the required nodal coordinates but also other necessary instructions specific to a CNC or a 3D printer to operate it safely. As such, G-codes vary depending on the machine it is programming to run, that is, machine-specific. Using MATLAB coding, specifically targeted to the G-Code type (CAM/AM), the tool end point coordinates are extracted discarding all other information. At first, only the ‘‘exhaust’’ component is studied, and then, the full powerplant model is used which gives a dense plot of $36525\times3$ double points, that includes all the tool end point coordinates (Fig. 12(a)).

Figure 12. Coordinates obtained from a G-Code (AM method) of the power plant using MATLAB.

3.3.3. Toolpath/Trajectory modification

Although the targeted programming successfully extracts the coordinates, the number points is very dense as is expected from a 3D printer slicer generated G-Code (Fig. 12(a)). There is also variation of point cloud density, for example, the flat horizontal surfaces have denser point clouds than the other parts. To reduce the number of points, to get rid of noise and also to get a uniform external trajectory pattern without any conflict with the plant surface mesh, MATLAB code is modified that gives an updated trajectory (Fig. 12(b)). However, the arbitrary jump from one layer to another give trajectory lines through the surface mesh itself. Also, the support generation at the bottom of the structure remains which is unnecessary.

Trajectory is further modified, and an improved trajectory is achieved that has significantly smaller number of points (373x3 double) as well as a conflict-free trajectory (Fig. 13). Also, the significant reduction of point cloud (almost 1/100th) will result in better trajectory and computationally less expensive when fed to the inspection UAV (Fig. 14).

Figure 13. Final modified trajectory the power plant obtained using MATLAB programming.

Figure 14. Significant reduction in the number of points from the original Slicer G-code (toolPath1) (Fig. 12(a)), a modified trajectory (toolpath2) (Fig. 12(b)) to the final trajectory generation (toolPath3) (Fig. 13) using MATLAB programming.