*Article* **Selection of the Depth Controller for the Biomimetic Underwater Vehicle**

**Michał Przybylski**

Faculty of Mechanical and Electrical Engineering, Polish Naval Academy, 81-127 Gdynia, Poland; m.przybylski@amw.gdynia.pl; Tel.: +48-261-262-524

**Abstract:** The aim of this paper is to select a depth controller for innovative biomimetic underwater vehicle drives. In the process of optimizing depth controller settings, two classical controllers were used, i.e., the proportional–integral–derivative (PID) and the sliding mode controllers (SM). The parameters of the regulators' settings were obtained as a result of optimization by three methods of the selected quality indicators in terms of the properties of the control signal. The starting point for the analysis was simulations conducted in the MATLAB environment for the three optimization methods on three types of indicators for three different desired depth values. The article describes the methods and quality indicators in detail. The paper presents the results of the fitness function obtained during the optimization. Moreover, the time courses of the vehicle position relative to the desired depth, the side fin deflection angles, the calculated parameters of the control signals, and the observations and conclusions formulated in the research were presented.

**Keywords:** biomimetic underwater vehicle; depth controller; genetic algorithms; particle swarm optimization; Pareto optimization

#### **1. Introduction**

In the 21st century, there has been dynamic development of mobile robots. One of the fields is underwater robots, where we can distinguish between ROV (Remote Operated Vehicle) and AUV (Autonomous Underwater Vehicle). Vehicles based on ROV technology, which dates back to the 1960s, have been ideally developed and are successfully used in all kinds of underwater operations [1–3]. A significant limitation of the mobility of these vehicles is the use of cable tether, which leads scientists and engineers to develop AUVs [4]. Autonomy, advanced control and positioning systems allow the realization of many civilian and military tasks. The development of bionics resulted in a new trend in the construction of underwater mobile robots, whose main idea is to imitate underwater animals. These vehicles are called biomimetic underwater vehicles (BUVs) [5,6], which mimic construction and motion kinematics. Most often, the prototype vehicle is inspired by the shape and movement of aquatic creatures, although there are designed based on manta rays, penguins and many others. Nevertheless, it is essential to carefully analyze the animal movement and to develop an appropriate simplified mathematical model, which will be used as close as possible to the BUV [7]. Artificial fins, similar in shape and appearance to the real ones, propel the BUV. Such propulsion is called wave or undulating propulsion, and it can be placed in different parts of the BUV's hull, depending on the type of design, its maneuverability and the speeds it can achieve. An electric motor usually drives it to generate sufficient vehicle thrust through oscillating motion. A sinusoidal function usually describes oscillatory motion, but different parts can be used depending on the type and purpose of the undulating propulsion. To describe oscillatory motion, usually use parameters such as the neutral position, which is the zero position for oscillation, the amplitude of oscillation, which defines the maximum deformation of the fins, and the frequency of oscillation [8].

**Citation:** Przybylski, M. Selection of the Depth Controller for the Biomimetic Underwater Vehicle. *Electronics* **2023**, *12*, 1469. https:// doi.org/10.3390/electronics12061469

Academic Editor: Hamid Reza Karimi

Received: 7 February 2023 Revised: 9 March 2023 Accepted: 14 March 2023 Published: 20 March 2023

**Copyright:** © 2023 by the author. 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/).

Control of underwater vehicles has been investigated quite thoroughly. However, still many difficulties occur, among other things, due to environmental disturbance, highly nonlinear behavior of vehicles, the complexity of the vehicle hydrodynamics or the application of new innovative propulsion systems. The articles [9–12] present various systems for controlling the depth of underwater vehicles. They concern both controllers tuned classically and using artificial intelligence algorithms. However, the works presented above concern, to a large extent, ROV or AUV vehicles using a classic screw drive propeller. Biomimetic underwater vehicles can be equipped with an artificial swim bladder [13] similar to a fish bladder, which controls buoyancy and depth. Another solution used to control the depth of biomimetic underwater vehicles [14,15] is changing the locomotion primitives of the fins, usually by changing frequency, amplitude or side fins' phase shift. In more advanced applications, changing the depth can be done by adjusting the angles of attack of the vehicle's control surfaces, their stiffness [16], or surface area of side fins [17] while the vehicle is moving at a certain speed through its wave propulsion.

The main contribution of this paper is to present depth controllers' different tuning methods for innovative BUV wave propulsion drives. The proposed methods use artificial neural networks in the tuning process and the experimental verification system of controller gains to evaluate its performance. In the presented literature, one can notice the lack of use of neural network-based methods for tuning depth controllers based on bioinspired wave propulsion drive. This research could have significant implications for the development of biomimetic autonomous underwater vehicles for various applications such as exploration of the ocean ecosystems, environmental monitoring, or even search and rescue missions. The rest of the paper is organized as follows. The mathematical model of the vehicle is based on the Fossen model [18] with modifications to include the new wave drive used. Commonly used underwater vehicles are driven mainly by a set of screw propellers, whereby the thrust generated by the propellers can be calculated using specific mathematical formulas [19]. In the case of a new wave propulsion system that mimics the action of fish fins, the forces and moments of force acting on the vehicle were calculated using the author's method presented in the literature [20]. Then the applied depth controllers and their control methods are presented. The following section offers the research problem and the results obtained in the simulation process. The final section formulates conclusions and plans for future research.

#### **2. Materials and Methods**

This chapter presents the simulation model and methods used in the tuning process of the depth controller settings. It contains two main subsections. The first concern describes biomimetic underwater vehicles with an innovative propulsion system. A mathematical model based on the above-mentioned vehicle is described in the further part of this subsection. While the second main subsection deal with the description of controllers, optimization methods and fitness functions used in the simulation process.

#### *2.1. Control Object*

This subsection describes the construction of a biomimetic underwater vehicle (mini-Cyber Seal) with particular emphasis on the new propulsion system. In the next part, the mathematical model of the vehicle is presented, where special attention is paid to the proposed solutions for modelling the vector of vehicle forces and moments.

#### 2.1.1. Mini CyberSeal

The study used a physical model of the BUV mini CyberSeal vehicle shown in Figure 1, is a smaller prototype of the larger vehicle. The purpose of the downsized BUV was to test the performance of the larger vehicle's counterpart's propulsion system and to test a new type of control. Unlike its predecessors, the mini CyberSeal's propulsion system had two tail fins instead of one. They generate the major thrust and are additionally responsible for changing course. In addition, a side propulsion system is mounted in front of the fuselage,

which generates additional thrust and changes the depth. The rear and side fins are made of polycarbonate and rubber.

**Figure 1.** Model of the BUV mini CyberSel [own source].

The rear fins can move at a frequency of up to 3 Hz within a range of about 80 degs outward, and up to 12 degs to the vehicle's inside axis of symmetry. The side fins can also move at up to 3 Hz over a range of ±45 degress. All electronic components, sensors, servos and a 7.2 V 10 Ah li-ion battery were enclosed in a sealed tube made of POM-C material. The servos with a 1:1 gearbox, responsible for the movement of the fins, are controlled via POLOLU-1353 miniMaestro servo controllers via RS232. The base station communicates via a WiFi network supported by a TP-Link TL-WR702N access point and a WIZNET Wiz-145SR server port providing four serial ports. The vehicle has an internal artificial buoyancy bladder, which is responsible for the vehicle's static depth adjustment or buoyancy control. The mini CyberSeal vehicle has been equipped with an OS-5000 digital compass from Ocean Server and an A-10 depth sensor from Wika with a measurement range of 0–1 bar. The sensors mentioned above make it possible to read current depth values and vehicle motion parameters, i.e., angle of an inclination concerning individual axes of the coordinate systems.

#### 2.1.2. Mathematical Model

The model captures the underwater vehicle's rigid body dynamics, hydrostatics and hydrodynamic effects. Critical issues for the modelling and simulation of BUVs are model complexity, ease of implementation and accuracy of prediction. Therefore, for the mathematical description, some simplifications are adopted for the vehicle: it has three planes of symmetry, moves at low speed in a viscous fluid, and has six degrees of freedom. When analyzing the motion of an underwater vehicle, two reference systems are defined:


