Energy Efficient Path Planning for the Robird

(1)

Energy Efficient Path Planning for the Robird

L. (Luuk) Lomillos Rozeboom

BSc Report

Committee:

Dr.ir. G.A. Folkertsma Prof.dr.ir. S. Stramigioli W. Straatman, MSc Dr.ir. H. Wormeester

July 2018 022RAM2018 Robotics and Mechatronics

EE-Math-CS University of Twente

P.O. Box 217 7500 AE Enschede The Netherlands

(2)

(3)

Abstract

The Robird is an ornithopter designed for bird control. The system has a large energy consumption during operation. By reducing the energy consumption, a larger flight time can be achieved. The aim of this project is to design a program that is able to find the most energy efficient path from point A to point B. For that, a simplified dynamic model of the system is developed and optimized to describe the behavior of the aircraft. It is then used in a program that iterates through different combinations of the flight conditions to find those that produce the path that consumes the least amount of energy.

(4)

(5)

Introduction

Birds can fly freely through the sky without encountering many obstacles. They can fly any- where, cover large distances in short amounts of time, are able to easily reach places that other animals would have a hard time getting to, etc. The search for this freedom is what lead to the development of airplanes but is this same freedom what causes humans to have a hard time keeping birds out of some areas.

There are several reasons why people would like keep birds out of an area. One reason is the damage cause by birds in crops. Birds can fly into farms, eating or damaging the crops and causing a loss of money to the farmer. Another reason would be a sanitary concern, getting birds away from sensitive areas to prevent the spread of diseases. There are many other examples but, at the moment, the most problematic areas to have birds in are airports. Birds flying in airports can cause a bird strike, which is when birds crash with an airplane. A bird strike can damage the airplane and put the lives of passengers at risk. On top of that, the caused damage can be costly to repair. The crash can also cause delays that add up to the expenses caused by the bird strike.

There are many bird control techniques being used at the moment. Some of these techniques are: auditory techniques such as ultrasound and explosions, visual techniques such as scare- crow, chemical techniques, habitat modification and lethal techniques. However, most of these techniques have some inconveniences like ethical issues and environmental impacts, which are mostly present in the last three mentioned techniques and, most of the techniques, lose their effectiveness over time, when habituation occurs. One of the most effective techniques is falconry. This technique consists in the use of prey birds to scare the birds away. Since it exploits the natural instincts of birds, it does not lose its effectiveness over time. The main problems with this method are that it is expensive and the fact that, because animals are used, it can be somewhat unpredictable. Clear Flight Solutions is trying to develop bird control a product that is environmentally friendly, effective and cheap that is able to compete with other bird control methods, developing the Robird, a robotic bird of prey. [1]

1.1 Robird

The Robird is an ornithopter, an aircraft that uses flapping wings to produce lift, that mimics the looks and behavior of a bird of prey. There are two models, one imitates a peregrine falcon and the other on an eagle. In this paper, only the peregrine falcon model will be considered.

The ornithopter was first designed by a Dutch falconry and airplane model hobbyist by trial and error. This means that even though a working model was developed, there is not a dynamic model describing the behavior of the system. Since the creation of the Robird, it has been analyzed and modify to improve the characteristics of the product, leading to a modified version

(8)

Figure 1.1: Photo of the Robird.

that has a better performance than the initial model and will be analyzed in this paper. [2]

At the moment, one of the problems with the Robird is that it is very difficult to maneuver.

Because of this, only highly trained pilots can fly it. To make the robot easy to fly so that everyone can use it, an autopilot is being developed. A take off controller has been developed and the landing controller is in a test phase. For the mid-flight controller, a basic program, able to go from point A to B and rotate around one point, has been developed. However, a more complex controller needs to be designed to be able to fly autonomously.

1.2 Aims and approach

The goal of the assignment is to make a simple model that is able to accurately predict the behavior of the Robird and use it to find a path that is energy efficient based on the input of the user and sensors.

To make the assignment, first a literature research must be done to find the needed theory and the work that has already been done on this topic. The research is focused in the physics behind the flight of birds, airplanes and ornithopters. Besides that, research on energy efficient path planners for aircrafts is also done.

With the gained knowledge, the model is developed. To make the model, data from test flights is used to determine the variables and parameters involved in the behavior of the system. The model is restricted to a 2D scenario. To describe the dynamics of the Robird, the model is divided into two separate ones, one that determines the behavior during gliding and another one that describes the powered flight.

With the model designed, a program that finds the most energy efficient path is developed. The

(9)

CHAPTER 1. INTRODUCTION Natural Path Planning for the Robird model takes the initial conditions of the flight and an initial and final position and calculates the path that consumes the least amount of energy.

(10)

(11)

Chapter 2

Theory

The theory behind the flight of birds is quite complex. They relay on flapping wings to generate both lift and thrust. The downstrokes are usually refered to as power strokes since they provide most of the forces involved in the flight and the upstroke is called the recovery stroke.

Regardless of the complexity of the problem, two basic laws play a main role: Newton’s third law, the wings create a donwards momentum on the air with which lift and thrust is achieved by reaction and Newton’s second law which states that the lift and thrust have to counterbalance both the weight and drag.

