3.1. Target Probability Density Function
3.1.1. Fire Hazard Risk
The concept of uniform coverage-based control laws obtained in
Section 2.3 is the convergence of the time average along a trajectory into the desired spatial measure. The spatial measure is then one of the inputs of the algorithm.
As previously mentioned in
Section 2.1, the only restriction applied to this measure is that it needs to be bounded, without infinite values in the space domain that needs to be surveyed.
For the purpose of this work, it is convenient to think of this measure as something that relates to a fire hazard risk. A measure that would make sense is the Fire Hazard Risk Index, a value calculated daily by Instituto Português do Mar e da Atmosfera (IPMA). According to the methodology document that is publicly available [
13], this value is a number from one to five, encompassing the different risk levels from reduced to maximum, calculated using various data such as weather prediction and the amount of fuel on the ground.
It is noted that the pixel resolution in IPMA’s calculations is only one square kilometer, and would not be enough for the purpose of monitoring small- to medium-sized areas, if needed. However, the classification method can be utilized to make custom risk maps based on the knowledge of the owner of the surveillance system.
On a more practical level, there is no information regarding the fire hazard probability or fire rate and its relation to the risk index besides the fact that a moderate index would have a higher probability than a reduced index, for example. Consequently, relative numbers are the easiest to use and a very practical choice that can be easily tweaked.
The fact that the Fourier Transform is applied to the measure and only its coefficients are used in the algorithm makes it so the relative numbers must not be so close together as to be indistinguishable in the reconstruction of the original function. The closer the relative numbers, the easier it is for the sinusoidal nature of the reconstruction to mix the different risk zones.
For the purpose of the simulation, the values for relative risk from
Table 1 are proposed and were used.
To understand the effect of these numbers in the context of this work, assigning the classifications Reduced and Moderate to two area zones of the same size would mean that, ideally, the UAV network would spend twice as much time monitoring the Moderate area as it would monitoring the Reduced area. Additionally, the results would be the same if instead of those classifications, any two adjacent classifications were used because each classification has a relative risk twice that of the last classification.
Depending on the number of classifications used in the entire surveillance area, the relative risk can be seen as a measure of time spent per area.
3.1.2. Designing the Target Probability Density Function
Considering that further in the work there will be a probability analysis, it is convenient that the measure for the Fire Hazard Risk is expressed as a probability density function with a volume of one in the entire domain.
Let
be the square target domain for surveillance and
be the target probability density function. Using an auxiliary function,
, which assigns a relative risk number to every point in the domain, we can have the target probability density function defined as:
for every
. The function is zero for every
outside of the domain
U.
The Fourier coefficients usable by the control law are recovered from (
50):
The number of harmonic wave vectors is an important factor and depends on the desired resolution and can heavily affect computation time. With the highest wave number
K, the number of integrations to obtain the Fourier coefficients from (
50) alone is
for each time step.
The resolution, on the other hand, is heavily influenced by the characteristics of the UAV used for surveillance. There is no point in having a one meter resolution if the UAV can only make coordinated turns with a 10 meter radius, for example.
For simulation purposes, the domain picked is a 2000 m × 2000 m domain with a cell resolution of 100 m × 100 m.
The risk distribution used is shown in
Figure 3.
The function used has only four levels, with the maximum being Very High.
3.1.3. Disadvantages Brought by Fourier Symmetry and How to Solve Them
While
Figure 3 shows a very reasonable probability risk map, the Fourier Transform needs to be applied to it to have any use in control law.
However, one of the characteristics of the Fourier Transform when the Fourier Basis in (
4) are used is its periodicity, or rather the symmetry on the edges. This is showcased in
Figure 4a. Additionally, both control laws heavily depend on the gradient of the basis function to decide what the next input is. The result is that if no other modification is made, the UAV would very quickly leave the intended surveillance domain and keep working “properly” according to the decrease in metric, as the metric does not differentiate between the periodic areas.
In [
18], this is not a problem for one of the models developed in the article, the single integrator model from
Section 2, a model in which the control input directly controls the linear velocity of the particle, because there is no lag in the response and the boundary conditions of the Fourier base function would make it so the closer to the boundary, the closer to zero the gradient would be. In this case, there would never be a linear velocity value that would point to the outside of the domain while the particle is on the boundary.
In the case of this work, the control input effect on the position is delayed, and even if the gradient vector gets close to zero in the boundary of the domain, there would often be cases where the UAV would not be able to turn swiftly enough and the maneuver would lead it outside of the domain. When outside of the domain, the metric and gradient would keep functioning as if nothing is wrong and the UAV would stay in that periodic zone until it eventually left again.
Figure 4b illustrates what happens if no countermeasures are employed when running the simulations detailed in
Section 3.3. In the figure, the plotted squares represent the different reflections of the 2000 × 2000 m domain to be covered.
There are two proposed solutions to this issue:
Forcing the vehicle to turn back to the initial domain, effectively overpowering the control law,
Extending the domain of the Fourier Transform with the desired distribution surrounded by zero values.
The first proposed solution, while functioning, would result in that every time the UAV is outside the domain, the trajectory would no longer be optimal in the sense that it would not be controlled by the optimal control law anymore. This would result in an undesirable increase in the metric for uniform coverage. Testing this solution alone showed that often times the vehicle would go outside the domain again right after returning, creating a loop that ensured the UAV stay in the boundary of the domain for long periods of time.
The second proposed solution requires an increase in the number of harmonics used in the Fourier Transform to keep the shape of the desired probability density function. This is not a problem with current hardware, either on-board or off-board, as the increase in harmonics increases the computation time linearly and not exponentially, as shown in
Section 3.3.1.
There is another issue regarding the second solution, in the sense that the Fourier Transform is only an approximation of the risk probability density function and not the exact function, meaning that a value of zero on the outside would never be quite zero. Depending on the parameters of the model, there could still be a small chance that the vehicle would find itself outside the domain, albeit a small one when the parameters or the model are well chosen, something that is discussed in
Section 3.4.1.
Consequently, it makes sense that a combination of both is used. The domain of the distribution is extended to 3000 m × 3000 m with the previous function centered in this domain, and should the UAV leave the extended domain, it is forced inside. Since the probability on the boundaries of the extended domain is approximately zero and increasingly approaching the actual coverage domain before the extension, it is very unlikely that the UAV would be stuck in the extended domain border in a loop.
Figure 5 is the graphical representation of the extension of the domain.
In
Figure 6, the Fire Hazard Probability and its Fourier Transform with
harmonics are represented. The small rectangle with the minimum cell resolution was inserted to make an example of how small changes influence the Fourier Transform. It is obvious that without a considerable amount of harmonics, the 100-m width rectangle would not be very influential in the final result.
With this in mind, it would be useful in a real life example to ignore very small zones of smaller risk near relatively bigger high risk zones and just assign a similarly high risk to the smaller zone. In contrast, a very small zone of high risk is still very important and should be extended in the risk function to account for the Fourier Transform, otherwise the transform might decrease the risk the algorithm takes into account.
It is, in conclusion, imperative to understand that the Fourier approximation of the Fire Hazard probability density function is the important factor in the algorithm, not the pre-transform distribution.
3.2. Simulink Simulations
The Simulink simulations will only take into account the models presented in
Section 2.2.1 and
Section 2.2.2. While they do not represent the full workings of an aircraft, the simulations will not take into account any inner-loop controllers as there is no set aircraft to be used in the project.
The design of the inner-loop controller is outside of the scope of this thesis, but the algorithm can output either heading angles, heading rate, linear velocities, or position waypoints. With this variety of variables, many off-the-shelf solutions can be used, such as a PX4 controller [
26].
The block diagram in
Figure 7 shows the basic version of the full controller when taking into account the aircraft inner-loop controllers. The section Simulink simulations will focus on the highlighted area. Every obtained result uses the ode4 Simulink solver [
27] (fourth-order Runge–Kutta method) with a fixed step size of 0.1 s.
The Simulink model is a simple integrator block preceded by an Interpreted MATLab function block containing all the calculations necessary to obtain the control law.
The Fourier coefficients of the desired distribution as well as the inner product between the Fourier basis function and itself were calculated using two nested one dimensional integral functions.
Regarding basic vehicle characteristics for parameters such as trim speed, heading turn rate, and flight time, in-depth market research, compiled in
Table 2, showed that commercial UAVs for surveillance have an endurance of hours. One specific case is the hybrid VTOL/Fixed-Wing aircraft Jabalis from Helvetis, with an endurance of up to eight hours flying at 39 m/s. Turn rates, considering the type of aircraft, are expected to go up to around 0.5–0.8 rad/s.
3.3. Dubin’s Vehicle Model
For the Simulink simulation for the Dubin’s Path Model, we considered an ASN composed of 3 UAVs with a maximum heading rate of rad/s and a constant vehicle speed of Va = 30 m/s based on the UAV market research. The simulation was run for a period of T = 3600 s, i.e., one hour.
The obtained trajectory is highlighted in
Figure 8, with the different colored lines depicting different UAVs.
It is evident from
Figure 8 that the solutions for the symmetry problem were effective. In fact, the fail-safe measure in the case that the UAV was found outside of the bigger domain was never executed. The UAVs still spent some time on the outside of the actual surveillance domain, but only in the beginning of the simulation due to the nature of the control law prioritizing the lower frequency harmonics first, due to their weighting decreasing as the harmonic’s number increases (
decreases as the index
k increases, from Equation (
8)).
The density of the lines is also higher in higher risk zones, and much lower in lower risk zones. The two corners in the diagonal are evidence of this, as the lower corner has the lowest Fire Hazard risk and the higher corner has the highest Fire Hazard risk.
It is also important to note that while the control input can only be one of two values, the vehicles do not move in perfect circles every time. This makes sense as the control law is calculated every 0.1 s. For example, if in the period of 2 s the inputs were always the opposite of the previous one, the end result would have been something very similar to a straight line.
In
Figure 9, the control input sequence for all UAVs is presented in a small excerpt of the simulation. The control input is, as can be expected, very erratic, and changing values very quickly, from the positive to negative end of the spectrum. Measures to try and overcome this issue were studied in
Section 3.3.4Some controllers can control a UAV solely given a reference turn rate, as it is very easy to translate to a roll angle reference in a fixed-wing UAV. Regardless of the controller, any reference input would be hard to follow due to the discontinuities of the reference. This is not the only drawback, of course. Given the way the input is defined, the servos in any aircraft would need to be constantly changing their positions from one extreme to the other, leading to a reduced vehicle lifespan, higher maintenance, and lower endurance due to constant use of the servo motors.
In regards to the metric for uniform coverage in
Figure 10, the first thing that needs to be said is that the graph is normalized against the metric in the starting point,
. This makes it easier to compare the metric data against other models and variables, as the starting metric is decided by the size of the domain, the risk function, and the initial vehicle positions.
In terms of the shape of the graph, it is clear the function is not monotonous. That would be impossible using this uniform coverage algorithm. On average, the metric appears to decrease exponentially in
Figure 10a. This result is in line with the results obtained in [
18], so it can be concluded that there was not any anomaly in the application of the algorithm when changing the model.
The noise has clear peaks which immediately decrease, almost in a recurring fashion. This is the result of the time average being used. This is better explained when thinking of a uniform distribution in the entire domain. Even if the algorithm was perfect, and the distance between the functions was 0, the next movement of the UAV would increase the distance between the distributions, and it would only reach 0 again once the rest of the domain was visited.
While the metric for uniform coverage is important, the shape of the obtained distribution is even more so. In
Figure 11, the desired distribution and the obtained distribution are compared. The values do not perfectly align with the values from
Figure 5, the distribution pre-Fourier Transform, but overall the different risk zones can be differentiated in the image.
3.3.1. Effect of the Number of Harmonics in the Computational Time
The process of running the simulations can be divided into two different stages: The setup stage, in which the Fourier coefficients of the desired distribution are obtained, as well as defining the initial conditions, and the simulation stage, where the Simulink model is run.
The number of harmonics used is directly related to computation time, both in setup and in simulation. For the setup stage, each harmonic adds two 2D integrals, one for the denominator and one for the numerator in (
50). For the simulation, each harmonic effectively adds one state, which means one additional integral.
In
Figure 12, the results for the computational time using different values of
K are shown. The increase in computational time in the setup stage is linearly related to
K, which is surprising because the number of harmonics is equal to
. This means that increasing the maximum number of harmonics is not such a detriment to computational time. The same can be said to the time it takes for the simulation to run, even if there is a slight inconsistency on the lower number of harmonics. This can be attributed to the way Simulink calculates the integrals, which depends heavily on threshold values and numerical precision.
3.3.2. Effect of the Constant Speed and Heading Turn Rate
In order to see how the basic parameters affected the metric, the simulation was run for different values of each of the parameters. Values from
m/s to
m/s were tested, as well as values from
rad/s to
rad/s. The results are shown in
Figure 13.
The obtained results were very similar. There were no drastic differences in the metric for uniform coverage’s behavior, only that with each increment added to each variable, the averages seem to decrease slightly. This is more evident in the case of the linear speed.
In conclusion, the higher the linear speed and/or heading turn rate, the better the results. This is to be expected, of course, as with a higher linear speed, more of the domain can be covered in a shorter time span.
In the case of the heading turn rate, it can be concluded that a higher rate is better, as the minimization of the Hamiltonian was given for the bounded control input and the control law is directly proportional to the boundary values. The higher the absolute value of the maximum heading turn rate, the more negative the derivative of the metric for uniform coverage is allowed to be, resulting in a better performance.
Ultimately, the values for the linear speed and heading turn rate will be determined by the maximum the selected vehicle is able to provide.
3.3.3. Effect of the Number of UAVs
The difference in the number of UAVs used in the simulation is very impactful.
In
Figure 14 the results for varying the number of aircrafts between 1 and 5 in terms of the metric for uniform coverage are shown.
3.3.4. Regularization of the Control Input
Invariably, increasing the number of UAVs increases the decay rate of the metric, although it is clear that there are two different groups of results, and . The clear difference between the results can be attributed to the area of the domain in combination with the speed at which the UAVs are allowed to move. The results clearly show that the difference in the metric between 2 and 3 UAVs is considerably bigger than the difference between 3 and 4 UAVs. It can be concluded that for these particular combinations of parameters, the ideal size of the fleet would be 3 UAVs, as the performance would decrease heavily with less than 3 and stay relatively the same after 3 UAVs.
The difference between the results in the second group is the height of the cyclical peaks. The higher the number of UAVs, the lesser the height of the peak. This is to be expected, as the metric is averaged over the number of aircrafts in the simulation and the cyclical increase in the metric is due to the problem constraints. The increase in UAVs also means that the frequency of the up and down cyclical behavior is higher, as more aircrafts would lead to a decrease in the time it would take to reach the minimum value in the metric for uniform coverage due to triple the coverage, even if it is average.
In order to try and solve the issue of the irregular input, the control law from (
36) was modified. Instead of using the
function, a saturation function was used, so that the regularized control law is defined as follows:
where
is a constant to tune the slope of the linear part of the saturation function.
The constant
is then a parameter to be changed, however given the way the quantity
is defined in (
35), it is nearly impossible to predict even its magnitude order. Furthermore, if the target distribution
changes, or the domain changes its size, the order of magnitude of this quantity will also change because
is directly related to the Fourier coefficients of
. The ideal value for
would be one that would allow for the control input to be almost always in the linear zone but have its values as close to the saturation limits as possible; meaning the higher the slope, the closer to optimal the regularized control law is.
While the order of magnitude is hard to predict, the equation for
, (
35), resembles the equation for the metric for uniform coverage, (
7), so as a first approximation,
was set to the inverse of the order of magnitude of the metric at
. In the particular case of this simulation, this value is
.
The trajectory of the different agents and their regularized control input throughout the simulation is showcased in
Figure 15. The most noticeable effect of this specific regularization is that the agents spend considerable time outside of the domain to be covered, even going so far as to frequently leave the extended domain. This is because the
parameter is too low, causing the system to not have a control input high enough to make the needed turn. With regularization at a low slope, the control input is much closer to zero than it is to the saturation point, which would be the optimal control input. This can be seen in
Figure 15b, where the control input is now always smooth and rarely ever goes above
, when the saturation point is
. The values in the saturation zone are not due to the regularization but to the method employed to keep the agents inside the extended domain.
With poor results for the trajectory, it is expected that the obtained metric for uniform coverage is also inferior when compared to the optimal result in
Figure 10. This is evident from
Figure 16, with the normalized metric never dropping below
when the optimal result almost reaches the
values.
In order to try to improve this, the slope was increased to because there was still the option to do so, with the control input being so far from the saturated values.
The trajectory result for this value of
is showcased in
Figure 17a, and it seems that the agents spend considerably less time outside of the domain to be covered, even if it is longer than if the input was optimal, and can still leave the extended domain. The regularized control input for the simulation is showcased in
Figure 17b, and it can be seen that it reaches the saturation values very often, meaning that increasing the value of
would only make it so the control input saturates more often, which would defeat the point of adding saturation.
Regarding the metric for uniform coverage, shown in
Figure 18, the values dropped below the
mark, but were not quite as low as the ones obtained from running the simulation with the optimal control law, as shown in
Figure 10. Therefore, the proposed regularization fixes chattering of the input signal but exhibits some degradation in the performance. It would be counter productive to increase the value of
any more, as it would saturate more often and defeat the purpose of its implementation in the first place.
Furthermore, the metric is very sensitive to changes in but there is no formula that can generate the ideal value. The values used were obtained after analyzing the initial metric, but the metric’s order of magnitude depends on the number of harmonics used, the risk function, and the domain size. To obtain the ideal value, extensive experimentation is required on a case-by-case basis.
3.4. Adapted Dubin’s Vehicle Model
When the Adapted Dubin’s Vehicle Model is introduced in
Section 2.2.2, the conclusion that the
y coordinate of the point in the body frame (
) does not affect the determinant of the final matrix is very important, as it can be set to zero and the system will still work perfectly fine. As such, with this
value, the point is always situated in the
axis of the body frame.
For this simulation, the base speed of the vehicle was set to m/s and its maximum heading turn rate was set to rad/s, which are the same values used in the previous simulation.
For the values specific to this model, the distance
was set to 2 m, and the speed differential was set to
m/s. The effects of varying the distance of the point will be studied further into this work, in
Section 3.4.1.
Choosing the value for needs to be done carefully because it directly relates to the minimum speed of the UAV. In this case, the minimum speed would be m/s. The value needs to be above the aircraft stall speed. The time span of the simulation was, again, one hour, s.
Simulating three UAVs using this model produced the trajectories in
Figure 19, with the different colored lines depicting the different UAVs.
The trajectories still have an acceptable shape, without any sharp edges, unobtainable in the model. With the chosen parameters, the UAV stayed inside the extended domain and was never forced inside the domain with a non-optimal input.
The density of the lines also clearly shows the difference in surveillance time between the different areas on the domain, with the top right being the densest zone, corresponding to the area with the highest Fire Hazard Risk.
A small sample of the control input for one UAV over the course of 10 s is shown in
Figure 20. The control input is smooth, even if the plot chosen was a
stairs function instead of the regular
plot to account for the discrete solver used.
Given that represents the heading turn rate, it is easy to see why the input would result in trajectories that are very easy to follow by a fixed-wing UAV. These values are small for the most part, with the transition being slow. It does not have the steep transitions from the previous model and so there should not be any agility imposition in the UAV, meaning a wider variety of aircrafts can be picked for usage.
The evolution of the metric for uniform coverage in this simulation, shown in
Figure 21, exhibits an exponential behavior in logarithmic scale, confirming that the system is working as expected. The result is, once again, not a monotonous function but decreases on average as the simulation time increases.
The metric once again shows cyclical behavior after an initial transitory period, likely due to having covered the higher weighted Fourier bases, the lower frequency ones.
The heatmap of the probability density function obtained at the end of one hour of simulation is shown in
Figure 22. The heatmap shows that the
Very High and
High areas were appropriately covered, even if the edges of the
High Risk zone do not exactly have the desired shape, at around
.
The Moderate areas, on the other hand, are considerably messier, with the obtained distribution having plenty of darker blue spots, associated with a Reduced classification. The edges can clearly be noticed, however the content is not entirely uniform. This is a case where a slightly higher simulation time would be very beneficial, as the finer details only get worked on after the first basis converges to acceptable results due to the weighting.
The shape of the surveillance square domain is also very blurry, due to the vehicle spending a considerable amount of the beginning of the simulation doing a short sweep on the extended domain. This is due to the parameters of the system, and will be discussed next.
Overall, the obtained result is still well within the acceptable domain. Even if the shapes are not perfectly defined and uniform in detail, it is clear the higher risk areas are being surveilled appropriately and the rest would, with a higher running time, converge to desirable uniformity.
3.4.1. Effect of Parameter Variation
The system has one defining characteristic that makes it so multiple inputs are allowed, and that is the tracking of a point on the body frame and not the center of the body frame itself.
As a result, the position of the point directly influences the amount of control over the UAV. While the control of the linear speed differential is important, it is inevitably the behavior of the heading turn rate input that makes the trajectory.
The control input component related to the heading turn rate is the second column of matrix , which is directly proportional to , the x coordinate of the point in the body frame. As a result, the closer to the origin, the less effect the control law equation will have on the heading turn rate component due to the input vector being unitary. If the vector norm of the columns of is very different, one component will have a bigger weighting.
In conclusion, it is better if and are close, or at the very least in the same order of magnitude.
To exemplify the lack of control, the same simulation was run, however for
m instead of the previous 2 m. The resulting trajectory is shown in
Figure 23, and it is clear that the behavior is completely unacceptable.
Before any further analysis,
Figure 23b shows the UAV going outside the extended domain, applying the second fail-safe measure of forcing the aircraft back inside. This figure proves the designed algorithm is working properly.
Both simulations use the same parameters, however have a different distance, with rad/s and m/s. To expand on the previous concept of similarity in values, the product between the heading turn rate and the distance is, in the case of simulation (a) and for simulation (b). Comparing both of these values to the value of makes it clear why the parameter is not at all well chosen. In this case, for there to be a significant change in the control input related to the turning rate in (b), the value in the multiplying vector from the control law would need to be much higher than the value related to the velocity due to the normalization of the vector.
The result of the bad choice of parameter is clear from
Figure 23b, in which the UAV has difficulty turning and is less agile. The effective consequence of bad parameterization is almost an attenuation in the turn rate, which explains the wider curves produced by the model.
The direct effect of a bad trajectory was translated to the metric for uniform coverage, shown in
Figure 24. While the graphs look similar at first glance, the order of magnitude on the
axis is of note, as the value in (a) is almost an order of magnitude lower than that of (b).
This figure shows that, while the control law still works for bad parameterized systems, the metric for uniform coverage would stabilize at a higher value. All the results shown up to this point suggest there is a point after which the metric for uniform coverage starts decaying at a much slower rate. The results from
Figure 24 suggest that this point is influenced by the freedom the UAV is afforded on the model, something that was never previously showcased in this work. In this case, while the metric is still low, the obtained trajectory is not acceptable. In conclusion, a decrease in the metric for uniform coverage, while necessary, is not sufficient to ascertain if the model is well parameterized or not.
While a bigger distance can mean better results, it can not be so distant that it does not make sense to calculate the distribution on an unrealistic point nor should it overpower the velocity component, as the input would quickly depreciate into a graph similar to the Dubin’s Path as would be reaching zero after normalizing.
3.5. Model Performance Comparison
Since both the previous models had the same initial condition, same domain, and same risk function, it is possible to overlap the metric for uniform coverage graphs.
In
Figure 25, the two metrics are overlapped and look very similar. The Dubin’s Vehicle model has a slightly lower average value, but the results could have easily been passed as each other’s, as they have no noticeable differences. With no noticeable differences in performance, it is up to a matter of preference as to which model to use.
The Dubin’s Path model is very simple to parameterize. Every parameter that the model has, the adapted model has as well. On the other hand, the control input is not smooth and can lead to issues with the aircraft and possibly bring difficulty in the control of the aircraft.
The Adapted Dubin’s Path model has the flexibility of having more parameters and it is not out of the realm of possibility that there is some combination of values that would perform significantly better than the Dubin’s Path mode, although none were found during the development of this work.
The biggest difference between the two models is then the fact that the control input in the Adapted Dubin’s Model is smooth. This is exactly the issue that needed to be solved and the reason this model was applied, and the fact that there is no drop in performance is a positive aspect of the Adapted Dubin’s Vehicle model.