The moving coordinate system is commonly called the "vehicle reference system", and its origin corresponds to the geometric center of the vehicle. The different axes of this coordinate system correspond to the following:


Changes in the position of the moving *xoyozo* coordinate system are described relative to the adopted *xyz* coordinate system associated with the Earth. Due to the low velocity of the vehicle, the acceleration of points on the Earth's surface due to its spin is neglected, and the *xyz* system is considered stationary. It is suggested that angular and linear velocities be described in the reference system associated with the vehicle while the vehicle's orientation is described in a stationary coordinate system. The quantities describing the vehicle's movement are defined according to the SNAME notation in Table 1 .

$$\mathbf{M}\dot{\boldsymbol{\nu}} + \mathbf{D}(\boldsymbol{\nu})\boldsymbol{\nu} + \mathbf{g}(\boldsymbol{\eta}) = \boldsymbol{\tau} \tag{1}$$

where:

*ν*—vector of linear and angular velocities, i.e., *ν* = [*u*, *v*, *w*, *p*, *q*,*r*];

*η*—vector of vehicle position and Euler angle coordinates in the stationary system;

**M**—inertia matrix (equal to the sum of the rigid body mass matrix **MRB** and the associated masses matrix **MA**);

**D**(*ν*)—hydrodynamic damping matrix;

g(*η*)—matrix of restoring forces (gravity forces *P* and buoyancy forces *B*);

*τ*—vector of forces and moments acting on the vehicle.

**Table 1.** Notation used in describing the movement of underwater vehicles.


Given the assumptions mentioned above, a nonlinear model of motion in six degrees of freedom is adopted for simulating the movement of the mini CyberSeal. The action of the vehicle is described by six differential equations, which, presented in matrix form, have the following format: The right side of Equation (1) represents the vector of forces and moments of force acting on the vehicle generated by the vehicle's propulsion system (2).

$$\pi = [X, Y, Z, K, M, N] \tag{2}$$

where:

*X*, *Y*, *Z* —forces acting on the vehicle in the longitudinal, transverse and vertical symmetry axis, respectively;

*K*, *M*, *N*—moments of forces acting in relation to the longitudinal, transverse, and vertical symmetry axis, respectively.

The vector of forces and moments of forces generated by the wave propulsion can be calculated by considering the propulsion system set-up in each design. Figure 2 shows the mini CyberSeal propulsion model consisting of two counter-phased tail fins and two independently controlled side fins. The thrust produced by each fin should be conveyed to the center of gravity *O* (Figure 2) using simple vector transformation formulas:

$$X = X\_{tl} + X\_{tp} + X\_l + X\_p \tag{3}$$

$$Y = Y\_{tl} + T\_{tp} \tag{4}$$

$$Z = Z\_1 + Z\_p \tag{5}$$

$$K = 0\tag{6}$$

$$M = M\_l + M\_p \tag{7}$$

$$N = N\_{tI} - N\_{tp} + N\_{I} - N\_{p} \tag{8}$$

where:

*tl*, *tp*, *l* and *p*—subscripts referring to the action of the left rear fin, right rear fin, left-side fin and right-side fin, respectively.

**Figure 2.** Mini CyberSeal propulsion model [own source].

The individual components of the vector, e.g., *Xtl*, *Ytl*, *Ntl* can be calculated using the position of these fins with respect to the centre of gravity according to the equations:

$$X\_{t\!I} = \cos(\beta\_{\!I}) \* T\_{t\!I} \tag{9}$$

$$X\_{tp} = \cos(\beta\_p) \ast T\_{tp} \tag{10}$$

$$Y\_{l\bar{l}} = \sin(\beta\_{\bar{l}}) \, \* \, T\_{l\bar{l}} \tag{11}$$

$$Y\_{tp} = \sin(\beta\_p) \* T\_{tp} \tag{12}$$

$$N\_{tl} = r\_2 \ast X\_t l + r\_1 \ast Y\_{tl} \tag{13}$$

$$N\_{tp} = r\_2 \ast X\_t p + r\_1 \ast Y\_t p \tag{14}$$

$$X\_l = \cos(\alpha\_l) \, \* \, T\_l \tag{15}$$

$$Z\_I = \sin(\alpha\_I) \ast T\_I \tag{16}$$

$$M\_l = r\_3 \ast Z\_l \tag{17}$$

$$N\_l = r\_4 \* X\_l \tag{18}$$

More information on the mathematical dependencies of the novel drive used is contained in the [21]. As shown in Figure 1 or Figure 3, the Mini CyberSeal has two side fins and two tail fins, which generate time-varying thrust.

**Figure 3.** Mini CyberSeal thrust measurement stand [own source].

The value of the thrust *Ttl*, *Ttp*, *Tl*, *Tl* depends on the control parameters, including the frequency and amplitude of the deflection of each fin. The thrust values generated by each fin for different frequencies and amplitudes were determined experimentally. Finally, the vehicle's speed is variable and dependent on the frequency of fin oscillations. The method of measuring the thrusts generated by the mini CyberSeal used in the mathematical equations is presented in the literature [20]. The thrust *T* generated by the fin is the sum of two components:

$$T = T\_{av} + T\_{osc} \tag{19}$$

where:

*Tav*—constant thrust component at a specific fin oscillation frequency;

*Tosc*—variable component modelled by a sinusoidal wave with a specific amplitude (at a specific fin oscillation frequency).

At the same time, the test stand is shown in Figure 3. In addition, the right side of Equation (1) considers the effect of environmental disturbances such as wind, waves and sea currents, which significantly impact the BUV. The left side of Equation (1) describes the forces and moments of force caused by physical phenomena, such as rigid body inertia and the inertia of masses accompanying viscous fluid, the hydrodynamic drag exerted by water, and the balance of gravity and buoyancy forces. Using the mathematical relationships included in the literature [18], the matrix parameters describing the left side of Equation (1) can be calculated.

#### *2.2. Methods*

This subsection deals with the main part of the work concerning the selection of depth controllers' settings. The first subsection describes the controllers used, i.e., PID and SM. The following subsections describe the optimization methods and the fitness functions based on which the aforementioned methods realized optimization of used controllers settings.

#### 2.2.1. Depth Controllers

Two classic controllers were used as depth controllers for the mini CyberSeal vehicle. The first one is about a PID controller that calculates the error value *e*(*k*) as the difference between the set depth value and the value received from the depth sensor and applies a correction based on proportional, integral and derivative terms (denoted P, I, and D, respectively). Hence the name [22,23]. The action of the PID controller is described by the following formula presented in the discrete form:

$$
\mu(k) = k\_P e(k) + k\_i \sum\_{k=1}^{k\_{max}} e(k) + k\_d \Delta e(k) \tag{20}
$$

where:

*u*(*k*) is a control signal in *k* step of simulation;


*kp*, *ki* and *kd* are constant quantities called gain factors.

The second controller used is the sliding mode controller (SM) [24,25], where sliding mode control is achieved by controlling nonlinear systems, which changes the dynamics of a nonlinear system by applying a discontinuous control signal, which forces the system to "slide" along upstream of the expected behavior of the system. It is calculated using the following formulas, also presented in the discrete form:

$$s(k) = \frac{\lambda e(k) + \Delta e(k)}{q} \tag{21}$$

$$\text{if} \quad |s(k)| > 1, \quad \text{then} \quad s(k) = \text{sign}(s(k)) \tag{22}$$

$$
\mu(k) = k\_s s(k) \tag{23}
$$

where:

*u*(*k*) is a control signal in *k* step of simulation;

*s*(*k*) is a normalized control signal in *k* step of simulation;

*e*(*k*) is an error signal in *k* step of simulation;

Δ*e*(*k*) is a change of error signal in *k* step of simulation, i.e., *e*(*k*) − *e*(*k* − 1);

*λ*, *ϕ*, *ks* is a constant settings of SM controller.

#### 2.2.2. Optimization Methods

To tune controller settings, used the Global Optimization toolbox of the MATLAB environment [26,27]. The following three optimization methods were used: Genetic Algorithm (GA), Particle Swarm Optimization (PSO) and Pareto Simulation (PSA).