For lift and forward movement to be achieved, birds rely on multiple mechanisms to maximize the force exerted by the wings and minimize drag like changing the area and shape of the wings, opening the feathers and changing the angle of attack. Some of these mechanisms are quite complex and difficult to accurately reproduce. In the Robird, it is much more simple. The wing is made out of foam with a laminated resin skin, making it light and flexible. Since it cannot change the area, shape or permeability of the wing, it relays only on the change of angle of attack to generate lift and thrust. [2, 3]

2.1 Reference frames

Before starting any derivation, it is important to clearly state the frames that will be used.

Since Newton-Euler will be used to derive the equations of motion, an inertial frame has to be defined. As an inertial frame, the Earth is used and its axes are denoted by the terms X⁰, Y⁰ and Z⁰. This inertial frame is necessary for the derivation of the equations of motions however, it is convenient to introduce another frame that rotates with the aircraft and with the center in its center of mass, denoted by the terms X, Y and Z. Both frames are shown in figure 2.1 Keeping track of the relative displacement and rotation of the axes is important. For measuring the displacement, a vector is used that describes the center of the body-fixed frame with respect to the inertial frame using the coordinates Xi, Yi and Zi. For the velocities in the inertial frame the terms U_i, V_i and W_i will be used and the U_b, V_b and W_b will be used for the velocities in the body fixed frame. For describing the rotation, Euler angles are used, which involve the introduction of different reference frames with successive rotations. Since the order in which the rotations are made will have an impact on the final result, it is important to pick an order and be consistent with it. In this paper, a widely used order known as the Z, Y, X right-hand rotation sequence will be used. It consist in rotating first around the Z-axis (Ψ), then around the Y-axis (Θ) and, finally, around the X-axis (Φ), as shown in figure 2.2

The angular velocities of the body-fixed axes are important in the Newton-Euler equations. P ,

(12)

Q and R will be used for the rotations around the X, Y and Z-axis respectively with the positive direction defined by the right-hand rule as shown in figure 2.2. It is important to note that this angular velocities are not of the Euler Angles and a relation between them will be introduced later.

Some important terms are the pitch angle (θ), the the yaw angle (Ψ) and the roll angle(Φ).

These terms coincide with the Euler angles so a redefinition is unnecessary (See figure 2.2).

2.2 Dynamic equations

Despite the relative simplicity of the flight of a robotic bird with respect to a real bird, the aerodynamic phenomena involved in the movement of the ornithopter is quite complex containing many non-linear effects. Since in this project a model of a real ornithopter has to be made, it is not necessary to make a model that solves the Navier-Stokes equation to accurately describe the aerodynamic forces involved in the flight. Instead, a model that describes the moment and force balances of the ornithopter can be developed and, using experiments, the model can be adjusted to describe the ornithopter that is being modeled.[5][6]

When deciding which approach to take to get the Newton-Euler equations, two main options were considered: use multi-body dynamics or use the equations for a rigid body. The dilemma was that the multi-body dynamic model would describe the flapping movement more accurately but the model would be more complex. The equations for a rigid body instead, would not be able to model the oscillating forces due to flapping motion but the model would be simpler.

Since the overall trajectory is of interest, it was decided that the equations for a rigid body would be used with the assumption that every stroke would have a total moment and only that would have to be considered to describe the trajectory.[7]

2.2.1 Assumptions

• The Earth is an inertial frame. In reality, the earth is moving so it is not an inertial frame.

However, since the Earth is really big compared with the Robird and they move together, it can be assumed that the earth is an inertial frame.

• The equations for a rigid body hold. It is assumed that the moments and forces generated due to the flapping wings can be substitute by a total moment or force. This is considered valid because only the overall movement has to be considered so, the oscillations due to the flapping movement can be neglected and the forces and moments during a period can be substituted by a total force and moment.

2.2.2 Newtown-Euler equations

The dynamics of the ornithopter will be described using the well known Newton-Euler equations for a rigid body shown below. This equations are widely used in the description of aircraft dynamics. There are six equations, one for each degree of freedom of the problem. The first three equations describe the force balance in the X, Y and Z axis respectively and the three last ones describe the moment balance around the different X, Y and Z axis respectively. It is important to note that, since the body-fixed frame in which the equations are described is not an inertial frame, some additional terms appear due to the relative movement of the body-fixed frame with respect to the inertial frame. [4]

m ˙Ub+ QbW − RVb

= mgX + FAX + FTX (2.1)

(13)

CHAPTER 2. THEORY Natural Path Planning for the Robird

m ˙V_b+ U_bR − P W_b

= mg_Y + F_A_Y + F_T_Y (2.2)

m

W˙_b+ P V_b− QU_b

= mg_Z+ F_A_Z+ F_T_Z (2.3)

P I˙ _XX− ˙RI_XZ− P QI_XZ+ RQ(I_ZZ − I_{Y Y}) = L_A+ L_T (2.4)

QI˙ _{Y Y} − P R(I_XX− I_ZZ) + (P²− R²)I_XZ = M_A+ M_T (2.5)

RI˙ ZZ− ˙P IXZ + P Q(IY Y − I_XX) + QRIXZ = NA+ NT (2.6)

These equations will be solved numerically and the aerodynamic forces will be estimated from experimental data. In the Robird, there are several controllers incorporated to stabilize the flight and, besides that, the pilot is contantly giving feedback control based on the response of the aircreaft. All tis feedback can be hard to simulate, so, in order to simplify the model, it was assumed that the feedback control for the rotations of the aircraft is perfect. With this assumptions, the moment equations 2.4, 2.5 and 2.6 can be dropped out since any moment on the flight would be met with an opposite moment that would counteract it to maintain the angular velocity and angle set by the user.

