2.1

System dynamics

For the simulation, a model derived from the Newton and Euler equations presented in [20] was used. Assuming the vector [x y z ϕ θ ψ]^T as the position and orientation vector in the F^e system, as well as the vector [u v w p q r]^T for the vector of linear and angular velocities in the F^b system, the dynamic UAV model takes the following form:

with the following definitions: s(α) := sin(α), c(α) := cos(α), t(α) := tan(α), and m as the UAV’s mass. Inertia matrix I, due to platform construction symmetry, takes the diagonal form of

Minor damping reflecting aerodynamic drag has been added to the model [22]. Diagonal matrices of displacement and rotation aerodynamic coefficients are expressed by the following relationships:

By organizing individual variables, a state vector is obtained:

along with the control input vector

This vector contains, respectively, the value of the resultant lift generated by the engines set f_t, torques τ_x and τ_y in the Ox and Oy axes of the UAV platform, and torque τ_z in the Oz UAV axis. Finally, the form of a nonlinear dynamic system is expressed as

2.2

Linearization of nonlinear model

The four-rotor model is a highly non-linear system with many coupled variables. In order to enable linear analysis of the system, the model should be linearized, assuming an appropriate operation point. The most frequently chosen operation point is the equilibrium condition obtained when the four-rotor model hovers. This means no UAV movement, so its dynamics are described by the following differential equation:

where x* = [0^1×9 x y z] and u* = [mg 0 0 0] are nominal values for the given operation point. This condition occurs when all four engines generate a net thrust equal to the weight mg of the UAV in the opposite direction. Under ideal conditions, there are no other forces or moments present.

The system was linearized using the Control System Toolbox of MATLAB/Simulink software, thanks to which the LTI (Linear Time Invariant) state space system described by two vector-matrix equations was obtained:

where A∈ℜ^12×12 is the system matrix, B∈ℜ^4×12 is the control matrix, C∈ℜ^12×12 is the output matrix, and D∈ℜ^12×4 is the transmission matrix. For the calculations, the parameters contained in Table 1 were used and the following results were obtained:

Table 1.

Parameters used to create UAV model.

The parameters used for UAV modeling are presented in Table 1.

3.1

Model discretization

This application uses a discrete controller with a sampling time of T_c = 0.1s using the PSO algorithm for cost function minimalization. The task of the algorithm, considering the assumed quality indexes, is to select such a value of control vector u(k) at a given moment k that the difference between the predicted state x(k+1|k) and corresponding reference point x_ref(k+1) is minimal. Due to its discrete nature, the dynamic model (9-13) has also been transformed to the discrete form:

with c2d() command in MATLAB. For sampling time Ts = 0.05s, the following results were obtained:

3.2

PSO algorithm for control optimization

The PSO algorithm searches for solutions in the n-dimensional space of the function to be optimized [18]. This is accomplished by iteratively moving the population of p particles to find coordinates for which the value of a given cost function is minimal. Each i-th particle has the following known parameters:

Figure 3.

PSO algorithm block diagram.

The particle with best global value g_best ∈ ℜⁿ is also determined for the entire population. The evolution of the algorithm is described by expressions:

where w is the inertia weight, c₁ is the value of the so-called personal coefficient, c₂ is the value of the social coefficient, and rand() is a pseudo-random function which selects a number from a range [0, 1]. The movement of each particle depends on its current velocity, its best own position so far, and the global best position for time (iteration) t. The c₁ and c₂ coefficients determine the tendency to choose the direction of particle movement. The rand() function gives the algorithm a pseudo-random character, while increasing its searching performance.

Due to its stochastic properties, the PSO is a metaheuristic algorithm [18, 23]. It does not guarantee an optimal solution, but only determines a sub-optimal variant that meets given assumptions, as close as possible to the optimal solution. This allows for the reduction of the computational load of the solution, at the expense of its accuracy. This is especially useful when reaction speed is more important than precision. The algorithm ends its execution after reaching the iteration limit or meeting the given criterion, returning the best solution so far g_best.

3.3

Search space and cost function optimization.

During the algorithm execution, particles explore the U⊂ℜ⁴ space of the control vector u values. At the beginning of each iteration (after the movement of the entire population), a cost function is calculated for each particle, considering its current position. In the proposed application, the function being optimized takes a quadratic form:

where J^t_i is the cost function value for the i-th particle in the iteration t, u_i(k+1) means the coordinates of the i-th particle (selected values for control signals), x_i means the state that the UAV would achieve by applying the control determined by the coordinates of the particle i, and u_nom is the nominal value of the control vector (for the equilibrium point). Disturbances and measurement errors are not considered [21]. For known initial values of the state vector, the values of x_i(k+1) are determined directly from the expression (14).

Diagonal weighting matrices Q ≥ ℜ^12×12 and R, R_Δu ≥ 0∈ℜ^4×4 determine the impact of each factor in total value of the cost function. The Q matrix is responsible for the influence of error magnitude between state variables and their reference values (fitness criterion). The R matrix determines the reaction to the difference between the nominal values of the control signals and the values selected in a given solution (energy criterion). The R_Δn matrix is responsible for the impact of the difference between the selected control and the current signal level (signal dynamics criterion).