The genetic algorithm, first formalized as an optimization method by Holland, is a global optimization technique for multi-dimensional, nonlinear, and noisy problems and a stochastic search technique based on the mechanism of natural selection and natural genetics [28,29]. A genetic algorithm (GA) solves both constrained and unconstrained optimisation problems based on a natural selection process that mimics biological evolution. The algorithm repeatedly modifies a population of individual solutions. At each step, the genetic algorithm randomly selects individuals from the current population and uses them as parents to produce the children for the next generation. Changes are introduced into the offspring through mutation, crossover and other genetic operators. The procedure ends when satisfactory genotypes (a set of traits of an individual) are obtained, which are matched by phenotypes with a high fitness function (an individual from the population). Over successive generations, the population "evolves" toward an optimal solution. An initial population of 40 individuals was generated using a MATLAB random generator during optimization. Individuals in the current generation are estimated using one of the three fitness functions described in the following subsection. After calculating the fitness function, the reproduction algorithm creates children for the next generation. The following operators are used in reproduction: fitness rank scaling, Stochastic uniform selection function, Crossover fraction equal to 0.8, and Gaussian mutation function. The GA optimization was stopped when the maximum number of 100 generations was reached and/or when no change in the best fitness function value for new generations was detected during the next 50 steps [30].

The PSO algorithm is based on a simplified social model closely tied to swarming theory. It solves a problem by having a population of candidate solutions, here dubbed particles, and moving them around the search space according to simple mathematical formula over the particle's position and velocity. Each particle's movement is influenced by its local best-known position. Still, it is also guided toward the most notable positions in the search space, updated as better places are found by other particles [31]. A physical analogy might be a swarm of bees searching for food sources. In this analogy, each bee (referred to as a particle here) uses its memory and knowledge gained by the swarm to find the best available food sources. This is expected to move the swarm toward the best solutions. Based on the literature [32,33], the following PSO parameters were assumed: (1) MaxStallIterations (relative change in the value of the best objective function): 20, (2) MinNeighborsFraction (setting both the initial neighborhood size for each particle and the minimum neighborhood size): 1, (3) SwarmSize: 200. As in the case of GA, one of the critical problems is to properly define the fitness function to get the correct optimization rates. The fitness functions used were analogous to GA and are presented in the following subsection.

Pareto optimization (PSA) is a field of multi-objective decision-making that deals with mathematical optimization problems involving more than one objective function for simultaneous optimization. Multi-objective optimization has found application in many scientific areas, including engineering, economics and logistics, where optimal decisions must be made during trade-offs between two or more conflicting objectives. When we have several objective functions that we want to optimize simultaneously, these solvers find optimal trade-offs between competing objective functions. This method can also be applied to a single-objective problem. PSA uses pattern search on a set of points to iteratively search for non-dominated points. It should fulfil all constraints and linear constraints in each iteration. Theoretically, the algorithm converges to points near the true Pareto front [34]. In the algorithm used, in the beginning, the PSA creates an initial set of 200 randomly selected points and then checks whether these points are feasible concerning the bounds and linear constraints. If impossible, the algorithm projects the initial points into a linear subspace of linearly feasible points by solving a linear programming problem and removing

duplicate points. The PSA then divides the points into two sets named "archive" and "iterative". The archive set contains non-dominated points associated with a mesh size less than 10−<sup>6</sup> and satisfying all constraints within 10<sup>−</sup>6. PSA checks the location of each point in the 'iterative' set. Success is achieved if the polled points yield at least one dominated point. PSA then extends the probing in successful directions multiple times, doubling the 6e grid to find a dominant point. If any non-dominated point is obtained, the grid size is halved. The algorithm stops when: (1) the mesh size exceeds the value (+Unity), (2) the fitness function decreases to the value (-Unity), (3) it reaches the maximum number of iterations equal to 400.

#### 2.2.3. Fitness Function

It is necessary to determine proper fitness functions to obtain appropriate optimization results for the depth controller. For this paper, three functions commonly used in mathematics were formulated [35,36]. The first is Integral Absolute Error (ISA), the sum of the absolute values of the error signals *e*(*k*) in all simulation steps. Its task is to select the controller parameters so that the error rate between the desired depth and the present depth value is as low as possible throughout the simulation period. The ISA in discrete form is shown in Equation (24).

$$f\_{fit1} = \sum\_{k=1}^{k\_{max}} |e(k)|\tag{24}$$

The second fitness function is an Integral of Squared Error (ISE) in the discrete form Equation (25). The ISE integrates the square of the error over time. ISE will penalize for significant errors more than smaller ones (since the square of a large error will be much bigger). Control systems specified to minimize ISE will tend to eliminate significant errors quickly but will tolerate minor errors persisting for an extended period. Often this leads to fast responses but with considerable, low amplitude oscillation.

$$f\_{fit2} = \sum\_{k=1}^{k\_{max}} e(k)^2\tag{25}$$

The third proposed fitness function is based on combining two direct control quality indexes and is presented in the following Equation (26).

$$f\_{fit3} = \sum\_{i=1}^{i\_{max}} t\_r(i) + k\_M \sum\_{i=1}^{i\_{max}} M\_p(i) \tag{26}$$

It takes into consideration rising times *tr* in [s] and first overshoots *Mp* in [rad] for all *imax* changes of the desired course. Because of the small value of *Mp* in [rad] compared to *tr* in [s] additional gain factor of the sum of first overshoots was introduced *kM* = 25.

#### **3. Research Problem and Results**

The research problem of this work is to find a solution for the most effective depth control of the mini CyberSeal vehicle. With a mathematical model of the vehicle, two types of controllers (PID, SM), three different methods for optimizing controller settings, and three fitness functions, an attempt was made to select a suitable controller and its parameters. All regulators, methods and functions should be examined in the tuning process and verified for each combination (Figure 4). The tuning process was carried out for three different immersion changes: (1) 0.2 [m]—shallow immersion, (2) 0.5 [m]—deep immersion, (3) and then two immersion changes—first 0.5 [m] then second 0.2 [m] depth change after 20 s of simulation. In contrast, the verification process was based on the tuning process for each combination of 25 randomly selected depths from 0 to 0.7 [m] (in 0.1 [m] increments). The average value of the fitness function was conducted from 25 verification tests.

**Figure 4.** Flowchart optimization process to receive the best controllers settings.

#### *Results and Discussion*

The test results are shown in Tables 2–4 for the corresponding fitness functions 1–3. The table includes three different combinations of depth change for both the testing (T) and verification (V) processes. Also included in the table are combinations of controller type and its optimization method, e.g., SM-GA means SM controller was optimized by GA method, PID-PSA means PID controller was optimized by PSA method, etc. The initial study's objective was to find the optimal barriers for the controllers to avoid going into local minima. The research was carried out by an expert using a designed model of the CyberSeal, designed depth controllers and various tuning methods. As a result of the preliminary study, the upper and bottom barriers assume the following values: (a) for PID controller [*kp*, *kd*, *ki*]—bottom barrier equals [200, 10,000, −1] and upper barrier equals [600, 40,000, 1], (b) for SM controller [*λ*, *ks*, *ϕ*]—bottom barrier equals [−5, −100, −2] and upper barrier equals [5, 100, 2]. The main objective of the primary research was to compare two classic controllers, three methods optimized by fitness functions in response to changing the desired immersion value. The only quality control criterion optimized during tuning depth controllers was the minimized value of three different fitness functions. It means that all

phrases describing the "best results", "most effective", etc., which are used later in the article, refer to the smallest fitness function values.

Analyzing the results contained in Table 2 obtained for fitness function no. 1, the following conclusions can be drawn: (1) in the verification process, the best results for both controllers and all methods were obtained for a large change in immersion value, (2) in the tuning process, the best results for both controllers for all methods were received for a slight change in the immersion value, (3) PID-GA obtained the best result in the verification process. In contrast, SM-PSO and SM-PSA got the best result in the tuning process considering all of the depth changes, (4) optimization methods present similar efficiency for all of the desired depth changes, (5) the smallest values of the *f fit*<sup>1</sup> were obtained during tuning and verification for the SM controller than for the PID controller.

**Table 2.** Values of fitness function no. 1 for tuning (T) and verifying (V) BUV's depth controllers for three changes of desired immersion: shallow, deep and two following depth changes.