2.2.3 Frame transformations

Some of the parameters that are represented in the body-fixed frame have to be represented in the body-fixed frame and vice-versa. To do that, frame transformations are used. One of them is the one that relates the change in time of the of the pitch, yaw and roll angles, represented in the inertial frame with the angular velocities of the body-fixed frame. This relation is given by equation 2.7. [4]





 Φ˙ Θ˙ Ψ˙







=





1 sin Φ tan Θ cos Φ tan Θ

0 cos Φ − sin Φ

0 sin Φ sec Θ cos Φ sec Θ









 P Q R







(2.7)

Another rotation that has to be made is one that represents gravity in the body-fixed frame. In the inertial frame, gravity only has a component in the positive Z-direction. To represent gravity in the body-fixed frame, a rotation using Euler Angles have to be made, obtainin equations 2.8, 2.9 and 2.10. [4]

gx= −g sin Θ (2.8)

gy = g cos Θ sin Φ (2.9)

g_z = g cos Θ cos Φ (2.10)

The final rotation that will be used is to go from the velocity in the body-fixed frame to the velocity in the inertial frame. With this rotation, the displacement in the inertial frame can be calculated from the velocities that have been calculated in the body-fixed frame. The rotation is shown in equation 2.11. After the rotation has been done, the obtained velocity in the inertial frame can be integrated to calculate the displacement of the aircraft. [4]

(14)





 X˙_i

Y˙i

Z˙i







=





cos Ψ cos Θ − sin Ψ cos Φ + cos Ψ sin Θ cos Φ sin Ψ sin Φ + cos Ψ sin Θ cos Ψ sin Ψ cos Θ cos Ψ cos Φ + sin Ψ sin Θ sin Φ − sin Φ cos Ψ + sin Ψ sin Θ cos Φ

− sin Θ cos Θ sin Φ cos Θ cos Φ









 U V W





 (2.11) The rotation of the angular velocity and gravity have to be done while simultaneously with the dynamic equations whereas the rotation of the linear velocity can be done after the simulation has been completed.

2.2.4 Aerodynamic forces

The aerodynamic forces can be divided into three different components, one for each body-fixed axis. The force acting in the x-axis is called the drag force. The form of the drag force is as follows:

F_drag = 1

2ρAC_dU² (2.12)

The drag coefficient depends on many parameters, however, a good approximation is C_d = C_d₁+ C_d₂α.

Similarly, the lift force can be expressed as follows:

F_{lif t}= 1

2ρAC_lU² (2.13)

As with the drag coefficient, the lift coefficient can be approximated as C_l = C_L₁ + C_L₂α. In the Z-direction, there is another force called the lift-induced drag. This drag force that opposes the motion in the body-fixed Z-axis and can be represented with the following expression:

F_L_drag = 1

2ρAC_d_LW² (2.14)

It should be noted that none of the above velocities, U and W , contain a subscript. This is because this speed refers to the airspeed which is a combination of the ground speed and the wind speed. [4, 8]

(15)

CHAPTER 2. THEORY Natural Path Planning for the Robird

Figure 2.1: Definition of the positive force, moment, linear velocity, angular velocity and gravity.

[4]

(16)

Figure 2.2: Definition of the Euler angles Θ, Ψ and Φ. [4]

(17)

Chapter 3

Dynamic model

The goal of the model is to be able to accurately describe the dynamic behavior of the Robird and use that to find out what the most energy efficient path is from point A to point B. To make the model, the differential equations shown in the theory are solved using ode45, a function that uses the Runge-Kutta method of 4^th and 5^th degree. The program solves the differential equations using both methods and compares the results. If the difference between the results is larger than a set tolerance, it decreases the time step and perform the simulation again. It repeats this step until the difference is smaller than the tolerance.

To be able to describe the system, the forces involved have to be described. To achieve that, the results of the simulation are compared with data from test flights and the equations are optimized to fit the data by changing the coefficients.

The model is done for 2D flight. Two different models have been developed, one for the gliding mode and the other one for the climbing mode. The equations used in both models are the same with the only difference being that the aerodynamic and trust forces are different.

In the 2D case, there is no movement in the Y axis and only rotation around the Y-axis is allowed. This means that: P = R = Φ = Ψ = 0. Using these constraints, the Newton-Euler equations for 2D motion become:

m ˙U_b+ QW_b

= mg_X+ F_A_X+ F_T_X (3.1)

m

W˙_b− QU_b

= mg_Z+ F_A_Z+ F_T_Y (3.2)

Using equations 2.8 and 2.9 with the angle constraints of the 2D model and the equations 2.12, 2.13 and 2.14 in equation 3.1 and 3.2 and rearranging them so that the linear accelerations are in the left-hand side, we obtain equations 3.3 and 3.4.

U˙_b= −QW_b− g sin θ − ρA1

2m (C_d₁+ C_d₂α) U²+ FT_X (3.3) W˙_b= QU_b+ g cos Θ −ρA₂

2m(C_L₁+ CL₂α)U²−ρA₃

2mC_d_LW · |W | (3.4) In equation 3.4 the lift induced drag term has the linear velocity in the Z axis multiplied by its absolute value insetead of the velocity square. This is because the drag term always opposes the movement of the aircraft. By doing the multiplication in this way, the value of the square operation is obtained without loosing the information of the direction of movement. This is not done with the velocity in the X-direction because the Robird can only move forward. An event function is added to the program that stops the simulation when the velocity in the X-direction

