dark0707
/
px4-firmware
forked from Archive/PX4-Autopilot
1
0
Fork 0
Code Pull Requests Activity
px4-firmware/apps/position_estimator/position_estimator_main.c

426 lines
14 KiB
C
Raw Blame History

/****************************************************************************
*
* Copyright (C) 2008-2012 PX4 Development Team. All rights reserved.
* Author: Tobias Naegeli <naegelit@student.ethz.ch>
* Thomas Gubler <thomasgubler@student.ethz.ch>
* Julian Oes <joes@student.ethz.ch>
* Lorenz Meier <lm@inf.ethz.ch>
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
*
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in
* the documentation and/or other materials provided with the
* distribution.
* 3. Neither the name PX4 nor the names of its contributors may be
* used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*
****************************************************************************/
/*
* @file Model-identification based position estimator for multirotors
*/
#include <nuttx/config.h>
#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <fcntl.h>
#include <v1.0/common/mavlink.h>
#include <float.h>
#include <nuttx/sched.h>
#include <sys/prctl.h>
#include <termios.h>
#include <errno.h>
#include <limits.h>
#include <math.h>
#include <uORB/uORB.h>
#include <uORB/topics/vehicle_status.h>
#include <uORB/topics/vehicle_attitude.h>
#include <uORB/topics/vehicle_gps_position.h>
#include <uORB/topics/vehicle_global_position.h>
#include <poll.h>
#include "codegen/position_estimator.h"
#define N_STATES 6
#define ERROR_COVARIANCE_INIT 3
#define R_EARTH 6371000.0
#define PROJECTION_INITIALIZE_COUNTER_LIMIT 5000
#define REPROJECTION_COUNTER_LIMIT 125
__EXPORT int position_estimator_main(int argc, char *argv[]);
static uint16_t position_estimator_counter_position_information;
/* values for map projection */
static double phi_1;
static double sin_phi_1;
static double cos_phi_1;
static double lambda_0;
static double scale;
/**
* Initializes the map transformation.
*
* Initializes the transformation between the geographic coordinate system and the azimuthal equidistant plane
* @param lat in degrees (47.1234567°, not 471234567°)
* @param lon in degrees (8.1234567°, not 81234567°)
*/
static void map_projection_init(double lat_0, double lon_0) //lat_0, lon_0 are expected to be in correct format: -> 47.1234567 and not 471234567
{
/* notation and formulas according to: http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html */
phi_1 = lat_0 / 180.0 * M_PI;
lambda_0 = lon_0 / 180.0 * M_PI;
sin_phi_1 = sin(phi_1);
cos_phi_1 = cos(phi_1);
/* calculate local scale by using the relation of true distance and the distance on plane */ //TODO: this is a quick solution, there are probably easier ways to determine the scale
/* 1) calculate true distance d on sphere to a point: http://www.movable-type.co.uk/scripts/latlong.html */
const double r_earth = 6371000;
double lat1 = phi_1;
double lon1 = lambda_0;
double lat2 = phi_1 + 0.5 / 180 * M_PI;
double lon2 = lambda_0 + 0.5 / 180 * M_PI;
double sin_lat_2 = sin(lat2);
double cos_lat_2 = cos(lat2);
double d = acos(sin(lat1) * sin_lat_2 + cos(lat1) * cos_lat_2 * cos(lon2 - lon1)) * r_earth;
/* 2) calculate distance rho on plane */
double k_bar = 0;
double c = acos(sin_phi_1 * sin_lat_2 + cos_phi_1 * cos_lat_2 * cos(lon2 - lambda_0));
if (0 != c)
k_bar = c / sin(c);
double x2 = k_bar * (cos_lat_2 * sin(lon2 - lambda_0)); //Projection of point 2 on plane
double y2 = k_bar * ((cos_phi_1 * sin_lat_2 - sin_phi_1 * cos_lat_2 * cos(lon2 - lambda_0)));
double rho = sqrt(pow(x2, 2) + pow(y2, 2));
scale = d / rho;
}
/**
* Transforms a point in the geographic coordinate system to the local azimuthal equidistant plane
* @param x north
* @param y east
* @param lat in degrees (47.1234567°, not 471234567°)
* @param lon in degrees (8.1234567°, not 81234567°)
*/
static void map_projection_project(double lat, double lon, float *x, float *y)
{
/* notation and formulas accoring to: http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html */
double phi = lat / 180.0 * M_PI;
double lambda = lon / 180.0 * M_PI;
double sin_phi = sin(phi);
double cos_phi = cos(phi);
double k_bar = 0;
/* using small angle approximation (formula in comment is without aproximation) */
double c = acos(sin_phi_1 * sin_phi + cos_phi_1 * cos_phi * (1 - pow((lambda - lambda_0), 2) / 2)); //double c = acos( sin_phi_1 * sin_phi + cos_phi_1 * cos_phi * cos(lambda - lambda_0) );
if (0 != c)
k_bar = c / sin(c);
/* using small angle approximation (formula in comment is without aproximation) */
*y = k_bar * (cos_phi * (lambda - lambda_0)) * scale;//*y = k_bar * (cos_phi * sin(lambda - lambda_0)) * scale;
*x = k_bar * ((cos_phi_1 * sin_phi - sin_phi_1 * cos_phi * (1 - pow((lambda - lambda_0), 2) / 2))) * scale; // *x = k_bar * ((cos_phi_1 * sin_phi - sin_phi_1 * cos_phi * cos(lambda - lambda_0))) * scale;
// printf("%phi_1=%.10f, lambda_0 =%.10f\n", phi_1, lambda_0);
}
/**
* Transforms a point in the local azimuthal equidistant plane to the geographic coordinate system
*
* @param x north
* @param y east
* @param lat in degrees (47.1234567°, not 471234567°)
* @param lon in degrees (8.1234567°, not 81234567°)
*/
static void map_projection_reproject(float x, float y, double *lat, double *lon)
{
/* notation and formulas accoring to: http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html */
double x_descaled = x / scale;
double y_descaled = y / scale;
double c = sqrt(pow(x_descaled, 2) + pow(y_descaled, 2));
double sin_c = sin(c);
double cos_c = cos(c);
double lat_sphere = 0;
if (c != 0)
lat_sphere = asin(cos_c * sin_phi_1 + (x_descaled * sin_c * cos_phi_1) / c);
else
lat_sphere = asin(cos_c * sin_phi_1);
// printf("lat_sphere = %.10f\n",lat_sphere);
double lon_sphere = 0;
if (phi_1 == M_PI / 2) {
//using small angle approximation (formula in comment is without aproximation)
lon_sphere = (lambda_0 - y_descaled / x_descaled); //lon_sphere = (lambda_0 + atan2(-y_descaled, x_descaled));
} else if (phi_1 == -M_PI / 2) {
//using small angle approximation (formula in comment is without aproximation)
lon_sphere = (lambda_0 + y_descaled / x_descaled); //lon_sphere = (lambda_0 + atan2(y_descaled, x_descaled));
} else {
lon_sphere = (lambda_0 + atan2(y_descaled * sin_c , c * cos_phi_1 * cos_c - x_descaled * sin_phi_1 * sin_c));
//using small angle approximation
// double denominator = (c * cos_phi_1 * cos_c - x_descaled * sin_phi_1 * sin_c);
// if(denominator != 0)
// {
// lon_sphere = (lambda_0 + (y_descaled * sin_c) / denominator);
// }
// else
// {
// ...
// }
}
// printf("lon_sphere = %.10f\n",lon_sphere);
*lat = lat_sphere * 180.0 / M_PI;
*lon = lon_sphere * 180.0 / M_PI;
}
/****************************************************************************
* main
****************************************************************************/
int position_estimator_main(int argc, char *argv[])
{
/* welcome user */
printf("[multirotor position_estimator] started\n");
/* initialize values */
static float u[2] = {0, 0};
static float z[3] = {0, 0, 0};
static float xapo[N_STATES] = {0, 0, 0, 0, 0, 0};
static float Papo[N_STATES * N_STATES] = {ERROR_COVARIANCE_INIT, 0, 0, 0, 0, 0,
ERROR_COVARIANCE_INIT, 0, 0, 0, 0, 0,
ERROR_COVARIANCE_INIT, 0, 0, 0, 0, 0,
ERROR_COVARIANCE_INIT, 0, 0, 0, 0, 0,
ERROR_COVARIANCE_INIT, 0, 0, 0, 0, 0,
ERROR_COVARIANCE_INIT, 0, 0, 0, 0, 0
};
static float xapo1[N_STATES];
static float Papo1[36];
static float gps_covariance[3] = {0.0f, 0.0f, 0.0f};
static uint16_t counter = 0;
position_estimator_counter_position_information = 0;
uint8_t predict_only = 1;
bool gps_valid = false;
bool new_initialization = true;
static double lat_current = 0;//[°]] --> 47.0
static double lon_current = 0; //[°]] -->8.5
//TODO: handle flight without gps but with estimator
/* subscribe to vehicle status, attitude, gps */
struct vehicle_gps_position_s gps;
struct vehicle_status_s vstatus;
struct vehicle_attitude_s att;
int vehicle_gps_sub = orb_subscribe(ORB_ID(vehicle_gps_position));
int vehicle_status_sub = orb_subscribe(ORB_ID(vehicle_status));
/* subscribe to attitude at 100 Hz */
int vehicle_attitude_sub = orb_subscribe(ORB_ID(vehicle_attitude));
/* wait until gps signal turns valid, only then can we initialize the projection */
while (!gps_valid) {
struct pollfd fds[1] = { {.fd = vehicle_gps_sub, .events = POLLIN} };
/* wait for GPS updates, BUT READ VEHICLE STATUS (!)
* this choice is critical, since the vehicle status might not
* actually change, if this app is started after GPS lock was
* aquired.
*/
if (poll(fds, 1, 5000)) {
/* Wait for the GPS update to propagate (we have some time) */
usleep(5000);
/* Read wether the vehicle status changed */
orb_copy(ORB_ID(vehicle_status), vehicle_status_sub, &vstatus);
gps_valid = vstatus.gps_valid;
}
}
/* get gps value for first initialization */
orb_copy(ORB_ID(vehicle_gps_position), vehicle_gps_sub, &gps);
lat_current = ((double)(gps.lat)) * 1e-7;
lon_current = ((double)(gps.lon)) * 1e-7;
/* publish global position messages only after first GPS message */
struct vehicle_global_position_s global_pos = {
.lat = lat_current * 1e7,
.lon = lon_current * 1e7,
.alt = gps.alt
};
orb_advert_t global_pos_pub = orb_advertise(ORB_ID(vehicle_global_position), &global_pos);
printf("[multirotor position estimator] initialized projection with: lat: %.10f, lon:%.10f\n", lat_current, lon_current);
while (1) {
/*This runs at the rate of the sensors, if we have also a new gps update this is used in the position_estimator function */
struct pollfd fds[1] = { {.fd = vehicle_attitude_sub, .events = POLLIN} };
if (poll(fds, 1, 5000) <= 0) {
/* error / timeout */
} else {
orb_copy(ORB_ID(vehicle_attitude), vehicle_attitude_sub, &att);
/* got attitude, updating pos as well */
orb_copy(ORB_ID(vehicle_gps_position), vehicle_gps_sub, &gps);
orb_copy(ORB_ID(vehicle_status), vehicle_status_sub, &vstatus);
/*copy attitude */
u[0] = att.roll;
u[1] = att.pitch;
/* initialize map projection with the last estimate (not at full rate) */
if (counter % PROJECTION_INITIALIZE_COUNTER_LIMIT == 0) {
map_projection_init(lat_current, lon_current);
new_initialization = true;
} else {
new_initialization = false;
}
/*check if new gps values are available */
gps_valid = vstatus.gps_valid;
if (gps_valid) { //we are safe to use the gps signal (it has good quality)
predict_only = 0;
/* Project gps lat lon (Geographic coordinate system) to plane*/
map_projection_project((double)(gps.lat) * 1e-7, (double)(gps.lon) * 1e-7, &(z[0]), &(z[1]));
/* copy altitude */
z[2] = (gps.alt) * 1e-3;
gps_covariance[0] = gps.eph; //TODO: needs scaling
gps_covariance[1] = gps.eph;
gps_covariance[2] = gps.epv;
} else {
/* we can not use the gps signal (it is of low quality) */
predict_only = 1;
}
// predict_only = 0; //TODO: only for testing, removeme, XXX
// z[0] = sinf(((float)counter)/180.0f*3.14159265f); //TODO: only for testing, removeme, XXX
// usleep(100000); //TODO: only for testing, removeme, XXX
/*Get new estimation (this is calculated in the plane) */
//TODO: if new_initialization == true: use 0,0,0, else use xapo
if (true == new_initialization) { //TODO,XXX: uncomment!
xapo[0] = 0; //we have a new plane initialization. the current estimate is in the center of the plane
xapo[2] = 0;
xapo[4] = 0;
position_estimator(u, z, xapo, Papo, gps_covariance, predict_only, xapo1, Papo1);
} else {
position_estimator(u, z, xapo, Papo, gps_covariance, predict_only, xapo1, Papo1);
}
/* Copy values from xapo1 to xapo */
int i;
for (i = 0; i < N_STATES; i++) {
xapo[i] = xapo1[i];
}
if ((counter % REPROJECTION_COUNTER_LIMIT == 0) || (counter % (PROJECTION_INITIALIZE_COUNTER_LIMIT - 1) == 0)) {
/* Reproject from plane to geographic coordinate system */
// map_projection_reproject(xapo1[0], xapo1[2], map_scale, phi_1, lambda_0, &lat_current, &lon_current) //TODO,XXX: uncomment!
map_projection_reproject(z[0], z[1], &lat_current, &lon_current); //do not use estimator for projection testing, removeme
// //DEBUG
// if(counter%500 == 0)
// {
// printf("phi_1: %.10f\n", phi_1);
// printf("lambda_0: %.10f\n", lambda_0);
// printf("lat_estimated: %.10f\n", lat_current);
// printf("lon_estimated: %.10f\n", lon_current);
// printf("z[0]=%.10f, z[1]=%.10f, z[2]=%f\n", z[0], z[1], z[2]);
// fflush(stdout);
//
// }
// if(!isnan(lat_current) && !isnan(lon_current))// && !isnan(xapo1[4]) && !isnan(xapo1[1]) && !isnan(xapo1[3]) && !isnan(xapo1[5]))
// {
/* send out */
global_pos.lat = lat_current;
global_pos.lon = lon_current;
global_pos.alt = xapo1[4];
global_pos.vx = xapo1[1];
global_pos.vy = xapo1[3];
global_pos.vz = xapo1[5];
/* publish current estimate */
orb_publish(ORB_ID(vehicle_global_position), global_pos_pub, &global_pos);
// }
// else
// {
// printf("[position estimator] ERROR: nan values, lat_current=%.4f, lon_current=%.4f, z[0]=%.4f z[1]=%.4f\n", lat_current, lon_current, z[0], z[1]);
// fflush(stdout);
// }
}
counter++;
}
}
return 0;
}
View Git Blame Copy Permalink