Considering the results shown in Table 3 for fitness function no. 2, the following conclusions can be formulated: (1) in the verification process, the best results for both controllers and all methods were obtained for two changes in immersion value, (2) similar to earlier, the best results for both controllers for all methods in the tuning process were obtained for a slight change in the immersion value, (3) PID-GA obtained the best result in the verification process, while SM-PSO and SM-PSA obtained the best result in the tuning process taking into account all of the depth changes, (4) optimization methods present similar efficiency in the tuning process for all of the desired depth change.

**Table 3.** Values of fitness function no. 2 for tuning (T) and verifying (V) BUV's depth controllers for three changes of desired immersion: shallow, deep and two depths changes.


Table 4 shows the results for fitness function no. 3. Analyzing its results, one can deduce: (1) as earlier, the best tuning process of both controllers was obtained for shallow immersion, while the best results for the verification process were received for two depth changes, (2) the smaller *f fit*<sup>3</sup> was obtained during tuning and verification process by PID controller than SM, (3) PID-PSO received the smallest value of *f fit*<sup>3</sup> in tuning for shallow immersion and verifying process for two depth changes, (4) optimization methods present similar efficiency for all of the desired depth change.


**Table 4.** Values of fitness function no. 3 for tuning (T) and verifying (V) BUV's depth controllers for three changes of desired immersion: shallow, deep and two following depths changes.

Analyzing the above results in Tables 2–4, it can be concluded that the best results were obtained by PID-GA using objective function no. 3 for two depth changes. The result of *f fit*<sup>3</sup> indicates that the selected controller settings in the optimization process performed best in the verification process for random desired depth (for 25 random changes in depth from 0 to 0.7 m). From Figures 5–7 show the simulation results for the controller settings obtained for the PID-GA method, using the *f fit*<sup>3</sup> function and two depth changes. Figures 5–7 show, respectively, the results for each desired depth, i.e., 0.2 m, 0.5 m and two depth changes, first of 0.5 and then of 0.2. The simulation was carried out for 35 s. Considering the third case, the second desired depth signal occurred in 20 s of simulation. The timing was chosen so that the following depth changes occurred when there were no fluctuations in depth after the first change. The depth change was realized by the parallel swing of the side fins by the angle set by the control signal. In all the simulations, the oscillation frequency of the vehicle's side fins was constant and equal to 2 [Hz]. In the simulation, the side fin deflection angles were limited to 60 [deg] to avoid too excessive a vehicle trim angle. Each figure shows plots of the various parameters as a function of time: (1) graph of the dependence of the current vehicle depth on the desired, (2) values of the angular deflections of the side fins responsible for changing the depth of the vehicle based on the control signal along with the value of the vehicle trim (pitch), (3) error of depth over time. Analyzing the graphs presented, it can be noticed that selected controller gains are appropriate, and their selection is important for achieving good performance and stability of the system. The above graphs show that the object obtains satisfactory stability not only for the depths for which it was optimized (Figure 7) but also for other desired depths (Figures 5 and 6). The controller, in each case, reaches the required immersion quickly, and there is minimal overshoot. The controller's gains are not too large, so the system avoids exhibiting chattering, which is rapid switching between different control modes, which can lead to wear and tear on the actuators. The vehicle achieves stability within 12 s for all desired depths. To confirm the correctness of the adopted solution, i.e., tuning the controller for a specific value and then checking its settings for 25 random values, the simulation results for other controller settings are presented in Figure 8. For this figure, the controller settings were obtained for the PID-GA method, using the *f fit*<sup>3</sup> function and small depth changes. Comparing the simulation results for the same depth change (in Figures 5 and 8), but obtained for two different controller settings, it can be observed that the controller tuned for the specific depth reached the desired depth faster, i.e., 8 s (Figure 8). This is confirmed by the results in Table 4, wherein the tuning process for a small change in depth, the PID-GA obtained the value of the fit function 8.21, while for two changes in depth is 20.1. However, better results were obtained in the verification process for two depth changes for the entire depth spectrum (from 0 to 0.7 m). The value of *f fit*<sup>3</sup> for two depth changes was better by 50 per cent than for a small depth.

**Figure 5.** Changes of immersion and side fins deflection in time in response to desired depth: 0.2 [m] obtained for PID controller settings using GA method and fitness function no. 3 for two depth changes.

**Figure 6.** Changes of immersion and side fins deflection in time in response to desired depth: 0.5 [m] obtained for PID controller settings using GA method and fitness function no. 3 for two depth changes.

**Figure 7.** Changes of immersion and side fins deflection in time in response to subsequent desired depths: 0.5, 0.2 [m] obtained for PID controller settings using GA method and fitness function no. 3 for two depth changes.

**Figure 8.** Changes of immersion and side fins deflection in time in response to desired depth: 0.2 [m] obtained for PID controller settings using GA method and fitness function no. 3 for small immersion.

#### **4. Conclusions**

This paper presents the tuning process of two classical depth controllers (PID, SM) of the biomimetic underwater vehicle with undulating propulsion. The tuning of the parameters of the controllers was carried out by applying three optimization methods and three fitness functions (GA, PSO, PSA), which provided a quality criterion. The optimization process was conducted for three different desired depth values, as shown in Figure 4.

Summarizing the research results obtained, several conclusions can be reached. The best results were not received with a slight change in immersion value. The results shown in Tables 2–4 confirm that despite the best values of the fitness function in the tuning process, the values for the verification process are several times higher. It can be recognized that the values of controller setting gains work efficiently only for small values. At the same time, they are not necessarily optimal for larger values used in the verification process. Also, the values of the fit function for different controller and optimization methods for the tuning and verification process are similar. Therefore, it is expected that the tuning of the controllers should be carried out for a larger desired depth or two or more changes of depth.

Furthermore, it can be deduced that, in most cases, better control quality was obtained for the PID controller than for the sliding controller. The differences can be seen in Table 4, where in all cases, both in the simulation and verification process, the fitness function obtained better values for the PID controller. The three optimization methods used to adjust the controller settings achieve comparable efficiency. Although the PSO and PSA methods achieved better tuning, the GA method achieved the best results in the verification process for each of the three different fitness functions.

Comparing the quality indicators, we can see that the ISE (*f fit*2), which strongly penalizes any large deviation, obtained the worst results for a large change or two depth changes. In contrast, the best results were obtained for fitness function no. 3 (*f fit*3) where used direct quality indicators, i.e., the rise time and the value of the first overshoots. It is because the fitness function obtained the most repetitive results for the tuning and verification process, and the values of fitness function no. 3 were the best compared to the other functions (*f fit*1) and (*f fit*2) using the classic indicators.

The presented results could be more comprehensive. An interesting issue would be further systematic testing of new multi-criteria indicators. In this way, it would be possible to use, for example, an LMS indicator or the like to limit sudden changes in the control signal and a more aggressive indicator such as ISE or ISE to increase the speed and accuracy of the process control. As part of future research, it is also planned to implement the applied regulators on a real object (mini CyberSeal) and verify them in natural conditions.

**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 author declares no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:



#### **References**


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

**Edyta Ładyzy ´ ˙ nska-Kozdra´s 1,\*, Anna Sibilska-Mroziewicz 1, Krzysztof Sibilski 2, Danyil Potoka <sup>1</sup> and Andrzej Zyluk ˙ <sup>2</sup>**


**Abstract:** Most aircraft launchers exhibit a rapid acceleration of the launching aircraft, often exceeding ten times the acceleration due to gravity. However, only magnetic launchers offer flexible control over the propulsion force of the launcher cart, enabling precise control over the aircraft's acceleration and speed during its movement on the launcher. Consequently, extensive research is being conducted on magnetic launchers to ensure the repeatability of launch parameters, protect against aircraft overloads, and ensure operator safety. This article describes the process of modeling and analyzing the dynamical properties of a launch cart of an innovative prototype launcher, which employs a passive magnetic suspension with high-temperature superconductors, developed under the GABRIEL project. The developed mathematical model of the magnetic catapult cart was employed to conduct numerical studies of the longitudinal and lateral movement of the cart, as well as the configuration of the UAV–cart system during UAV takeoff under variable atmospheric conditions. An essential aspect of the research involved experimentally determining the magnetic levitation force generated by the superconductors as a function of the gap. The results obtained demonstrate that the analyzed catapult design enables safe UAV takeoff. External factors and potential vibrations resulting from uneven mass distribution in the UAV–cart system are effectively balanced by the magnetic forces arising from the Meissner effect and the flux pinning phenomenon. The primary advantage of the magnetic levitation catapult, in comparison to commercial catapults, lies in its ability to provide a reduced and consistent acceleration throughout the entire takeoff process.