(18)

becomes negative. It also stops the simulation if the angle of attack becomes larger than 0.6 rads since the Robird would stall at that angle.

Applying the constraints for a 2D case in equations 2.7 equation 3.5 is obtained.

Θ = Q˙ (3.5)

The transformation from the body-fixed to the inertial frame can be reduced using the same constraints that were used to reduce the differential equations to obtain equation 3.6

X˙⁰ Y˙⁰

= cos(Θ) sin(Θ)

−sin(Θ) cos(Θ)

U W

(3.6)

3.0.1 Gliding

For gliding mode in 2D, the previous equations 3.1 and 3.2 hold when the thrust force is set to zero. The equations are put into MATLAB to be solved simultaneously with the ODE solver.

The results will be obtained in the body-fixed frame so, after getting the results, equation 3.6 is used to get the velocities in the inertial frame. Finaly, the velocities in the inertial frame can be integrated to get the displacement of the aircraft. Section 7.2 shows the function of the gliding mode.

3.0.2 Powered flight

To model the flapping mode, the same approach as in gliding mode was taken but adding an additional force term in equation 3.1 that represent the thrust force. This force was assumed to be constant during the simulation. The pitch angle and angular velocity were assumed to be constant as well. As in the gliding mode, the velocity is obtained in the body fixed frame, it can be transformed to the inertial frame using equation 3.6 and then, the velocity can be integrated to obtain the displacement in the inertial frame. Section 7.1 shows the function for the climbing mode.

3.1 Optimization

Optimizing a differential equation is the process of finding out the values of the coefficients of the equation that make the equation fit the data more accurately. To optimize the equations, a program was made. The program uses the function lsqcurvefit to change the lift and drag coefficients to minimize the error between the simulated data and the real one.

The experimental data used in the optimization process was the velocities in the inertial frame.

The velocity in the X-direction was obtained from the GPS sensor and the velocity in the Z-direction was obtained by differentiating the altitude. Before processing the data, it was in- terpolated so that the time steps of the different data sets were the same and calculations could be done between them. For the pitch angle, the average of the pitch angles was computed and assumed constant. As the initial conditions for the simulation, the first 100 values of the inter- polated experimental data were averaged and that value was used. This was done to reduce the effect of any disturbance or noise.

To be able to objectively compare the results of the simulations, the mean absolute error (MAE) was used. This error is calculated by adding the absolute value ot the error between the simulation and the experimental data at each point and dividing by the number of points as shown

(19)

CHAPTER 3. DYNAMIC MODEL Natural Path Planning for the Robird in equation 3.7. The absolute value is computed so that the positive and negative errors do not cancel each other out.

M AE = PN

i=1|y_i− x_i|

N (3.7)

In equation 3.7, yi is the experimental data, xi is the simulated value and N is the number of data points. The data used for the error calculation was the velocity in the X and Z-inertial frames.

3.2 Results

3.2.1 Gliding

Two different data sets from test flights were used. The data sets were of the Robird in gliding mode moving in a 2D plane. There were some deviations and rotations in other axes but they were considered small enough to be neglected. The simulations were individually optimized to fit the different data sets and then, the results were compared. However, the sets of coefficients obtained to match the different data sets were not the same. Figure 3.1 shows a comparison between the experimental data and the simulated one with the coefficients optimized for that data set. It can be seen that the results of the simulation are similar to the experimental data but without the deviations. However, this deviations are caused by sudden changes in the ex- ternal conditions, changes in angles that were not considered in the simulation and by the fact that the real flight had some components in the Y-axis that were neglected.

Figure 3.1: Comparison of the experimental data of the first set with the simulated data with the constants optimized for the first set.

Figure 3.2 shows the comparison of the data of the second data sets with the simulation results done with the coefficients for the first data sets. It can be seen that the results do not match.

This is because the constants obtained for each data set were different. A set of constants that satisfied both experimental results was not found. It was tried to optimize both data sets at the same time, minimizing the error for both functions but, by doing this, the data did not match any of the experimental results.

To obtain an objective value with which the results of the simulations could be compared, the mean absolute error (MAE ) was computed. To have a useful simulation, not only the data set

(20)

Figure 3.2: Comparison of the experimental data of the second set with the simulated data with the constants optimized for the first set.

C1 C2

Ui Wi Ui Wi

Set 1 0.4824 0.6470 3.3219 1.8975 Set 2 1.9629 0.5982 0.1800 0.5716

Table 3.1: Mean absolute error of simulations experimental data for the gliding mode. They give a measure of the validity of the simulation.

with which the constants were calculated has to be fitted, it should also fit any oder data set that meets the conditions for which the simulation was designed. So, to decide which constants work better, the MAE for the velocities in the X and Z-inertial frames were computed for the two sets of coefficients in the two different data sets. The results of all the MAEs can be seen in table 3.1. The lower the value of the MAE, the better the fit for a specific data set. It can be see that the constants for the second data set fit better the data set with which they were obtained than the constants obtained with the first data set. However, looking at how the simulations perform with the data set that was not used to calculate them, the first constants work better than the second ones by a large margin. So, it was decided to use the constants obtained with the first data set since, overall, they perform better. The selected coefficients were C_d₁ = 0.0424, Cd2 = −0.2395, CL1 = −0.0794, CL2 = 0.3813 and Cd_L = 28.4118.

