PX4 Firmware
PX4 Autopilot Software http://px4.io
ecl_l1_pos_controller.cpp
Go to the documentation of this file.
1 /****************************************************************************
2  *
3  * Copyright (c) 2013 Estimation and Control Library (ECL). All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  *
9  * 1. Redistributions of source code must retain the above copyright
10  * notice, this list of conditions and the following disclaimer.
11  * 2. Redistributions in binary form must reproduce the above copyright
12  * notice, this list of conditions and the following disclaimer in
13  * the documentation and/or other materials provided with the
14  * distribution.
15  * 3. Neither the name ECL nor the names of its contributors may be
16  * used to endorse or promote products derived from this software
17  * without specific prior written permission.
18  *
19  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
22  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
23  * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
24  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
25  * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
26  * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
27  * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
28  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
29  * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  * POSSIBILITY OF SUCH DAMAGE.
31  *
32  ****************************************************************************/
33 
34 /**
35  * @file ecl_l1_pos_controller.h
36  * Implementation of L1 position control.
37  * Authors and acknowledgements in header.
38  *
39  */
40 
41 #include <float.h>
42 
43 #include "ecl_l1_pos_controller.h"
44 
45 using matrix::Vector2f;
46 using matrix::wrap_pi;
47 
48 
50 {
51  float roll_new = atanf(_lateral_accel * 1.0f / CONSTANTS_ONE_G);
52  roll_new = math::constrain(roll_new, -_roll_lim_rad, _roll_lim_rad);
53 
54  if (_dt > 0.0f && _roll_slew_rate > 0.0f) {
55  // slew rate limiting active
57  }
58 
59  if (ISFINITE(roll_new)) {
60  _roll_setpoint = roll_new;
61  }
62 
63 }
64 
66 {
67  /* following [2], switching on L1 distance */
68  return math::min(wp_radius, _L1_distance);
69 }
70 
71 void
73  const Vector2f &vector_curr_position, const Vector2f &ground_speed_vector)
74 {
75  /* this follows the logic presented in [1] */
76  float eta = 0.0f;
77  float xtrack_vel = 0.0f;
78  float ltrack_vel = 0.0f;
79 
80  /* get the direction between the last (visited) and next waypoint */
81  _target_bearing = get_bearing_to_next_waypoint((double)vector_curr_position(0), (double)vector_curr_position(1), (double)vector_B(0), (double)vector_B(1));
82 
83  /* enforce a minimum ground speed of 0.1 m/s to avoid singularities */
84  float ground_speed = math::max(ground_speed_vector.length(), 0.1f);
85 
86  /* calculate the L1 length required for the desired period */
87  _L1_distance = _L1_ratio * ground_speed;
88 
89  /* calculate vector from A to B */
90  Vector2f vector_AB = get_local_planar_vector(vector_A, vector_B);
91 
92  /*
93  * check if waypoints are on top of each other. If yes,
94  * skip A and directly continue to B
95  */
96  if (vector_AB.length() < 1.0e-6f) {
97  vector_AB = get_local_planar_vector(vector_curr_position, vector_B);
98  }
99 
100  vector_AB.normalize();
101 
102  /* calculate the vector from waypoint A to the aircraft */
103  Vector2f vector_A_to_airplane = get_local_planar_vector(vector_A, vector_curr_position);
104 
105  /* calculate crosstrack error (output only) */
106  _crosstrack_error = vector_AB % vector_A_to_airplane;
107 
108  /*
109  * If the current position is in a +-135 degree angle behind waypoint A
110  * and further away from A than the L1 distance, then A becomes the L1 point.
111  * If the aircraft is already between A and B normal L1 logic is applied.
112  */
113  float distance_A_to_airplane = vector_A_to_airplane.length();
114  float alongTrackDist = vector_A_to_airplane * vector_AB;
115 
116  /* estimate airplane position WRT to B */
117  Vector2f vector_B_to_P_unit = get_local_planar_vector(vector_B, vector_curr_position).normalized();
118 
119  /* calculate angle of airplane position vector relative to line) */
120 
121  // XXX this could probably also be based solely on the dot product
122  float AB_to_BP_bearing = atan2f(vector_B_to_P_unit % vector_AB, vector_B_to_P_unit * vector_AB);
123 
124  /* extension from [2], fly directly to A */
125  if (distance_A_to_airplane > _L1_distance && alongTrackDist / math::max(distance_A_to_airplane, 1.0f) < -0.7071f) {
126 
127  /* calculate eta to fly to waypoint A */
128 
129  /* unit vector from waypoint A to current position */
130  Vector2f vector_A_to_airplane_unit = vector_A_to_airplane.normalized();
131  /* velocity across / orthogonal to line */
132  xtrack_vel = ground_speed_vector % (-vector_A_to_airplane_unit);
133  /* velocity along line */
134  ltrack_vel = ground_speed_vector * (-vector_A_to_airplane_unit);
135  eta = atan2f(xtrack_vel, ltrack_vel);
136  /* bearing from current position to L1 point */
137  _nav_bearing = atan2f(-vector_A_to_airplane_unit(1), -vector_A_to_airplane_unit(0));
138 
139  /*
140  * If the AB vector and the vector from B to airplane point in the same
141  * direction, we have missed the waypoint. At +- 90 degrees we are just passing it.
142  */
143 
144  } else if (fabsf(AB_to_BP_bearing) < math::radians(100.0f)) {
145  /*
146  * Extension, fly back to waypoint.
147  *
148  * This corner case is possible if the system was following
149  * the AB line from waypoint A to waypoint B, then is
150  * switched to manual mode (or otherwise misses the waypoint)
151  * and behind the waypoint continues to follow the AB line.
152  */
153 
154  /* calculate eta to fly to waypoint B */
155 
156  /* velocity across / orthogonal to line */
157  xtrack_vel = ground_speed_vector % (-vector_B_to_P_unit);
158  /* velocity along line */
159  ltrack_vel = ground_speed_vector * (-vector_B_to_P_unit);
160  eta = atan2f(xtrack_vel, ltrack_vel);
161  /* bearing from current position to L1 point */
162  _nav_bearing = atan2f(-vector_B_to_P_unit(1), -vector_B_to_P_unit(0));
163 
164  } else {
165 
166  /* calculate eta to fly along the line between A and B */
167 
168  /* velocity across / orthogonal to line */
169  xtrack_vel = ground_speed_vector % vector_AB;
170  /* velocity along line */
171  ltrack_vel = ground_speed_vector * vector_AB;
172  /* calculate eta2 (angle of velocity vector relative to line) */
173  float eta2 = atan2f(xtrack_vel, ltrack_vel);
174  /* calculate eta1 (angle to L1 point) */
175  float xtrackErr = vector_A_to_airplane % vector_AB;
176  float sine_eta1 = xtrackErr / math::max(_L1_distance, 0.1f);
177  /* limit output to 45 degrees */
178  sine_eta1 = math::constrain(sine_eta1, -0.7071f, 0.7071f); //sin(pi/4) = 0.7071
179  float eta1 = asinf(sine_eta1);
180  eta = eta1 + eta2;
181  /* bearing from current position to L1 point */
182  _nav_bearing = atan2f(vector_AB(1), vector_AB(0)) + eta1;
183 
184  }
185 
186  /* limit angle to +-90 degrees */
187  eta = math::constrain(eta, (-M_PI_F) / 2.0f, +M_PI_F / 2.0f);
188  _lateral_accel = _K_L1 * ground_speed * ground_speed / _L1_distance * sinf(eta);
189 
190  /* flying to waypoints, not circling them */
191  _circle_mode = false;
192 
193  /* the bearing angle, in NED frame */
194  _bearing_error = eta;
195 
197 }
198 
199 void
200 ECL_L1_Pos_Controller::navigate_loiter(const Vector2f &vector_A, const Vector2f &vector_curr_position, float radius,
201  int8_t loiter_direction, const Vector2f &ground_speed_vector)
202 {
203  /* the complete guidance logic in this section was proposed by [2] */
204 
205  /* calculate the gains for the PD loop (circle tracking) */
206  float omega = (2.0f * M_PI_F / _L1_period);
207  float K_crosstrack = omega * omega;
208  float K_velocity = 2.0f * _L1_damping * omega;
209 
210  /* update bearing to next waypoint */
211  _target_bearing = get_bearing_to_next_waypoint((double)vector_curr_position(0), (double)vector_curr_position(1), (double)vector_A(0), (double)vector_A(1));
212 
213  /* ground speed, enforce minimum of 0.1 m/s to avoid singularities */
214  float ground_speed = math::max(ground_speed_vector.length(), 0.1f);
215 
216  /* calculate the L1 length required for the desired period */
217  _L1_distance = _L1_ratio * ground_speed;
218 
219  /* calculate the vector from waypoint A to current position */
220  Vector2f vector_A_to_airplane = get_local_planar_vector(vector_A, vector_curr_position);
221 
222  Vector2f vector_A_to_airplane_unit;
223 
224  /* prevent NaN when normalizing */
225  if (vector_A_to_airplane.length() > FLT_EPSILON) {
226  /* store the normalized vector from waypoint A to current position */
227  vector_A_to_airplane_unit = vector_A_to_airplane.normalized();
228 
229  } else {
230  vector_A_to_airplane_unit = vector_A_to_airplane;
231  }
232 
233  /* calculate eta angle towards the loiter center */
234 
235  /* velocity across / orthogonal to line from waypoint to current position */
236  float xtrack_vel_center = vector_A_to_airplane_unit % ground_speed_vector;
237  /* velocity along line from waypoint to current position */
238  float ltrack_vel_center = - (ground_speed_vector * vector_A_to_airplane_unit);
239  float eta = atan2f(xtrack_vel_center, ltrack_vel_center);
240  /* limit eta to 90 degrees */
241  eta = math::constrain(eta, -M_PI_F / 2.0f, +M_PI_F / 2.0f);
242 
243  /* calculate the lateral acceleration to capture the center point */
244  float lateral_accel_sp_center = _K_L1 * ground_speed * ground_speed / _L1_distance * sinf(eta);
245 
246  /* for PD control: Calculate radial position and velocity errors */
247 
248  /* radial velocity error */
249  float xtrack_vel_circle = -ltrack_vel_center;
250  /* radial distance from the loiter circle (not center) */
251  float xtrack_err_circle = vector_A_to_airplane.length() - radius;
252 
253  /* cross track error for feedback */
254  _crosstrack_error = xtrack_err_circle;
255 
256  /* calculate PD update to circle waypoint */
257  float lateral_accel_sp_circle_pd = (xtrack_err_circle * K_crosstrack + xtrack_vel_circle * K_velocity);
258 
259  /* calculate velocity on circle / along tangent */
260  float tangent_vel = xtrack_vel_center * loiter_direction;
261 
262  /* prevent PD output from turning the wrong way */
263  if (tangent_vel < 0.0f) {
264  lateral_accel_sp_circle_pd = math::max(lateral_accel_sp_circle_pd, 0.0f);
265  }
266 
267  /* calculate centripetal acceleration setpoint */
268  float lateral_accel_sp_circle_centripetal = tangent_vel * tangent_vel / math::max((0.5f * radius),
269  (radius + xtrack_err_circle));
270 
271  /* add PD control on circle and centripetal acceleration for total circle command */
272  float lateral_accel_sp_circle = loiter_direction * (lateral_accel_sp_circle_pd + lateral_accel_sp_circle_centripetal);
273 
274  /*
275  * Switch between circle (loiter) and capture (towards waypoint center) mode when
276  * the commands switch over. Only fly towards waypoint if outside the circle.
277  */
278 
279  // XXX check switch over
280  if ((lateral_accel_sp_center < lateral_accel_sp_circle && loiter_direction > 0 && xtrack_err_circle > 0.0f) ||
281  (lateral_accel_sp_center > lateral_accel_sp_circle && loiter_direction < 0 && xtrack_err_circle > 0.0f)) {
282  _lateral_accel = lateral_accel_sp_center;
283  _circle_mode = false;
284  /* angle between requested and current velocity vector */
285  _bearing_error = eta;
286  /* bearing from current position to L1 point */
287  _nav_bearing = atan2f(-vector_A_to_airplane_unit(1), -vector_A_to_airplane_unit(0));
288 
289  } else {
290  _lateral_accel = lateral_accel_sp_circle;
291  _circle_mode = true;
292  _bearing_error = 0.0f;
293  /* bearing from current position to L1 point */
294  _nav_bearing = atan2f(-vector_A_to_airplane_unit(1), -vector_A_to_airplane_unit(0));
295  }
296 
298 }
299 
300 void ECL_L1_Pos_Controller::navigate_heading(float navigation_heading, float current_heading,
301  const Vector2f &ground_speed_vector)
302 {
303  /* the complete guidance logic in this section was proposed by [2] */
304 
305  /*
306  * As the commanded heading is the only reference
307  * (and no crosstrack correction occurs),
308  * target and navigation bearing become the same
309  */
310  _target_bearing = _nav_bearing = wrap_pi(navigation_heading);
311 
312  float eta = wrap_pi(_target_bearing - wrap_pi(current_heading));
313 
314  /* consequently the bearing error is exactly eta: */
315  _bearing_error = eta;
316 
317  /* ground speed is the length of the ground speed vector */
318  float ground_speed = ground_speed_vector.length();
319 
320  /* adjust L1 distance to keep constant frequency */
321  _L1_distance = ground_speed / _heading_omega;
322  float omega_vel = ground_speed * _heading_omega;
323 
324  /* not circling a waypoint */
325  _circle_mode = false;
326 
327  /* navigating heading means by definition no crosstrack error */
328  _crosstrack_error = 0;
329 
330  /* limit eta to 90 degrees */
331  eta = math::constrain(eta, (-M_PI_F) / 2.0f, +M_PI_F / 2.0f);
332  _lateral_accel = 2.0f * sinf(eta) * omega_vel;
333 
335 }
336 
338 {
339  /* the logic in this section is trivial, but originally proposed by [2] */
340 
341  /* reset all heading / error measures resulting in zero roll */
342  _target_bearing = current_heading;
343  _nav_bearing = current_heading;
344  _bearing_error = 0;
345  _crosstrack_error = 0;
346  _lateral_accel = 0;
347 
348  /* not circling a waypoint when flying level */
349  _circle_mode = false;
350 
352 }
353 
355 {
356  /* this is an approximation for small angles, proposed by [2] */
357  Vector2f out(math::radians((target(0) - origin(0))),
358  math::radians((target(1) - origin(1))*cosf(math::radians(origin(0)))));
359 
360  return out * static_cast<float>(CONSTANTS_RADIUS_OF_EARTH);
361 }
362 
364 {
365  _L1_period = period;
366 
367  /* calculate the ratio introduced in [2] */
368  _L1_ratio = 1.0f / M_PI_F * _L1_damping * _L1_period;
369 
370  /* calculate normalized frequency for heading tracking */
371  _heading_omega = sqrtf(2.0f) * M_PI_F / _L1_period;
372 }
373 
375 {
376  _L1_damping = damping;
377 
378  /* calculate the ratio introduced in [2] */
379  _L1_ratio = 1.0f / M_PI_F * _L1_damping * _L1_period;
380 
381  /* calculate the L1 gain (following [2]) */
382  _K_L1 = 4.0f * _L1_damping * _L1_damping;
383 }
constexpr _Tp constrain(_Tp val, _Tp min_val, _Tp max_val)
Definition: Limits.hpp:66
void set_l1_damping(float damping)
Set the L1 damping factor.
static constexpr double CONSTANTS_RADIUS_OF_EARTH
Definition: geo.h:61
float _dt
control loop time in seconds
float _roll_slew_rate
roll angle setpoint slew rate limit in rad/s
float _K_L1
L1 control gain for _L1_damping.
bool _circle_mode
flag for loiter mode
float switch_distance(float waypoint_switch_radius)
Get the switch distance.
void update_roll_setpoint()
Update roll angle setpoint.
float _target_bearing
the heading setpoint
void navigate_loiter(const matrix::Vector2f &vector_A, const matrix::Vector2f &vector_curr_position, float radius, int8_t loiter_direction, const matrix::Vector2f &ground_speed_vector)
Navigate on an orbit around a loiter waypoint.
#define FLT_EPSILON
static constexpr float CONSTANTS_ONE_G
Definition: geo.h:51
Implementation of L1 position control.
Vector normalized() const
Definition: Vector.hpp:103
float _nav_bearing
bearing to L1 reference point
float _L1_ratio
L1 ratio for navigation.
constexpr T radians(T degrees)
Definition: Limits.hpp:85
Vector< float, 6 > f(float t, const Matrix< float, 6, 1 > &, const Matrix< float, 3, 1 > &)
Definition: integration.cpp:8
Vector2< float > Vector2f
Definition: Vector2.hpp:69
float _L1_distance
L1 lead distance, defined by period and damping.
float _roll_lim_rad
maximum roll angle in radians
float get_bearing_to_next_waypoint(double lat_now, double lon_now, double lat_next, double lon_next)
Returns the bearing to the next waypoint in radians.
Definition: geo.cpp:320
void set_l1_period(float period)
Set the L1 period.
constexpr _Tp min(_Tp a, _Tp b)
Definition: Limits.hpp:54
Type wrap_pi(Type x)
Wrap value in range [-π, π)
void navigate_waypoints(const matrix::Vector2f &vector_A, const matrix::Vector2f &vector_B, const matrix::Vector2f &vector_curr_position, const matrix::Vector2f &ground_speed)
Navigate between two waypoints.
constexpr _Tp max(_Tp a, _Tp b)
Definition: Limits.hpp:60
float _heading_omega
Normalized frequency.
float _bearing_error
bearing error
float _crosstrack_error
crosstrack error in meters
float _roll_setpoint
current roll angle setpoint in radians
#define M_PI_F
Definition: ashtech.cpp:44
void navigate_level_flight(float current_heading)
Keep the wings level.
float _L1_damping
L1 damping ratio.
void navigate_heading(float navigation_heading, float current_heading, const matrix::Vector2f &ground_speed)
Navigate on a fixed bearing.
#define ISFINITE(x)
Definition: ecl.h:100
float _L1_period
L1 tracking period in seconds.
matrix::Vector2f get_local_planar_vector(const matrix::Vector2f &origin, const matrix::Vector2f &target) const
Convert a 2D vector from WGS84 to planar coordinates.
float _lateral_accel
Lateral acceleration setpoint in m/s^2.