Module/guidance vector field

From PaparazziUAV
Jump to: navigation, search


We have written an algorithm for solving the problem of tracking smooth curves by an unmanned aerial vehicle travelling with a constant airspeed and under a wind disturbance. The algorithm is based on the idea of following a guiding vector field which is constructed from the implicit function that describes the desired (possibly time-varying) trajectory.

For fixed-wings the output of the algorithm can be directly expressed in terms of the bank angle of the UAV in order to achieve coordinated turns. Furthermore, the algorithm can be tuned offline such that physical constraints of the UAV, e.g., the maximum bank angle, will not be violated in a neighborhood of the desired trajectory.

The implementation for rotorcraft is still in the TO DO list.

This work is based on the paper Guidance algorithm for smooth trajectory tracking of a fixed wing UAV flying in wind flows, presented at ICRA 2017. A more rigorous and detailed analysis can be found in the paper A guiding vector field algorithm for path following control of nonholonomic mobile robots.

We can further exploit the convergence properties of the fixed-wing aircraft to the desired trajectory for controlling formations of several aircraft. Check the subsection python Apps. A detailed analysis can be found in the paper Circular formation control of fixed-wing UAVs with constant speeds (submitted to IROS 2017). A video about how the circular formation algorithm works .

Here we show an animation of an airplane following an ellipse and a sinudoidal track.


How does the algorithm work?

We start from the implicit equation of the trajectory, for example for a circumference we have that \varphi(x,y) = x^2 + y^2 - R^2, where x, y and R are the x, y coordinates with respect to HOME (O in the following figure) and R is the radius of the circumference. Note that when the vehicle is over the desired trajectory, then \varphi(x,y) = 0, otherwise it will be different from zero. We will use this concept of "level set" to define the notion of distance to the trajectory, namely e := \varphi(x,y). It has to be noted that this is different from the Euclidean distance (see figure below).

For the sake of clarity we will consider just 2D trajectories and parallel to the ground. Let us stack the x and y coordinates in p := \begin{bmatrix}x & y\end{bmatrix}^T and consider the following kinematical model for our vehicle

\begin{cases} \dot p &= sm(\psi) + w \\ \dot\psi &= u, \end{cases}

where s\in\mathbb{R}^+ is a constant that can be considered as the airspeed, m = \begin{bmatrix}\cos(\psi) & \sin(\psi)\end{bmatrix}^T with \psi\in(-\pi, \pi] being the attitude yaw angle, w\in\mathbb{R}^2 is a constant with respect to O representing the wind and u is the control action that will make the vehicle to turn. We also notice that the course heading \chi\in(-\pi, \pi], i.e. the direction the velocity vector \dot p is pointing at, in general is different from the yaw angle \psi because of the wind.


Given a desired trajectory, we compute its normal n(p) := \nabla \varphi(p) and its tangent \tau(p) = En(p), \quad E=\begin{bmatrix}0 & 1 \\ -1 & 0\end{bmatrix}. The idea is to steer the vehicle such that it follows the direction given by the unit vector calculated from

\dot p_d(p) := \tau(p) - k_e e(p)n(p),

where k_e is a positive gain. At each point p we can build a unit vector from \dot p_d(p). This collection of vectors is our Guidance Vector Field. Note that when the error is zero, then we are just tracking the tangent to the trajectory.

We assume that the trajectory is twice differentiable (so its gradient and Hessian exists) and that it is regular in a neighborhood of \mathcal{P}, i.e.

\nabla\varphi(p) \neq 0, \quad p\in \mathcal{N}_{\mathcal{P}},\quad \mathcal{N}_{\mathcal{P}} := \{p \, : \, |\varphi(p)| \leq c^* > 0\}.

The algorithm needs for work the measurements of \dot p (typically derived from GPS), position p w.r.t. HOME (typically derived from GPS) and yaw angle \psi or sideslip angle \beta = (\psi - \chi). The algorithm still works quite well in practice with only \dot p and p, in other words, you will only need to have installed a GPS. The latter is the implemented version in pprz/master.

Using the GVF in the flight plan

First of all, do not forget to include in your airframe or flightplan the module. I would recommend to do it in the airframe.xml, since each trajectory has its own sets of gains and these gains depend on the aircraft.

<module name="gvf_module">

and in your telemetry file (already in default_fixedwing_gvf.xml)

<message name="GVF" period="1"/>

List of functions available for the user ready to be called from the flightplan.


<call fun="gvf_ellipse_wp(WP_mywaypoint, gvf_ellipse_par.a, gvf_ellipse_par.b, gvf_ellipse_par.alpha)"/>
<call fun="gvf_ellipse_XY(float x, float y, gvf_ellipse_par.a, gvf_ellipse_par.b, gvf_ellipse_par.alpha)"/>

Straight lines