Figure 4.

a-d) PSO solution searching process.

The algorithm allows for the explicit limitation of the values of the control signals. By imposing a condition on the allowed positions of particles in the search space

it can be ensured that no control value outside the specified range is selected.

4.1

Initial parameter configuration

The following parameters for the PSO algorithm were adopted for the simulation tests: number of particles p = 40, iterations limit it_lim = 150, inertia coefficient w = 0.9298, and personal and social coefficients c₁ = c₂ = 1.4986. The diagonals of the weight matrices were selected as follows:

The search space of the PSO algorithm has been limited by the following ranges:

where this range refers to the nominal control u_nom in which the value of the first coordinate mg has been compensated, so u_nom = [0 0 0 0].

In order to perform the simulations, a test trajectory was generated containing references for position, velocity and course (Figs. 5-6).

Figure 5.

Reference values for state variables of UAV.

Other state variables, without reference values, were compared with their equilibrium values x*, so

Figure 6.

3D path profile (T): S – starting point, D – destination point.

4.2

PSO algorithm improvements

The first attempts showed unstable operation of the algorithm. Control signal waveforms generated by the algorithm were characterized by strong distortion (Fig. 7).

Figure 7.

Control signal waveforms generated in initial configuration.

Despite the occurrence of undesirable distortion of control signals, the algorithm showed the ability to track the given trajectories of individual state variables (Figs. 8-9). Therefore, an improvement in the process of determining signal values could be achieved by the additional modification and tuning of algorithm parameters.

Figure 8.

State variable waveforms [u, v, w] for initial configuration.

Figure 9.

State variables waveform [x, y, z, ψ] for initial configuration.

After observing the execution of the PSO algorithm, it was noticed that the problem with stability is due to the difficulty in achieving convergence. By design, each particle should gradually approach the global best solution, while reducing its velocity in subsequent iterations. Thus, it is possible to exploit more closely neighborhood of the best solution. This effect was not achieved in the basic algorithm configuration (Fig. 10). The velocities of the particles did not decrease sufficiently, resulting in long-distance motion. The end result was a large spread around the globally best solution. Let S ⊂ ℜ⁺ be a set of Euclidean distances s_i of all particles from the position g_best, so

The average distance between the particles and g_best in the last iteration was assumed as the measure of population spread:

In order to improve the convergence of the algorithm, expression (19) was modified by implementing a factor λ that scales the velocity of particles in each iteration:

Table 2 presents the simulation results for the factor λ in the range <0.3, 0.8>. Without considering λ, the algorithm obtained the following results: average value of the dispersion index ε_avr = 1.3944, minimum value ε_min = 0.9380, maximum value ε_max = 2.0023, and total path cost J_c1 = 8.5574e + 05.

Figure 10.

Values of ε for initial configuration.

Table 2.

Comparison of experiment results for different values of λ.

Based on the simulation results, it can be concluded that the best variant is λ equal to 0.65. Fig. 11 shows a graph of changes in the index ε for the selected value λ. The difference can be clearly seen compared to the initial version of the algorithm – the magnitude order of the average value ε_avr is several times smaller, which coincides with the original assumptions.

The total cost J_c also improved, mainly due to the decrease in the cost of control signal dynamics, whose waveforms are shown in Fig. 12.

Figure 11.

Values of ε for λ = 0.65.

Figure 12.

Control signals generated for λ = 0.65.

After the stabilization of the control signals, the trajectories were smoothened (Figs. 13-14) but their shape and fitness changed only slightly. This indicates that the algorithm’s instability has no effect on its efficiency (in the context of matching to reference trajectories). The sum of squares of differences between real values and reference state variables [u, v, w, x, y, z, ψ] at all discrete moments k was taken as fitness criterion μ_e:

Table 3 compares the individual components of the indicator μ_e for λ_opt = 0.65 and λ₀ = 1.

Table 3.

Comparison of error values for the basic and modified variants of the PSO algorithm

Figure 13.

State variable waveforms [u, v, w] for λ = 0.65.

Figure 14.

State variable waveforms [x, y, z, ψ] for λ = 0.65.

4.3

Selection of weight matrices elements

Due to the lack of significant improvement in the adjustment of state variable trajectories to their reference values (despite the stabilization of the PSO algorithm), attempts have been made to improve the fitness indicator by changing the values of the weight matrices. Analyzing the results from Table 3, it can be seen that the largest errors are coupled with variables u, v, x, y, ψ with error e_ψ an order of magnitude greater than other errors. However, weighting factors for individual matrices were intuitively modified, considering that the weighting factor for e_ψ should be greater than the weighting factors corresponding to the remaining state variables being tracked. Finally, the weighting matrices took the following values:

The simulation results for the new values of the weight coefficients are presented in Figs. 15-17. Table 4 presents a comparison of individual components of the index μ_e for the modified values of R_opt and the unmodified R₀.

Figure 15.

Control signal waveforms for modified weight matrices.

Figure 16.

State variable waveforms [u, v, w] for modified weight matrices.

Figure 17.

State variable waveforms [x, y, z, ψ] for modified weight matrices.

Table 4.

Comparison of error values for the basic and modified variants of the PSO algorithm.