There are two main possible reasons why the sets of coefficients do not match. The first reason is that some of the data is wrong. It could be that one of the sensors was not giving the right measurement and, because of that, the data to be fitted should be discarded. This is unlikely but should be considered. This first reason also covers the fact that the wind speed was not considered at all. The lift and friction forces are both affected by the airspeed of the aircraft, this is the combination of the aircrafts speed with respect to the inertial frame and the wind speed. Since the Robird with which the measurements were done did not have an airspeed sensor, it was assumed that the wind at the time of the measurements was so low that it could be neglected. This was considered valid because the pilot had measured a low wind speed in the ground, however, that does not necessarily mean that the airspeed at a higher altitude is also low. The fact that the wind-speed was not considered is believed be the reason of the discrepancy in the obtained coefficients. More experimental data should be obtained including the wind-speed to try the optimization process again. Finally, the last possible reason is that the model that is being fitted is wrong. It could be possible that the equations that are being

(21)

CHAPTER 3. DYNAMIC MODEL Natural Path Planning for the Robird fitted is not of the form of the real dynamic equation and there is a parameter that has been neglected that should be considered. If this was the case, those parameters would have to be considered in the lift and drag forces. However, this is unlikely since the dynamic equations that were used are widely used to describe the movement of aircrafts.

3.2.2 Powered flight

The model for the climbing mode was obtained using the same approach as in the gliding mode but with one more constant to find which represents the thrust force. Two data sets were used to find the constant the dynamic system. The results were similar to the ones in the gliding mode, the constants obtained with one data set do not fit the other data set. Table 3.2 shows the MAE for all compared sets. It can be seen that the constants for the first data set perform better than the ones for the second one. The selected coefficients were Cd1 = −0.0082, Cd2 = −0.0830, C_L₁ = 0.6730, C_L₂ = 0.3987 , C_d_L = 25.8872 and F = 8.9438. This last value represents the force and is multiplied by the throttle level, an number between zero and one.

The reasons why the obtained constants do not fit both data sets are the same as those in the climbing mode. However, there is one more reason to add. During the derivation of the dynamic equation, it was assumed that the aircraft could be considered a rigid body but, during powered flight, the wings flap so it is not a rigid body anymore. It was assumed that, since the overall movement was of interest, the oscillatory behavior was not to be considered and instead could be replaced by a continuous constant force. However, it was never checked if this assumption holds. A way of solving this problem would be to use multi-body dynamics to simulate the behavior during the climbing mode. The aircraft would be divided into three different parts that move relatively to each other, then the kinematics of this movements could be derived as well as the dynamic behavior of each individual body. Once the equations are obtained, they can be combined to form a dynamic equations that describes the entire aircraft. [7]

C₁ C2

Ui Wi Ui Wi

Set 1 0.3170 0.7314 1.7930 0.7102 Set 2 4.1840 1.8199 0.2859 0.6495

Table 3.2: Mean absolute error of simulations experimental data for the climbing mode. They give a measure of the validity of the simulation.

(22)

(23)

Chapter 4

Path Planner

The path planner is a program that takes the models for the climbing and gliding modes that were developed earlier and, together with an initial point A a final one B and the initial velocity calculates what is the most energy efficient path that connects both points. The path planner is developed for 2D flight. The inertial frame is chosen such that both points are included in the XZ-plane and gravity only has a component in the Z-axis. The positive directions of the frame are the same as in the rest of the report. The following paragraphs in the report explain how the program decides what the most energy efficient path is. In the explanation, the lines of the program in which the explained action is taking place are specified. The program can be found in Appendix A ????.

The first thing that the program has to do is to check if it is possible to reach the final point only using 2D flight modes. For that, the slope of the line that goes through the initial and final points is calculated and compared with the maximum climbing and gliding slopes. To be able to use 2D flight modes, the slope has to be between those two values, otherwise, the point is only reachable using 3D flight mode. In the last case, the program will give an error message indicating that it is not possible to reach the point (lines 10-22). And if not it will continue with the program.

The models developed in this model can only go forward since the robird is not able to do loops or fly upside down. However, it was decided that the model should be able to go in both directions. To solve this problem, the program checks the direction in the X-axis to which it has to go and, if it is negative, it mirrors the initial and final points arround the origin of the axis to perform the simulations with a positive direction and, when the path is completed it mirrors the solution back to show the right path (lines 26-32).

Gliding does not consume any energy since it was assumed that the energy of consumed by the flight control boards during operation could be neglected and the same was true for the energy consumed by the servos. So, the most energy efficient mode is gliding. It was also discovered that there is one angle at which the gliding ratio is maximize. So it makes sense that, for an energy efficient flight, if possible, the gliding mode with the best angle is used.

So, to apply this in the path planner, an imaginary line is drawn that goes through the desti- nation point and has the slope of the best gliding ratio that the Robird can achieve (lines 34-50).

For the climbing mode, the energy consumption depends on the needed thrust force, the pitch angle and the time that it takes. To find the most energy efficient path, the program first selects N angles evenly spaced from the minimum to the maximum pitch angle and M forces also evenly spaced and between the minimum and maximum thrust forces (lines 52-68). Then, using two nested for loops it runs through all the combinations of the generated angles and thrusts (lines

(24)

83-156). The climbing mode runs until it reaches the line for gliding mode. It also checks if the maximum climbing angle is ever surpassed.