**Keywords:** magnetic launcher; levitation force; UAV; dynamics modeling; numerical simulation

#### **1. Introduction**

The designers of modern unmanned aerial vehicles (UAVs) must solve a series of technical problems related to the critical aspects of their operation in order to meet the needs of the rapidly developing market. One of the key aspects is the UAV takeoff procedure. The takeoff of a large portion of fixed-wing UAVs requires separate devices called launchers or aircraft catapults. Currently, rocket systems, rubber, pneumatic, and hydraulic launchers are most commonly used for this purpose. They are powered by various means, including steam, compressed air, energy stored in elastic elements (such as rubber bands), or even rocket boosters. The operation of these devices involves high acceleration for the launching UAV, often exceeding 10 times the Earth's gravitational acceleration, and a lack of control over the aircraft's trajectory during the launch. Magnetic launchers, which allow for significantly higher speeds compared to traditional catapult solutions while also enabling contactless operation, are an attractive alternative to the currently used launchers.

The first interest in electromagnetic catapults came from the navy, and during World War II, Westinghouse constructed the initial aircraft catapult system utilizing electromagnetic interactions [1]. However, due to the high costs of production and maintenance, as

**Citation:** Ładyzy ´ ˙ nska-Kozdra´s, E.; Sibilska-Mroziewicz, A.; Sibilski, K.; Potoka, D.; Zyluk, A. Dynamics of ˙ Separation of Unmanned Aerial Vehicles from the Magnetic Launcher Cart during Takeoff. *Electronics* **2023**, *12*, 2883. https://doi.org/10.3390/ electronics12132883

Academic Editors: Piotr Szymak, Stanisław Hozy ´ ˙ n and Paweł Piskur

Received: 13 May 2023 Revised: 26 June 2023 Accepted: 27 June 2023 Published: 29 June 2023

**Copyright:** © 2023 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/).

well as their large size and weight, these systems were not commercially viable. Electromagnetic technology was displaced by steam launchers in the past, and it was only at the turn of the century that the concept was revisited. Scientists were prompted by the miniaturization of electrical devices, which allowed for the construction of increasingly cheaper and smaller equipment, thus reducing the size and operating costs of electromagnetic launchers. There are several types of magnetic launchers, classified based on their design and motion generation methods. Rail launchers and coil launchers are particularly noteworthy. Among coil launchers, synchronous and asynchronous (inductive) launchers can be distinguished based on the method of generating magnetic force. Such launchers have been employed in recently commissioned US Navy aircraft carriers and in high-speed MAGLEV trains like the Inductrack. Rail launchers can utilize linear induction motors as the propulsion for the launch cart. Rail launchers are typically composed of two conductive rails made of electrically conductive material. A movable connector, also made of a conductor, is placed between the rails. The contact between the two rails and the connecting element must ensure good current flow. As a result of electric charge flow, following Ampere's and Biot–Savart's laws, a magnetic field is generated around the rails, which interacts with the charges moving in the connector. The movable element experiences a magnetic force that pushes it along the rails. The method described above is used in the high-speed MAGLEV trains of the Inductrack type. However, the Inductrack technique has a drawback as levitation forces occur after reaching a certain minimum speed. Additionally, due to their cost, they are difficult to accept in unmanned aircraft launch solutions. The solution presented in this article is significantly simpler and more cost-effective than the Inductrack approach.

The subject of the presented research is an innovative prototype of a catapult utilizing passive magnetic suspension with high-temperature superconductors [2,3]. The catapult prototype was constructed as part of the GABRIEL project (Integrated Ground and Onboard system for Support of the Aircraft Safe Take-off and Landing) [4,5]. The design of the prototype catapult features magnetic tracks attached to a stationary base, a levitating takeoff platform, and a linear drive (Figure 1a). During the UAV takeoff, it is mounted on the takeoff platform, which is an integral part of the cart made of a duralumin frame and supported on four supports (Figure 1b), in which high-temperature YBCO superconductors are placed (Figure 1c). After cooling the YBCO with liquid nitrogen, the takeoff cart lifts off and levitates above the source of the magnetic field.

**Figure 1.** Elements of a magnetic catapult utilizing the Meissner effect: (**a**) magnetic tracks with a launch platform installed on them; (**b**) a levitating support for the launch cart; and (**c**) hightemperature superconductors placed in the support.

One undeniable advantage of this solution is its simple design, which facilitates the transportation and installation with just two personnel required for operation. Moreover, this solution incorporates all the advantages of electromagnetic launchers, allowing for the adjustment of the launch cart speed and acceleration to minimize the applied overload on the aircraft during takeoff. An additional benefit of this solution is the ability to

utilize the system for the landing of the UAV, a feature that conventional or coil launchers cannot guarantee.

At the end of the GABRIEL project, tests were carried out on the designed unmanned aerial vehicle in Aachen for takeoff and landing. The tests were successful [6], but technical issues were encountered, such as the slowly damped vibrations of the cart and the loss of its stable position after disconnecting the linear drive. These issues served as inspiration for further research on the innovative magnetic catapult project.

This article presents numerical studies of the motion of the cart and the catapult system during the UAV launch under varying atmospheric conditions based on a developed model of the dynamics of the launch cart. An important aspect of the research was the experimental determination of the magnetic levitation force as a function of the gap generated by the superconductors.

#### **2. The Model of the Launch Cart for a Magnetic Launcher**

In this study, a mathematical model of the launch cart for the magnetic launcher was developed. This model takes into account key factors, such as the characteristics of the environment and the individual components of the launcher, which include the magnetic tracks and the cart that launches the UAV. To further describe and understand the system, various coordinate systems were introduced, which are shown in Figure 2.

**Figure 2.** Model of the magnetic launcher system along with the adopted reference systems.

An inertial system (Ofxfyfzf) was firmly tied to the surface of the Earth, providing a fixed reference point. The magnetic system (Omxmymzm) is used to describe the distribution of the magnetic field and the levitation force. This system is aligned with the symmetry axis of the launcher tracks. The cart system (Osxsyszs) is aligned with the longitudinal symmetry axis of the cart and illustrates the direction of the cart's movement. More detailed descriptions of these reference systems, along with a thorough analysis of the physical model and external influences, can be found in the referenced paper [7].

The dynamics of the launch cart for a magnetic launcher, which moves in three-dimensional space, is described by a system of mutually coupled nonlinear ordinary differential equations. These equations were derived for the launch cart modeled as a rigid body with six degrees of freedom, using the principle of momentum and angular momentum conservation. The mathematical model of the launch cart's motion dynamics was thoroughly developed in [8], and its main stages were discussed in the article [9]. The formalism for writing the equations of motion is based on a precisely defined method of indexing vector quantities and the properties of rotation matrices, particularly the representation of the derivative of the rotation matrix using a skew-symmetric matrix. The final form of the equations of motion is written in a local, moving reference frame associated with the launch cart Osxsyszs. The matrix form of the equations of motion is given by Formula (1). The first equation describes the relationship between the three-dimensional vector of forces acting on the system and the vector of change in the body's momentum. The second equation maps the vector of force moments to the vector of change in angular momentum.