<call fun="gvf_line_XY_heading(float x, float y, gvf_line_par.heading)"/>
<call fun="gvf_line_XY1_XY2(float x1, float y1, float x2, float y2)"/>
<call fun="gvf_line_wp1_wp2(WP_mywaypoint1, WP_mywaypoint2)"/>
<call fun="gvf_segment_XY1_XY2(float x1, float y1, float x2, float y2, float gvf_line_par.d1, gvf_line_par.d2)"/>
<call fun="gvf_segment_wp1_wp2(WP_mywaypoint1, WP_mywaypoint2, float d1, float d2)"/>
<call fun="gvf_line_wp_heading(WP_mywaypoint, gvf_line_par.heading)"/>

Watch out, `line' means that the aircraft tracks an infinite line, while `segment' means the aircraft will track the line between two points. The values 'd1' and 'd2' mean how far the aircraft will go beyond the end of the segment before turning back (useful when it is windy and you want to be sure that your aircraft tracks the segment by converging to it before been in the middle of the two points).


<call fun="gvf_sin_XY_alpha(float x, float y, gvf_sin_par.alpha, gvf_sin_par.w,, gvf_sin_par.A)"/>
<call fun="gvf_sin_wp_alpha(WP_mywaypoint, gvf_sin_par.alpha, gvf_sin_par.w,, gvf_sin_par.A)"/>
<call fun="gvf_sin_wp1_wp2(WP_maywaypoiny1, WP_mywaypoint2, gvf_sin_par.alpha, gvf_sin_par.w,, gvf_sin_par.A)"/>


For setting the tracking direction (the direction of the tangent vector \tau in the above figure)

<call fun="gvf_set_direction(s)"/>

where s is an integer number 1 or -1.

Gain tuning

Each trajectory has its own set of two gains k_n, k_e. The former one determines the convergence rate for aligning the vehicle to the vector field. A typical value for starting tuning k_n should be between 0.2 and 1. The latter gain determines how aggressive is the vector field. For example, in the following figure we have an ellipse with two different values of k_e, at the left we have a value of 3 and at the right 0.4. Note how in left one the vectors are "more aggressive" towards the trajectory. While a big value can make the vehicle to converge quickly to the trajectory, it can make the tracking unstable once the vehicle is close it. This is because the vector field might change so quick that physically the vehicle cannot follow it. Check the Section IV in the original paper in the introduction.

The gains can be set in the settings file of the module or in the airframe file. For example

<define name="GVF_ELLIPSE_KE" value="1"/>

Python Apps

Visualizing the desired trajectory and the aircraft within the vector field

Until the "visual integration" of the GVF into the ground station is ready, one can track the vehicle and the trajectory by using the python script at "$PAPARAZZI_HOME/sw/ground_segment/python/gvf/gvfApp ac_id", where ac_id is the ID of the fixedwing to be tracked by the App.

Circular formations

One can synchronize or arrange a group of fixed-wing aircraft (assuming they have equal ground speed) in a desired circle by using the python script at "$PAPARAZZI_HOME/sw/ground_segment/python/gvf/gvfFormation". For running the script, all your aircraft MUST follow a circle employing the guidance vector field, i.e.,

<call fun="gvf_ellipse_wp(WP_CIRCLE, gvf_ellipse_par.a, gvf_ellipse_par.b, gvf_ellipse_par.alpha)"/>

the values of gvf_ellipse_a, gvf_ellipse_b will be modified by the formation control script, and the value of gvf_ellipse_alpha is irrelevant. Note that the aircraft do not need necessarily to orbit around the same waypoint.

We need three text files in order to run the script:

  • ids.txt , which labels the aircrafts. For example, if the file contains '14 56 34', then the aircraft with ids 14, 56 and 34 are labeled with the tags 1, 2 and 3 respectively.
  • topology.txt , which defines the links between the aircraft. The file contains a matrix where the number of rows is the number of aircraft, and the number of columns is the number of desired links. A link is defined by setting 1 and -1 in a column. Following the example, the matrix \begin{bmatrix}1 & 0 \\ -1 & 1 \\ 0 & -1\end{bmatrix} defines two links. The first one between the aircraft 1 and 2, and the second one between the aircraft 2 and 3.
  • sigmas.txt , which defines the desired angles \sigma_k between aircraft for each link k. Following the example, if the file contains '0 90', then the desired angle for the first link is 0 degrees (e.g., the aircraft 1 and 2 will meet each other) and the second desired angle will be +90 degrees (positive is clock-wise).

Some useful tips:

  • How you place 1 and -1 in the topology matters in the following sense. In our example, for the second link, the +90 degrees is from the aircraft 2 to the aircraft 3. If the topology were \begin{bmatrix}1 & 0 \\ -1 & -1 \\ 0 & 1\end{bmatrix}, then the +90 degrees will be from the aircraft 3 to the aircraft 2.
  • Watch out for "impossible" configurations, e.g., do not set up loops in the topology such as 1 -> 2 -> 3 -> 1 \begin{bmatrix}1 & 0 & -1\\ -1 & 1 & 0 \\ 0 & -1 & 1\end{bmatrix}, so the desired sigmas could be asking for an impossible configuration in the circle. Avoid loops! even if the desired configuration is ok, they can introduce undesired equilibria in the system, i.e., different configurations.
  • The aircraft do not need to share the same center or waypoint. This script just control the sigmas.