After the climbing mode is finished, it is changed to the gliding mode. For this change it is assumed that it occurs instantaneously conserving the same velocity in the inertial frame. In reality, there is a transition between modes where a moment is applied that rotates the aircraft to go from the climbing angle to the gliding angle. However, this transition is really short compared to the total flight so it is assumed that it can be neglected. The gliding mode is run until the displacement in the X-axis of the simulation coincides with X-coordinate of the final point.

At that point the Y-coordinate is compared and, if it is larger than the one of the final point, the simulation is repeated but with the climbing mode finishing earlier and, if the opposite occurs, the simulation is run again but with the climbing mode lasting longer. This process is repeated until the Y-coordinate is above the final point with an allowed error. It was decided to only allow the error above the final point because then, in that case, the system has more energy than needed so, with the help of a controller, the target point could be reached without using more energy from the battery. (lines 158-254)

The final step is to assign a total energy consumption to each path. For that the power needed in the climbing mode is given based on the applied thrust and the gliding mode is assumed to not use any power. By multiplying the power by the time that it was used, the consumed energy can be calculated. (lines 278-296)

For the last step of adding the power consumption in the flight, a function is needed that relates some conditions of the flight to the energy consumption. To this end, some data was collected from test flights as shown in table 4.1. There are six different data sets with slightly different conditions. The data shown in the table is the average value over a time span with relatively constant conditions. The first four data sets have a thrust input of 100% and the last one of 75%. It can be seen that the first data sets with a maximum thrust input have different power consumptions ranging from 110W to 137W, however, it was assumed to be constant and the average of the five measurements was set as the power consumption at maximum thrust. The last measurement is done with a thrust input of 75% and has a power consumption of 79W.

So, by making the assumption that the power consumption only depends on the thrust input, three data points are obtained, the power at 100, 75 and 0% thrust which are 126, 79 and 0W respectively. This is not enough to determine the shape of the function so it was decided to use three different functions that determine the power consumption and see how they affect the final energy consumption. The three functions used where a linear function that goes through the zero and 126W points, a second degree polynomial that goes through the three experimental points and a third degree polynomial that goes through the three points as shown in figure 4.1

Data Set Climb rate Pitch angle Power Throttle

SET 1 2.06 0.1439 134 1

SET 2 1.65 0.0625 115 1

SET 3 2.37 0.1520 137 1

SET 4 2.25 0.1860 134 1

SET 5 -0.05 0.0532 110 1

SET 6 0.56 0.0989 79 0.75

Table 4.1: Average of different parameters in six data sets from test flights.

(25)

CHAPTER 4. PATH PLANNER Natural Path Planning for the Robird

Figure 4.1: Different functions for the power depending on the input thrust in climbing mode.

The functions used are a linear function, a second degree polynomial and a third degree polynomial.

4.1 Results

Even though the constants for the flight modes were not completely correct, they were still used in the path planner. The idea being that, even though the results of the path planner are not accurate, they should be somewhat close to the real solution. The program was designed in such a way that the models could be easily changed once better results were found. With a better model, the maximum climbing angle and best gliding slope would have to be changed in the path planner program. Besides that, the model to of the dynamic simulation can be changed to any model that has the thrust and pitch angle as an input together with some initial conditions and a time span for the simulation. The model would also have to give the velocities in the body fixed frame as the output.

Figure 4.2 shows the initial and final points represented with a blue and a red star respectively as well as eleven different paths that were obtained. Every path has a black cross that shows the point at which the flight mode was changed from climbing mode to gliding mode. The time at which the mode is changed differs in each path. It can be seen that the gliding mode at the beginning has different trajectories but, at the end, have the same slope and go parallel. This is because each path has different initial conditions for the gliding mode depending on how the climbing mode ended but, after some time, the forces involved balance each other out achieving a constant speed in both coordinates. This constant speed that is reached after some time always have the same value and result in parallel paths.

2 3 6 7 8 10 11 12 14 15 16

Pitch -0.133 0.133 -0.133 0.133 0.4 -0.133 0.133 0.4 -0.133 0.133 0.4

Thrust 0.2 0.2 0.467 0.467 0.467 0.733 0.733 0.733 1 1 1

Table 4.2: Thrust and pitch angle of the climbing mode for different simulations. The numbers on top are the number of the simulation that are shown in the legend of figures 4.2, 4.3 and 4.4.

(26)

Figure 4.2: Different paths obtained with the path planner. The numbers in the legend are related to the conditions of the climbing mode by table 4.2. The blue star represents the initial position and the red one the final position. The black crosses at each path show the point at which the mode was changed from climbing to gliding.

Figure 4.3 shows the velocity in the X-direction of the different paths. It is important to note that there is a nonlinearity in the velocity at the time when the modes have been change. This occurs because the inertial velocity is assumed to be constant, however, the angle of the aircraft is suddenly changed from the angle at which it was climbing to the fixed gliding angle. This nonlinearity could be solved by adding an intermediate mode to change from one mode to the other one. Another option would be to add a fixed angular velocity at the beginning of the gliding mode for a certain amount of time such that the gliding mode goes smoothly from the climbing angle to the gliding angle.