$$
\begin{bmatrix} m\_{(s)} \mathbb{I}\_{3 \times 3} & -m\_{(s)} \mathfrak{s} \begin{pmatrix} r^s \\ g\_{s'} \end{pmatrix} \\ m\_{(s)} \mathfrak{s} \begin{pmatrix} r^s \\ g\_{s'} \end{pmatrix} & \mathfrak{I}^s\_{(s)s} \end{bmatrix} \begin{bmatrix} d\_{s/f}^{s \langle f \rangle} \\ \varepsilon\_{s/f}^{s \langle f \rangle} \end{bmatrix} - \begin{bmatrix} -m\_{(s)} \mathfrak{s} \begin{pmatrix} \omega r\_{s/f}^s \\ \omega r\_{s/f}^s \end{pmatrix} \mathfrak{s} \begin{pmatrix} r\_{s/f}^s \\ g\_{s/f} \end{pmatrix} \omega\_{s/f}^s \\ \mathfrak{s} \begin{pmatrix} \omega\_{s/f}^s \end{pmatrix} \mathfrak{I}\_{(s)s}^s \omega\_{s/f}^s \end{bmatrix} = \begin{bmatrix} F\_{(s)}^s \\ M\_{(s)s}^s \end{bmatrix} \tag{1}
$$

where: <sup>I</sup>3×<sup>3</sup> is a 3 <sup>×</sup> 3 identity matrix, *<sup>m</sup>*(*s*) is the mass of the launch cart, <sup>I</sup>*<sup>s</sup>* (*s*)*<sup>s</sup>* is the inertia tensor of the cart, s *rs gs*/*s* is the skew-symmetric tensor of the displacement vector of the cart's center of mass, *a s*(*f f*) *<sup>s</sup>*/ *<sup>f</sup>* is the linear acceleration of the cart, *ε s*(*f*) *<sup>s</sup>*/ *<sup>f</sup>* is the angular acceleration of the cart, *ω<sup>s</sup> <sup>s</sup>*/ *<sup>f</sup>* is the angular velocity of the cart, *<sup>F</sup><sup>s</sup>* (*s*) and *<sup>M</sup><sup>s</sup>* (*s*)*<sup>s</sup>* are the forces and moments acting on the cart.

The launch cart levitating above magnetic tracks is subject to magnetic forces and moments of force *F<sup>s</sup>* (*Ls*) , *M<sup>s</sup>* (*Ls*) , aerodynamic forces and moments of force *F<sup>s</sup>* (*As*) , *M<sup>s</sup>* (*As*)*s* , gravitational forces and moments of force *F<sup>s</sup>* (*Gs*) , *M<sup>s</sup>* (*Gs*) as well as the propulsion force from the linear motor *F<sup>s</sup>* (*Ns*) ; therefore, the resulting force vector and moment of force have the following form

$$F\_{(s)}^{s} = F\_{(G\_s)}^{s} + F\_{(L\_s)}^{s} + F\_{(As)}^{s} + F\_{(N\_s)}^{s} \tag{2}$$

$$M\_{(s)s}^{s} = M\_{(G\_{\bullet})s}^{s} + M\_{(L\_{\bullet})s}^{s} + M\_{(As)s}^{s} \tag{3}$$

The mathematical model of the cart should be supplemented with kinematic relationships describing its angular velocity *ω<sup>s</sup> <sup>s</sup>*/ *<sup>f</sup>* and linear velocity *v s*(*f*) *<sup>s</sup>*/ *<sup>f</sup>* , as well as a geometric relationship describing the vector between the local coordinate system of the cart and the inertial coordinate system *r f <sup>s</sup>*/ *<sup>f</sup>* .

$$
\omega\_{s/f}^s = \int \left( R\_f^s \varepsilon\_{s/f}^{f(f)} \right) dt \tag{4}
$$

$$w\_{s/f}^{s(f)} = \int \left( a\_{s/f}^{s(f)} - \mathfrak{s} \left( \omega\_{s/f}^{s} \right) v\_{s/f}^{s(f)} \right) dt \tag{5}$$

$$r\_{s/f}^f = \int \left(\mathcal{R}\_s^f v\_{s/f}^{s(f)}\right) dt\tag{6}$$

where *R<sup>s</sup> <sup>f</sup>* is the rotation matrix describing the orientation of the inertial frame relative to the local frame of the cart.

The described mathematical model of the launch cart provides a basis for further research on its dynamic properties, including the numerical analysis of its motion. These studies were preceded by an analysis of the acting forces. The process of the theoretical modeling of forces and moments acting on the launcher system was discussed in detail in [7]. The analyses presented there indicated interactions between the loads and their influence on the UAV launch process. The key role is played by the magnetic levitation force. Its experimental identification process is included in Section 3.

#### **3. Experimental Identification of Magnetic Levitation Force**

The levitation force, which lifts the launcher cart above the magnetic tracks, is caused by the Meissner effect [10]. The Meissner effect is a phenomenon in which a superconductor expels magnetic field lines from its interior. It has been studied in various applications, including Maglev trains and energy storage systems [11,12], as well as in the context of superconducting materials used in aerospace engineering, such as in the development of high-performance motors and generators [13].

In the analyzed magnetic launcher, the levitation force is a result of the Meissner effect occurring between the tracks of the magnetic rails and the superconductors enclosed in containers—supports of the launch cart. After filling the supports with liquid nitrogen, the superconductors transition into the superconducting state, and as a result of the Meissner effect, they begin to levitate above the tracks.

To identify the value of magnetic levitation force, a measurement experiment was carried out using the MTS Bionics testing machine. A section of the magnetic track rail was attached to the lower, stationary jaws of the machine, and a support with four YBCO superconductors was mounted on the upper jaws (Figure 3).

**Figure 3.** Experimental setup for measuring magnetic levitation force.

Figure 4 presents the results of measurements aiming to determine the effect of the levitation gap height on the value of the magnetic levitation force. In the first series of measurements (red line), the superconductors were cooled with liquid nitrogen at a height of 20 mm above the magnets, where the magnetic field has zero value (ZFC—zero field cooling). Such a superconductor is only affected by the magnetic force resulting from the Meissner effect. In the second series of measurements, the superconductors were cooled with liquid nitrogen in a non-zero magnetic field at a height of 3 mm above the magnets (dFC = 3 mm). The measurements were repeated twice (blue and green lines), which showed the convergence of the obtained results. The upper graph shows the levitation force value registered by the force sensor of the strength machine, while the lower graph shows the size of the levitation gap, i.e., the distance between the levitating support and the launcher rails.

The conducted observations of the superconductor's behavior in a magnetic field showed that the levitation force value is a function of the levitation gap (distance of the superconductor from the magnetic field). Additionally, in the case of filling the superconductor with liquid nitrogen in a non-zero magnetic field, the maximum recorded value of the levitation force is almost half the value obtained when filling the superconductor in the zero magnetic field. This confirms the fact that there are two forces acting on the superconductor: a force repelling the superconductor from the source of the magnetic field, resulting from the Meissner effect *FLM*, and a force attracting the superconductor to the source of the magnetic field, resulting from the flux pinning effect *FLP*.

$$F\_{\rm Ls} = F\_{\rm LM} + F\_{\rm LP} \tag{7}$$

**Figure 4.** Results of the measurements illustrating the effect of the levitation gap (**lower graph**) on the value of the magnetic levitation force (**upper graph**) for different liquid nitrogen flooding heights (dFC = 3 mm; ZFC—zero field cooling).

As a result of the flux pinning phenomenon, the superconductors placed inside the supports of the cart become magnetized [14]. The widely used advanced mirror-reflection method [15] for modeling magnet–superconductor interactions suggests that the levitation force can be expressed by the following equation [16].

$$F\_{LP} = \boldsymbol{m} \cdot \nabla \mathcal{B} \tag{8}$$

where *B* is the external magnetic field and *m* is the vector of the magnetic moment frozen within the superconductor at the moment of the phase transition.

The approximating functions for the levitation forces used in the numerical model are presented in Figure 5a. The green and red lines indicate the experimentally measured levitation force. The green line corresponds to the superconductor being cooled with liquid nitrogen in a zero magnetic field (*FLP* = 0), while the red line represents the case where the superconductor is cooled with liquid nitrogen at a height of 3 mm. Thus, the red line represents the sum of the force resulting from the Meissner effect (green line) and the force arising from the flux pinning phenomenon.

**Figure 5.** (**a**) Forces acting on the support of the cart with superconductors; and (**b**) spatial distribution of the vertical component of the magnetic field gradient above the launcher tracks.