The script needs two additional scalar values:

  • radius , sets the desired radius of the circular formation, i.e., once the aircraft are synchronized, they will be orbiting the waypoint with this desired radius.
  • k, this is the gain of the control algorithm. As an example, for three aircraft with a desired circumference of 80 meters, a value of k=10 works ok.

The algorithm works by setting different radii to the aircraft. If two aircraft orbit a waypoint with different radius (assuming equal ground speeds), then the one with bigger radius travels more distance in a completed lap. We exploit this fact in order to control the different \sigma_k. We refer to Circular formation control of fixed-wing UAVs with constant speeds for details about the control signal and how to calculate an appropriated gain k.

In the following figure we show the rendezvous of three aircraft in the same circle by using this script.


Defining your own trajectory

The algorithm is placed in $PAPARAZZI_HOME/sw/airborne/modules/guidance/gvf . The trajectories are placed in $PAPARAZZI_HOME/sw/airborne/modules/guidance/gvf/trajectories .

Here there is an example (the straight line) about how to code your own trajectory.



// Straigh line
extern bool gvf_line_XY_heading(float x, float y, float heading);
extern bool gvf_line_XY1_XY2(float x1, float y1, float x2, float y2);
extern bool gvf_line_wp1_wp2(uint8_t wp1, uint8_t wp2);
int out_of_segment_area(float x1, float y1, float x2, float y2, float d1, float d2);
extern bool gvf_segment_XY1_XY2(float x1, float y1, float x2, float y2, float d1, float d2);
extern bool gvf_segment_wp1_wp2(uint8_t wp1, uint8_t wp2, float d1, float d2);
extern bool gvf_line_wp_heading(uint8_t wp, float heading);



bool gvf_line_wp_heading(uint8_t wp, float heading)
  heading = heading * M_PI / 180;

  float a = waypoints[wp].x;
  float b = waypoints[wp].y;

  gvf_line_XY_heading(a, b, heading);

  return true;

The idea is that while you can define a straight line in many different ways, you make all the transformations in the function called by the user such that the algorithm is always evaluating the trajectory in the same way. It helps for maintaining the code and for tuning the gains (you do not want to have different set of gains depending on how you call your trajectory).

We finally call the function gvf_line_XY_heading that will invoke the algorithm.



bool gvf_line_XY_heading(float a, float b, float heading)
  float e;
  struct gvf_grad grad_line;
  struct gvf_Hess Hess_line;

  gvf_trajectory.type = 0;
  gvf_trajectory.p[0] = a;
  gvf_trajectory.p[1] = b;
  gvf_trajectory.p[2] = heading;

  gvf_line_info(&e, &grad_line, &Hess_line); =;
  gvf_control_2D(1e-2 *,, e, &grad_line, &Hess_line);

  gvf_control.error = e;

  horizontal_mode = HORIZONTAL_MODE_WAYPOINT;
  gvf_segment.seg = 0;

  return true;

The three first lines are the \varphi(p), \nabla\varphi(p) and H(\varphi(p)), where H stands for the Hessian. These three guys will be populated by gvf_line_info (in a moment we will be see how to write it).

Then gvf_traj_type stands for the kind of trajectory. Right now we have 0 for lines, 1 for ellipses and 2 for sinusoidals. Feel free to choose an index that has not been taken. The rest gvf_param[x]'s are the parameters of the trajectory. You have up to seven for describing a trajectory.

The function gvf_control_2D executes the algorithm for calculating the desired turn for the vehicle. It is a good practice that all the trajectories (not only straight lines) have the same order of magnitude (between 0.00 and 5.00) for k_e. So if you have to add a scaling factor, here is were you have to do it.

The variables gvf_xx are sent to the ground station as telemetry, so you can draw the vector field, the trajectory, etc.

Finally for defining \varphi(p), \nabla\varphi(p) and H(\varphi(p)) you need to write the two files gvf_line.{c,h}



extern void gvf_line_info(float *phi, struct gvf_grad *, struct gvf_Hess *);

and in gvf_line.c you write the implicit math of your trajectory (quite explicit)



void gvf_line_info(float *phi, struct gvf_grad *grad,
        struct gvf_Hess *hess){

    struct EnuCoor_f *p = stateGetPositionEnu_f();
    float px = p->x;
    float py = p->y;
    float a = gvf_param[0];
    float b = gvf_param[1];
    float alpha = gvf_param[2];

    // Phi(x,y)
    *phi = -(px-a)*cosf(alpha) + (py-b)*sinf(alpha);

    // grad Phi
    grad->nx =  -cosf(alpha);
    grad->ny =   sinf(alpha);

    // Hessian Phi
    hess->H11 = 0;
    hess->H12 = 0;
    hess->H21 = 0;
    hess->H22 = 0;


There is a flightplan called 'demo_gvf.xml' in order to test the different settings in a simulation. Do not forget to choose default_telemetry_gvf.xml in your paparazzi center.

TO DO list

  • Better integration of the GVF with the ground station
  • Support to rotorcrafts