Figure 4.4 shows the energy consumption of every path that was obtained. The energy shown in this graph was calculated using the linear function of figure 4.1. It can be seen that the energy increases with four different slopes. This is because the power consumtion only depends on the thrust so, in the paths that used the same thrust, the energy will increase at the same rate. So, in paths with the same input thrust, the energy consumption will only depend on the time at which the climbing mode finishes. Since the only variable that remains is the pitch angle, this variable will determine the time at which it will stop. By looking at the simulations 6, 7 and 8 in figure 4.4, which have the same input thrust, it can be seen that even though the pitch angles are in an increasing order, the energy consumption is smaller in simulation 7. This means that between simulation 7 and 8, there is a minimum.

To show how the energy changes as a function of the thrust and pitch of the climbing mode, a surface plot was made. Three different plots were made with the different power functions of figure 4.1. Figure 4.5 shows the energy as a function of the pitch and thrust with the power linearly related to the thrust. Figure 4.6 shows the same but with the power being a related to the thrust with a second degree polynomial. And figure 4.7 is the same but with the power related to the thrust with a third degree polynomial. It can be seen that for figure 4.5 and 4.6 the region with the most energy efficient path is quite similar, with a low thrust, low pitch being the best option. However, in figure 4.7 the region with the most energy efficient path is at a higher thrust and pitch. The best thrust and pitch angle for the linear power and thrust relation are 0.2 and −0.061rad respectively with an energy consumption of 315.6J . For the

(27)

Figure 4.3: Velocity of the aircraft in the X-axis of the inertial frame during the different paths.

The numbers in the legedn are related to the throttle and pitch angle of the climbing mode as shown in table 4.2.

second degree polynomial the same thrust and pitch angles should be used but, this time, the energy consumption is of 192.2J . For the third degree polinomial, the thrust should be 0.648 with a pitch angle of 0.17rad to achieve an energy consumption of 546J .

By looking at figure 4.4it was understood that at every thrust level, there is a pitch angle that takes less time to finish the climbing mode. If the power is only dependent on the thrust, then, for every thrust level, there will be a pitch angle that will always produce the most energy efficient path at that thrust level. The pitch angle that takes the least amount of time to change to gliding mode for every give thrust can be seen in figure 4.8. Figure 4.9 shows the time that the climbing mode was on at the most effective pitch angle. By multiplying the function in figure 4.9 by the function that relates the power with the input thrust, the energy consumption is obtained. It was observed that the shape of the function in figure 4.9 does not change for paths in which the climbing mode is on for long enough, so those in which every path has time to reach the final velocity. So, even though the energy cannot be estimated without running the whole simulation, an estimation of the best pitch angle and thrust can be obtained combining figures 4.1, 4.8 and 4.9. To do that, the function that gives the power depending on the thrust would have to be multiplied by function that relates the duration of the simulation as a function of the time. This operation would result in the energy consumption as a function of the thrust.

The operation was done for the power functions in figure 4.1 and the time function in figure 4.9 obtaining the energy function in figure 4.10.

In order to check the a change in one of the parameters of the model had a large effect in the result, three different coefficients were increased by a 30% one at a time. the changed parameters were the first lift coefficient C_L₁ of the climbing mode, the first drag coefficient C_d₁ of the climbing mode and the second lift coefficient CL2 of the gliding mode. It was observed by looking at the surface plots that the changes were not significant. The flight conditions for the minimum energy consumption path do not change much as can be observed in table 4.3.

(28)

Figure 4.4: Energy consumption of every path. To calculate the energy of the graph, the linear power consumption and thrust relation in figure 4.1 was considered for the climbing mode and, for the gliding mode it was assumed to be zero.

Figure 4.5: Surface plot that shows the relation between the pitch angle and thrust during the climbing mode and the energy consumption. The power consumption is linearly related to the thrust and can be seen in figure 4.1.

(29)

Figure 4.6: Surface plot that shows the relation between the pitch angle and thrust during the climbing mode and the energy consumption. The power consumption related to the thrust by second degree polynomial and can be seen in figure 4.1.

Figure 4.7: Surface plot that shows the relation between the pitch angle and thrust during the climbing mode and the energy consumption. The power consumption related to the thrust by third degree polynomial and can be seen in figure 4.1.

(30)

Figure 4.8: Pitch angle that gives the shortest time in gliding mode for a given thrust level.

Figure 4.9: Shortest time of the climbing mode for a given thrust level.

(31)

Figure 4.10: Energy consumption as a function of the thrust for the different functions of figure 4.1

Thrust Pitch Energy

Linear 0.2 -0.0881 303.1

Climb C_L₁ Second degree 0.2 -0.0881 156.8

Third degree 0.68 0.1288 555.5

Linear 0.2 -0.0746 312.5

Climb C_d₁ Second degree 0.2 -0.0746 148.5

Third degree 0.634 0.1424 522.5

Linear 0.2 -0.061 331.4

Glide C_L₂ Second degree 0.2 -0.061 157.4

Third degree 0.688 0.17 544.7

Table 4.3: Values for the thrust and pitch angle for the most energy efficient path for the three different functions of the power consumption of figure 4.1 with the energy that is consumed.

There are three different sets corresponding to an increase of a 30% of the value of the coefficient that is shown at the right of the table.

(32)

(33)

Chapter 5

Discussion and conclusion

In this report, two main things were done, first a dynamic model describing the behavior of the Robird was developed and, then, a path planner was designed that takes the developed models and chooses what the most energy efficient path from point A to B is.