At a height of 5.2 mm, the force resulting from the Meissner effect (with a value of 5.4 N) and the force resulting from flux pinning balance each other. By determining the value of the magnetic field gradient along the vertical axis at a height of 5.2 mm (with an absolute value of ∇*B* = 20.6 N/T) (Figure 5b) based on finite element analysis [12], the induced magnetic moment in the superconductor was calculated to be *m* = 0.26 T using Equation (8). Thus, knowing the value of *m*, the force resulting from flux pinning was approximated when the superconductors were immersed in liquid nitrogen at a height of 3 mm (purple line).

The light blue line was obtained by summing the measurements of the levitation force when the superconductor was immersed in liquid nitrogen under zero magnetic field conditions (green line) and the approximated force resulting from flux pinning (purple line). Similarly, the force acting on the superconductor immersed in liquid nitrogen at any height can be approximated using a similar approach.

The identification of the magnetic levitation force provided a basis for further numerical studies of the magnetic launcher system. The results of the numerical simulations of the longitudinal and lateral motion of the cart as well as the configuration of the UAV–cart system during UAV launch under variable atmospheric conditions are presented below.

#### **4. Numerical Studies of the Dynamics of the Launch Cart in a Launcher**

The computer program designed for the analysis of dynamics and the visualization of the trajectory of the launch cart was implemented in MATLAB. The simulation model was formulated as an initial value problem of ordinary differential equations based on the mathematical model discussed above.

Two simulation scenarios were demonstrated, encompassing:


The longitudinal motion of the launch cart is generated by the linear motor. In the GABRIEL technology demonstrator, which was designed as part of the project, the motor stator is mounted in the middle of the magnetic tracks, and its structure restricts the lateral movements of the cart [17]. Therefore, in the simulation analyses of the investigated launcher system, the omission of the lateral movements of the launch cart is fully justified. In the second simulation, the behavior of the cart was additionally analyzed when it is detached from the motor stator, taking into account the effects on its motion in space caused by the trapping of magnetic flux inside the superconductors.

The authors' intention was to conduct numerical analyses based on a nonlinear model of the launch cart's dynamics, described by Equations (1)–(6). Therefore, the analysis of these simulation scenarios is a result of applying initial conditions and physical relationships derived from modeling the forces acting on the system.

In both presented simulation scenarios, it was assumed that the launcher tracks are in a horizontal position. The simulation time was 2 s. During the first second, the launch cart is driven by a constant force of 20 N, while during the second s, the cart is braked and the propulsion force changes direction.

#### *4.1. Longitudinal Motion of the Cart*

In the presented simulation scenario, the behavior of the cart was investigated, where its center of mass is consistently located on the axis of symmetry of the tracks, but it does not coincide with the geometric center. It was assumed that the center of mass of the cart was displaced along the *xs* longitudinal axis by *rgs*/*s*}*xs* = 0.1 m. At the initial time, the cart is located at a height of {*rs*/*m*}*zm* = 6.6 mm. Due to the displacement of the center of mass, the cart is subjected to the gravitational moment of force. The unbalanced moment of force causes the cart to tilt around the *ys* lateral axis. The difference in height between the front and rear supports of the cart results in a different value of levitation forces exerted on individual supports. As a result, a moment of force is created that balances the gravitational moment of forces and stabilizes the orientation of the cart. Figure 6 shows the trajectory of the cart. The blue line represents the motion of the center of the cart, the purple line marks the position of the launch rails, while the green and red lines visible on the left and right reflect the position of the cart supports. It can be seen that the cart has tilted forward and the front supports are lower than the rear supports.

**Figure 6.** The trajectory of the longitudinal motion of the launch cart.

After 2 s, the cart traveled a distance of 14.2 m (Figure 7a) and reached a maximum speed of 14 m/s, and then came to a stop (Figure 7b). The simulation showed small vibrations of the cart (Figure 8a) around the equilibrium position. The vertical velocity of the cart oscillates within ±0.0035 m/s (Figure 8b). The tilt angle of the cart stabilizes at −0.005 rad (Figure 9a), and then its angular velocity decreases to zero (Figure 9b).

**Figure 7.** The longitudinal component of position (**a**) and velocity (**b**) of the launch cart.

**Figure 8.** The vertical component of position (**a**) and velocity (**b**) of the launch cart.

**Figure 9.** Pitch angle (**a**) and angular velocity of pitch (**b**) of the launcher cart.

Powered by a linear motor, the cart moves with a constant horizontal acceleration of 14.2 m/s<sup>2</sup> during the first second, and then decelerates with the same magnitude of acceleration (Figure 10a). The vertical acceleration component oscillates around zero (Figure 10b). The vibrations are related to the cart's stabilization around the equilibrium position.

**Figure 10.** Acceleration of the cart with respect to (**a**) the horizontal axis *xf* and (**b**) the vertical axis *zf*.

The launch cart tends towards the equilibrium position (Figure 11). The gravity and levitation forces components along the longitudinal axis having opposite values. The oscillations are slowly damped, which is a characteristic feature of the magnetic suspension of the launcher.

**Figure 11.** The components acting on the launch cart: forces relative to the vertical axis of the cart (**a**) and pitching moments (**b**).

The results of the conducted simulations indicate that the equilibrium point at which the weight of the cart is balanced by the levitation force is 6.6 mm; the maximum height to which the cart rises due to the action of the levitation force is 7.16 mm; shifting the center of mass of the cart in the x axis direction results in a non-zero tilting moment, balanced by the moment of forces resulting from the difference in the value of the levitation force acting on the front and rear supports of the cart. The simulation studies confirmed that a characteristic feature of the launcher suspension is the low damping of vibrations, which was demonstrated in experimental studies [7].

An important attribute that distinguishes the levitation launcher from other unmanned aircraft launchers is the maintenance of a constant acceleration value for the launch cart during takeoff (Figure 10), ensuring a consistent acceleration throughout the entire runway. It is worth emphasizing that the value of this acceleration can be easily controlled by increasing or decreasing the thrust force of the linear motor.

#### *4.2. Lateral Motion of the Cart*

In the next simulation, the motion of the cart detached from the linear motor was considered, with its center of mass shifted by 2 mm along the lateral axis of the magnetic tracks. The simulation took into account the forces resulting from the phenomenon of flux trapping inside the superconductor.

It was assumed that the superconductor was immersed in liquid nitrogen at a height of dFC = 3 mm, similarly to the experiments described in Chapter 3. The superconductor is then subject to the resultant force *FLs* described by Equation (7), as well as a torque aiming to align the magnetic moment vector of the superconductor with the external magnetic field induction vector.

The trajectory of the cart's motion is shown in Figure 12. When detached from the motor stator, the cart moves along the rails of the launcher in a weakly damped periodic motion. The oscillations of the left and right supports are synchronized in antiphase, aiming to maintain the center of mass of the cart aligned with the center of the magnetic rails.

**Figure 12.** The 3D trajectory of the cart's movement after destabilization with respect to the lateral axis.

Figure 13 shows the course of the position and orientation components of the cart. It can be seen that the cart was set in an undamped vibrational motion relative to the center of the launcher tracks (Figure 13a). In the simulations, the oscillations of the cart deflection angle (around the *xf* axis) were observed (Figure 13b), which explains the spatial trajectory behavior (Figure 12)—at the moment when the right supports of the cart are the highest, the left supports reach the lowest position, and then the situation is reversed. The cart oscillations along the longitudinal *xf* and lateral *yf* axes are synchronized and are in antiphase. The oscillation period is 0.45 s, and the amplitude of lateral motion is 2 mm. The maximum tilt of the cart in periodic motion is 0.033 rad. As the cart moves away from the center of the rails, its height increases (Figure 13c). At zero displacement of the cart along the *yf* axis, its height is 5 mm, which is smaller than the equilibrium point determined for a cart solely subjected to the Meissner levitation effect.

**Figure 13.** Position and orientation components of the cart with respect to the lateral *yf* and vertical *zf* axes associated with the magnetic tracks: (**a**) position of the cart relative to the *yf* axis, (**b**) roll angle, (**c**) cart height, (**d**) cart height as a function of displacement along the lateral axis *yf*. The green line marks the trajectory during the first 0.18 s, when the motion of the cart is aperiodic.