To develop the dynamic model, the well known dynamic equations of a rigid aircraft were used.

Using the theoretical form of the aerodynamic forces and combining them with the dynamic equations, the equations for a 2D model were obtained. To simplify the model, it was assumed that the controller on board of the aircraft would perfectly obtained the pitch angle wanted by the user. With this last assumption, only two equations were left. Using experimental data, the drag and lift coefficients were determined.

In the gliding mode, two different data sets obtained in test were used to find the coefficients.

Using the optimizing function lsqcurvefit, two different sets of coefficients were found. The mean absolute error (MAE ) was used to compare the effectiveness of the different coefficients. It could be seen that neither of the coefficients was able to accurately fit both data sets, however, one of the data sets had a better fit and was used for the rest of the report. The reason why it was thought that the coefficients did not fit was due to a lack of data. It is believed that the air has a big influence in the dynamic behavior of the aircraft, however, the Robird with which the test flight was done did not have a air speed sensor so it was assumed that the wind velocity was zero.

The same approach was used to fit the climbing mode model. Two experimental data sets were again used to find the lift and drag coefficients as well as the thrust force. The MAE was again used to compare the accuracy of the results. The results again did not match both data sets, however, one of the constants set had a lower error for both experimental data sets and was used for the path planner. On top of the absence of the wind speed data as a cause of errors, in the climbing mode, the equation that is being fitted could not be right. The equations to be fitted were deriving by using a rigid-body assumption. However, if the wings move, it cannot be considered a rigid body. It was assumed that, for the overall trajectory, the oscillations could be neglected and the forces could be exchanged by a total force so the rigid-body model could be used. This assumption may not hold. It was proposed to use a multi-body dynamic model to make the climbing mode model. Such a model does not have the rigid body assumption and, thus, would hold for a flapping flight.

With the dynamic models already optimized, the path planner was developed. The program has an initial velocity and point and a final point as input. The program goes through a grid of input thrust and pitch angles to compute different climbing paths. Then it goes to the gliding mode and changes the time at which the modes have to be changed to find the most energy efficient path. Then, the program calculates the energy used during the flight. The relation between the

(34)

energy consumption and the thrust was not known so three different functions were tried. It was shown that if the power consumption only depend on the thrust, there is one pitch angle, the one for the path with the least time in climbing mode, for each thrust input that will always have the minimum energy consumption for that thrust level. If that is the case, the time that the climbing mode last can be shown as a function of the thrust level. It was observed that the form of this function does not change if all the calculated paths have enough time to reach the steady state velocity for climbing. So, by multiplying the power function by this last function, the minimum energy consumption at each thrust level can be found. The lowest value of this function will be the thrust at which the energy consumption is minimum. And that thrust can be related to a specific pitch angle.

As a conclusion, it can be said that a dynamic model was developed to describe the behavior of the Robird. However, the coefficients were not successfully matched. It is believed that for the gliding mode, the main reason for the mismatch is that there was not data available for the wind speed so it was not considered. It is believed that with that data, the right coefficients should be found. For the climbing mode, besides the lack of wind speed data as a cause of error, it could be that a more sophisticated model is needed to describe the behavior of the aircraft.

This is because a rigid-body assumption was made when deriving the model and, since there are two problems that can be causing the mismatch, it is not possible to say if the rigid-body assumption is valid during the flapping mode. The usage of a multi-body dynamic model was proposed if the rigid body assumption was proven not to be valid.

A program to find the most energy efficient path was developed. The obtained dynamic models were used in the path planner even though it is known that they are not completely correct.

Since the power consumption was not known it was assumed that it depended on the throttle level and three different functions were analyzed. A surface plot was used to show how the energy depended on the thrust level and pitch angle of the climbing mode. The surface plot was broken into three plots using the assumption that the power consumption only depends on the thrust, making it easier to visualize how the energy consumption depends on the thrust and pitch. It was also observed that, for long enough flights, the function that relates the energy consumption with the trust level has the same form, meaning that the same pitch angle and thrust level would give the most energy efficient flight. Finally, if was observed that the conditions for a efficient flight did not change much with relatively large changes in the coefficients. This means that, even though the models were not quite right, a similar result is expected for the right models.

(35)

Chapter 6

Bibliography

[1] D. Parrot J. Allan J. Bishop H. McKay. “Review of international research literature re- garding the effectiveness of auditory bird scaring techniques and potential alternatives”. In:

(2003).

[2] N. Nijenhuis C.H. Venner S. Sramigioli G. A. Folkertsma W. Straatman. “Robird: A Robotic Bird of Prey”. In: (2017).

[3] R. Dvoˇr´ak. “Aerodynamics of bird flight”. In: (2016).

[4] M.R. Napolitano. Aircraft Dynamics: From Modeling to Simulation. 2012.

[5] J. Tang D. Viieru H. Liu W. Shyy Y. Lian. Aerodynamics of low reynolds number flyers.

2008.

[6] S.K. Chimakurthi P. Trizila C.K. Kang C.E.S. Cesnik H. Liu W. Shyy H. Aono. “Recent progress in flapping wing aerodynamics and aeroelasticity”. In: (2010).

[7] Tahir Rashidr. The Flight Dynamics of a Full-Scale Ornithopterl. 1995.

[8] J. Grauer. Modeling and System Identification of an Ornithopter Flight Dynamics Model.

2012.

(36)

Energy Efficient Path Planning for the Robird