Figure 13d shows the height of the cart as a function of the displacement of the cart along the lateral axis *yf*. The green line marks the trajectory during the first 0.18 s, when the motion of the cart is aperiodic. After that, the cart is set in periodic vibrations marked on the graph in red. The graphs show clear nonlinearity and hysteresis of the cart motion.

Based on simulation studies, it was found that the unpowered launch cart, when its center of mass is slightly shifted, does not return to its equilibrium position. The cart vibrations are not damped. This shows that the lateral motion of the cart is very unstable and requires additional position stabilization, e.g., using a linear motor.

#### **5. Numerical Investigations of the UAV Launch Process from a Launcher**

The complement to the analyzed dynamics of the magnetic launcher cart is the analysis of the process of the UAV uncontrolled takeoff, which was conducted taking into account the influence of disturbances resulting from air movements in the form of constant wind and gusts. The test UAV considered was a Bullit airplane (Figure 14). The aerodynamic model of the UAV [18] was examined. The main data of the analyzed UAV are gathered in Table 1.

**Figure 14.** The Micro UAV Bullit by Topmodel [19].


**Table 1.** Parameters of the Micro UAV Bullit [18,19].

According to the takeoff procedure developed within the GABRIEL project, the aircraft is mounted on the launch cart prior to takeoff. A linear motor is used to accelerate the cart. When the system reaches the separation speed, the aircraft is detached from the cart and continues its independent flight.

The numerical simulations were conducted assuming that, at the initial moment, the UAV is located on the launcher cart above the tracks, and the aircraft's propulsion is turned off. The motion of the system is initiated by the launcher cart, which is accelerated by a force of 50 N. The takeoff is performed with the UAV's propulsion turned off. In the developed numerical model, the couplings between the cart and the UAV starting from it were analyzed, taking into account the influence of constant headwind and tailwind, as well as gusts on the process of the UAV takeoff from the magnetic launcher.

#### *5.1. Takeoff of UAV in Calm Atmosphere and with Headwind or Tailwind*

In the conducted numerical simulation analyses, the results of modeling the UAV takeoff process were compared under calm atmospheric conditions, and then with the consideration of a constant headwind or tailwind blowing at a speed of 5 m/s.

The motion parameters of the launcher cart and the UAV taking off from it in calm air are presented in Figure 15 by a solid line. In this case, the UAV separation occurs at t = 2.2 s with a speed of 24 m/s after covering a distance of 28 m from the cart. The headwind shortens the UAV takeoff distance, while the tailwind lengthens it, as shown in Figure 15a. This is caused by a faster increase in the lift force in the case of headwind (Figure 15d). The separation velocity measured in the inertial system will then be lower (Figure 15c). During the climbing phase, the UAV without propulsion loses its speed.

The mutual interaction between the motions of the starting cart and UAV is shown in Figure 16. After the UAV leaves the launcher, the oscillations of the velocities of both elements of the system are visible, with the high-frequency oscillations of the velocity vector of the starting cart being damped much more slowly. This is a characteristic feature of magnetic suspension, which has been confirmed in experimental studies [6].

#### *5.2. Gust Impact on UAV Takeoff*

In the analyzed cases of uncontrolled UAV takeoff in a gust, it was assumed that the gust of wind occurs along the horizontal axis and is described by the function

$$V\_w = \begin{cases} 0 & \text{if } t < t\_s \text{ or } t > t\_s + d\\ A \sin\left(\frac{\pi}{d}(t - t\_s)\right) & \text{if } t\_s \le t \le t\_s + d \end{cases} \tag{9}$$

where it was assumed that the wind gust occurs along the horizontal axis and is described by the function, where the gust amplitude *A* = 5 m/s, duration *d* = 1 s, and the gust onset time *ts* = 2.3 s coincides with the moment of UAV detachment from the launcher.

On Figure 17, the results of the analysis of the start without wind are compared with the results of the start with variable head and tailwinds. The simulations conducted show the effect of the gust on the change in UAV trajectory and speed. The headwind causes the UAV speed relative to the ground to decrease. Conversely, in the case of tailwind, its speed relative to the ground increases (and at the same time, the aerodynamic speed decreases), which in special cases can lead to exceeding the critical angle of attack.

**Figure 15.** Influence of constant wind on (**a**) UAV trajectory;(**b**) launch cart trajectory; (**c**) UAV velocity; and (**d**) its lift force.

**Figure 16.** (**a**) Climb rate and (**b**) pitch angle of the cart and UAV during takeoff in calm atmosphere.

During the analyzed gust, the kinematic parameters of the launching cart do not change because the gust starts at the moment when the UAV departs from the launcher. However, it does affect the dynamics of the cart, as shown by the levitation force graphs (Figure 17b). During the gust, the levitation force increases or decreases its value in a manner analogous to the function of the wind speed change.

In the conducted numerical analyses, the relationship between the wind and the motion parameters of the UAV launch system was shown. The impact of wind on the parameters of the launch system, such as the trajectory, speed, tilt angle, and angular velocity, as well as the lift force of the launch cart, was examined. Compared to the trajectory without wind, headwind shortens the takeoff distance while the tailwind lengthens it. The results of the conducted numerical studies confirm the correctness of the developed mathematical model, demonstrating that the design of the analyzed launcher enables a

safe takeoff. Based on these results, the optimal parameters for the magnetic launcher system can be selected. Primarily, the dimensions of the magnetic rails must correspond to the size of the launching UAV, ensuring that the generated levitation force is capable of lifting it. Therefore, in the analyzed simulation cases, the Micro UAV Bullit was taken into consideration.

**Figure 17.** The influence of the gust on (**a**) UAV trajectory; (**b**) the levitation force of the launching cart; (**c**) the longitudinal velocity of the UAV; and (**d**) its vertical velocity.

#### **6. Conclusions**

The article describes the process of modeling and analyzing the dynamic properties of the launch cart of an innovative prototype of an unmanned aircraft launcher utilizing passive magnetic suspension with high-temperature superconductors. The aim of modeling the dynamics of the launch cart of the analyzed launcher, constructed as part of the GABRIEL project, was to specify the operating conditions of the system. Two types of cart motion were analyzed. In the first case, the motion of a cart driven by a linear motor, which constrains lateral movements, was considered. Then, the behavior of the driven cart detached from the motor stator was analyzed, allowing for the examination of the spatial behavior of the cart on the launcher rails, taking into account its lateral movements and the influence of forces resulting from magnetic flux trapping inside the superconductors. Complementing the dynamics analysis of the magnetic launcher's launch cart were numerical simulations, in which variable atmospheric parameters and couplings resulting from the mutual interaction of individual elements of the launch system were considered. The analysis included the study of the process of uncontrolled UAV launch, taking into account the influence of constant headwind, tailwind, and gusts.

On the basis of the results presented in the paper, the following conclusions can be drawn:


The analyzed launcher is an innovative alternative to the available market options for unmanned aircraft launchers. The developed mathematical and numerical model, enabling the simulation of the launch cart's motion on the launcher rail during the takeoff and landing of the UAV, indicates that the design of the analyzed launcher allows for a safe takeoff under variable weather conditions. The longer takeoff path provided by the launcher reduces the overload on the UAV during takeoff, resulting in increased safety for the launcher, aircraft, and onboard equipment. Any vibrations of the UAV system, caused by uneven mass distribution or external factors, are dampened by the magnetic forces and moments arising from the Meissner effect and the trapping of magnetic flux within the superconductors.

The most significant advantage of electromagnetic launchers, especially magnetic levitation launchers, is the reduction in acceleration and its constant value during takeoff. Conventional catapults are characterized by very high accelerations during takeoff, which can damage the aircraft and equipment. Due to the fact that, in levitation launchers, the aircraft is accelerated at a constant acceleration throughout the entire length of the runway, the takeoff process is very smooth. It should be emphasized that the value of this acceleration can be easily controlled.

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

**Funding:** Research was funded by the Warsaw University of Technology within the Excellence Initiative: Research University (IDUB) program. The APC was funded by Air Force Institute of Technology.

**Data Availability Statement:** Data supporting reported results are contained within the article.

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

#### **References**


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

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

*Electronics* Editorial Office E-mail: electronics@mdpi.com www.mdpi.com/journal/electronics

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

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-8425-6