|
@@ -1,380 +0,0 @@
|
|
|
-//
|
|
|
-// Created by zx on 22-12-1.
|
|
|
-//
|
|
|
-
|
|
|
-#include "loaded_mpc.h"
|
|
|
-#include <cppad/cppad.hpp>
|
|
|
-#include <cppad/ipopt/solve.hpp>
|
|
|
-
|
|
|
-size_t N = 20; //优化考虑后面多少步
|
|
|
-
|
|
|
-size_t nx = 0;
|
|
|
-size_t ny = nx + N;
|
|
|
-size_t nth = ny + N;
|
|
|
-size_t nv = nth + N;
|
|
|
-size_t ndlt = nv + N;
|
|
|
-
|
|
|
-class FG_eval_half_agv {
|
|
|
- public:
|
|
|
- // Fitted polynomial coefficients
|
|
|
- Eigen::VectorXd m_coeffs; //曲线方程
|
|
|
- Eigen::VectorXd m_statu; //当前状态
|
|
|
- Eigen::VectorXd m_condition; //搜索条件参数
|
|
|
- FG_eval_half_agv(Eigen::VectorXd coeffs,Eigen::VectorXd statu,Eigen::VectorXd condition) {
|
|
|
- m_coeffs=coeffs;
|
|
|
- m_statu=statu;
|
|
|
- m_condition=condition;
|
|
|
- }
|
|
|
-
|
|
|
- typedef CPPAD_TESTVECTOR(CppAD::AD<double>) ADvector;
|
|
|
- void operator()(ADvector& fg, const ADvector& vars) {
|
|
|
- // Implementing MPC below
|
|
|
- // `fg` a vector of the cost constraints, `vars` is a vector of variable values (state & actuators)
|
|
|
- // The cost is stored is the first element of `fg`.
|
|
|
- // Any additions to the cost should be added to `fg[0]`.
|
|
|
- fg[0] = 0;
|
|
|
- double dt=m_condition[0];
|
|
|
- double ref_v=m_condition[1];
|
|
|
- double v=m_statu[0];
|
|
|
- double delta=m_statu[1];
|
|
|
- double wheelbase=m_statu[2];
|
|
|
-
|
|
|
- // Reference State Cost
|
|
|
- // Below defines the cost related the reference state and
|
|
|
- // any anything you think may be beneficial.
|
|
|
-
|
|
|
- // Weights for how "important" each cost is - can be tuned
|
|
|
- const int y_cost_weight = 12000;
|
|
|
- const int th_cost_weight = 1000;
|
|
|
- const int v_cost_weight = 20;
|
|
|
- const int vth_cost_weight = 100;
|
|
|
- const int a_cost_weight = 1;
|
|
|
- const int ath_cost_weight=10;
|
|
|
-
|
|
|
- // Cost for CTE, psi error and velocity
|
|
|
- for (int t = 0; t < N; t++) {
|
|
|
- CppAD::AD<double> xt=vars[nx+t];
|
|
|
- CppAD::AD<double> fx = m_coeffs[0] + m_coeffs[1] * xt + m_coeffs[2] * pow(xt, 2) + m_coeffs[3] * pow(xt, 3);
|
|
|
- fg[0] += y_cost_weight * CppAD::pow(vars[ny+t]-fx, 2);
|
|
|
-
|
|
|
- //朝向loss
|
|
|
- CppAD::AD<double> fth = CppAD::atan(m_coeffs[1] + 2*m_coeffs[2]*xt + 3*m_coeffs[3]*pow(xt,2));
|
|
|
- fg[0] += th_cost_weight * CppAD::pow(vars[nth + t]-fth, 2);
|
|
|
-
|
|
|
- //目标速度loss
|
|
|
- fg[0]+=v_cost_weight*CppAD::pow(vars[nv+t]-ref_v,2);
|
|
|
-
|
|
|
- }
|
|
|
-
|
|
|
- // Costs for steering (delta) and acceleration (a)
|
|
|
-
|
|
|
- for (int t = 0; t < N-1; t++) {
|
|
|
- //速度,加速度,前轮角 weight loss
|
|
|
- fg[0] += vth_cost_weight * CppAD::pow(vars[ndlt+t], 2);
|
|
|
- fg[0] += a_cost_weight * CppAD::pow(vars[nv + t+1]-vars[nv+t], 2);
|
|
|
- fg[0] += ath_cost_weight * CppAD::pow(vars[ndlt+t+1]-vars[ndlt+t], 2);
|
|
|
- }
|
|
|
- fg[0]*=0.01;
|
|
|
-
|
|
|
- /////////////////////
|
|
|
- fg[1 + nx] = vars[nx]-vars[nv]*dt;
|
|
|
- fg[1 + ny] = vars[ny];
|
|
|
- CppAD::AD<double> w0=vars[nv]/wheelbase*CppAD::tan(vars[ndlt]);
|
|
|
- fg[1 + nth] = vars[nth]-w0*dt;
|
|
|
-
|
|
|
- //位姿约束
|
|
|
- for (int t = 1; t < N; t++) {
|
|
|
- // State at time t + 1
|
|
|
- CppAD::AD<double> x1 = vars[nx + t];
|
|
|
- CppAD::AD<double> y1 = vars[ny + t];
|
|
|
- CppAD::AD<double> th1 = vars[nth + t];
|
|
|
-
|
|
|
- // State at time t
|
|
|
- CppAD::AD<double> x0 = vars[nx + t -1];
|
|
|
- CppAD::AD<double> y0 = vars[ny + t -1];
|
|
|
- CppAD::AD<double> th0 = vars[nth + t-1];
|
|
|
- CppAD::AD<double> v0 = vars[nv + t-1];
|
|
|
-
|
|
|
- CppAD::AD<double> w_0=vars[nv+t-1]/wheelbase*CppAD::tan(vars[ndlt+t-1]);
|
|
|
-
|
|
|
- // Setting up the rest of the model constraints
|
|
|
- fg[1 + nx + t] = x1 - (x0 + v0 * CppAD::cos(th0) * dt);
|
|
|
- fg[1 + ny + t] = y1 - (y0 + v0 * CppAD::sin(th0) * dt);
|
|
|
- fg[1 + nth + t] = th1 - (th0 + w_0 * dt);
|
|
|
-
|
|
|
- }
|
|
|
- //加速度和dlt约束
|
|
|
- fg[1 + nv]=(vars[nv]-v)/dt;
|
|
|
- fg[1+ndlt]=(vars[ndlt]-delta)/dt;
|
|
|
- for(int t=1;t<N;++t)
|
|
|
- {
|
|
|
- fg[1+nv+t]=(vars[nv+t]-vars[nv+t-1])/dt;
|
|
|
- fg[1+ndlt+t]=(vars[ndlt+t]-vars[ndlt+t-1])/dt;
|
|
|
- }
|
|
|
-
|
|
|
- }
|
|
|
-};
|
|
|
-
|
|
|
-LoadedMPC::LoadedMPC(double wheelbase){
|
|
|
- wheelbase_=wheelbase;
|
|
|
-}
|
|
|
-LoadedMPC::~LoadedMPC(){}
|
|
|
-
|
|
|
-bool LoadedMPC::solve(const Trajectory& trajectory, const Pose2d& target,Eigen::VectorXd statu,
|
|
|
- std::vector<double>& out,Trajectory& select_traj,Trajectory& optimize_trajectory)
|
|
|
-{
|
|
|
- auto start=std::chrono::steady_clock::now();
|
|
|
- // State vector holds all current values neede for vars below
|
|
|
- Pose2d pose_agv = Pose2d(statu[0], statu[1], statu[2]);
|
|
|
- double line_velocity = statu[3];
|
|
|
- double angular_velocity = statu[4];
|
|
|
- std::vector<Pose2d> filte_poses;
|
|
|
- if (filte_Path(pose_agv, target, trajectory, filte_poses, 10) == false)
|
|
|
- {
|
|
|
- printf( "filte path failed ...\n");
|
|
|
- return false;
|
|
|
- }
|
|
|
- select_traj= Trajectory(filte_poses);
|
|
|
-
|
|
|
- //将选中点移动到小车坐标系
|
|
|
- std::vector<Pose2d> transform_poses;
|
|
|
- for (int i = 0; i < filte_poses.size(); i++)
|
|
|
- {
|
|
|
- double x = filte_poses[i].x() - pose_agv.x();
|
|
|
- double y = filte_poses[i].y() - pose_agv.y();
|
|
|
- transform_poses.push_back(Pose2d(x,y,0).rotate(-pose_agv.theta()));
|
|
|
- }
|
|
|
-
|
|
|
- Eigen::VectorXd coef = fit_path(transform_poses);
|
|
|
-
|
|
|
- //优化
|
|
|
- bool ok = true;
|
|
|
- typedef CPPAD_TESTVECTOR(double) Dvector;
|
|
|
-
|
|
|
- //根据当前点和目标点距离,计算目标速度
|
|
|
-
|
|
|
- double ref_velocity = Pose2d::distance(pose_agv, target) / 2.5;
|
|
|
- //目标点与起点的连线 朝向与启动朝向 > M_PI/2.0
|
|
|
-
|
|
|
- float gradient=fabs(target.x()-pose_agv.x())>1e-6?(target.y()-pose_agv.y())/(target.x()-pose_agv.x()):500.0;
|
|
|
- double direction=gradient2theta(gradient,target.x()>=pose_agv.x());
|
|
|
-
|
|
|
- // pose_agv.theta()*k,direction*k,k*fabs(pose_agv.theta()-direction));
|
|
|
- if (fabs(direction-pose_agv.theta()) >M_PI/2.0)
|
|
|
- ref_velocity = -ref_velocity;
|
|
|
-
|
|
|
-
|
|
|
- double dt = 0.1;
|
|
|
- double min_velocity = 0.05;
|
|
|
- double max_velocity = 1.9;
|
|
|
- double max_dlt=5*M_PI/180.0;
|
|
|
- double max_acc_line_velocity=1.9/(N*dt);
|
|
|
- double max_acc_dlt=5*M_PI/180.0;
|
|
|
-
|
|
|
- size_t n_vars = N * 5;
|
|
|
- size_t n_constraints = N * 5;
|
|
|
-
|
|
|
- Dvector vars(n_vars);
|
|
|
- for (int i = 0; i < n_vars; i++)
|
|
|
- {
|
|
|
- vars[i] = 0.0;
|
|
|
- }
|
|
|
-
|
|
|
- Dvector vars_lowerbound(n_vars);
|
|
|
- Dvector vars_upperbound(n_vars);
|
|
|
-
|
|
|
- for (int i = 0; i < n_vars; i++)
|
|
|
- {
|
|
|
- vars_lowerbound[i] = -1.0e19;
|
|
|
- vars_upperbound[i] = 1.0e19;
|
|
|
- }
|
|
|
-
|
|
|
- //// limit v
|
|
|
- for (int i = nv; i < nv + N; i++)
|
|
|
- {
|
|
|
- vars_lowerbound[i] = min_velocity;
|
|
|
- vars_upperbound[i] = max_velocity;
|
|
|
- }
|
|
|
-
|
|
|
- ////limint dlt
|
|
|
- for (int i = ndlt; i < ndlt + N; i++)
|
|
|
- {
|
|
|
- vars_lowerbound[i] = -max_dlt;
|
|
|
- vars_upperbound[i] = max_dlt;
|
|
|
- }
|
|
|
-
|
|
|
- // Lower and upper limits for the constraints
|
|
|
- Dvector constraints_lowerbound(n_constraints);
|
|
|
- Dvector constraints_upperbound(n_constraints);
|
|
|
- for (int i = 0; i < n_constraints; i++)
|
|
|
- {
|
|
|
- constraints_lowerbound[i] = 0;
|
|
|
- constraints_upperbound[i] = 0;
|
|
|
- }
|
|
|
- //// acc v
|
|
|
- for (int i = nv; i < nv + N; i++)
|
|
|
- {
|
|
|
- constraints_lowerbound[i] = -max_acc_line_velocity;
|
|
|
- constraints_upperbound[i] = max_acc_line_velocity;
|
|
|
- }
|
|
|
- //// acc ndlt
|
|
|
- for (int i = ndlt; i < ndlt + N; i++)
|
|
|
- {
|
|
|
- constraints_lowerbound[i] = -max_acc_dlt;
|
|
|
- constraints_upperbound[i] = max_acc_dlt;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- Eigen::VectorXd statu_velocity(3);
|
|
|
- if(line_velocity>max_velocity)
|
|
|
- {
|
|
|
- line_velocity=max_velocity;
|
|
|
- printf(" input +v limited\n");
|
|
|
- }
|
|
|
- if(line_velocity<min_velocity)
|
|
|
- {
|
|
|
- line_velocity=min_velocity;
|
|
|
- printf(" input -v limited\n");
|
|
|
- }
|
|
|
-
|
|
|
- double omga=0;
|
|
|
- if(abs(line_velocity)>0.000001)
|
|
|
- {
|
|
|
- omga=angular_velocity*wheelbase_/tan(line_velocity);
|
|
|
- }
|
|
|
- if(omga>max_dlt)
|
|
|
- {
|
|
|
- omga = max_dlt;
|
|
|
- printf(" input +dlt limited\n");
|
|
|
- }
|
|
|
- if(omga<-max_dlt)
|
|
|
- {
|
|
|
- omga = -max_dlt;
|
|
|
- printf(" input -dlt limited\n");
|
|
|
- }
|
|
|
-
|
|
|
- statu_velocity << line_velocity, omga,wheelbase_;
|
|
|
-
|
|
|
- Eigen::VectorXd condition(2);
|
|
|
- condition << dt, ref_velocity;
|
|
|
- FG_eval_half_agv fg_eval(coef, statu_velocity, condition);
|
|
|
-
|
|
|
- // options for IPOPT solver
|
|
|
- std::string options;
|
|
|
- // Uncomment this if you'd like more print information
|
|
|
- options += "Integer print_level 0\n";
|
|
|
- options += "Sparse true forward\n";
|
|
|
- options += "Sparse true reverse\n";
|
|
|
- options += "Numeric max_cpu_time 0.5\n";
|
|
|
- // place to return solution
|
|
|
- CppAD::ipopt::solve_result<Dvector> solution;
|
|
|
-
|
|
|
- // solve the problem
|
|
|
- CppAD::ipopt::solve<Dvector, FG_eval_half_agv>(
|
|
|
- options, vars, vars_lowerbound, vars_upperbound, constraints_lowerbound,
|
|
|
- constraints_upperbound, fg_eval, solution);
|
|
|
-
|
|
|
- auto now=std::chrono::steady_clock::now();
|
|
|
- auto duration = std::chrono::duration_cast<std::chrono::microseconds>(now - start);
|
|
|
- double time = double(duration.count()) * std::chrono::microseconds::period::num / std::chrono::microseconds::period::den;
|
|
|
-
|
|
|
- ok &= solution.status == CppAD::ipopt::solve_result<Dvector>::success;
|
|
|
- if (ok == false)
|
|
|
- {
|
|
|
- printf(" mpc failed statu : %d input: %.4f %.4f \n",solution.status,line_velocity,
|
|
|
- angular_velocity);
|
|
|
- return false;
|
|
|
- }
|
|
|
-
|
|
|
- // Cost
|
|
|
- auto cost = solution.obj_value;
|
|
|
-
|
|
|
- out.clear();
|
|
|
-
|
|
|
- double solve_velocity=solution.x[nv];
|
|
|
- /*if(fabs(solve_velocity)<min_velocity)
|
|
|
- {
|
|
|
- solve_velocity= solve_velocity>0?min_velocity:-min_velocity;
|
|
|
- }*/
|
|
|
-
|
|
|
- out.push_back(solve_velocity);
|
|
|
- double omg=solve_velocity/wheelbase_*tan(solution.x[ndlt]);
|
|
|
- out.push_back(omg);
|
|
|
- printf(" input : %.4f %.4f output : %.4f %.4f ref : %.3f time:%.3f\n",line_velocity,
|
|
|
- angular_velocity,solve_velocity,omg,ref_velocity,time);
|
|
|
-
|
|
|
- optimize_trajectory.clear();
|
|
|
- for (int i = 0; i < N; ++i)
|
|
|
- {
|
|
|
- Pose2d pose(solution.x[nx + i], solution.x[ny + i], solution.x[nth + i]);
|
|
|
- optimize_trajectory.push_point(pose_agv+pose.rotate(pose_agv.theta()));
|
|
|
- }
|
|
|
-
|
|
|
- return true;
|
|
|
-}
|
|
|
-
|
|
|
-bool LoadedMPC::filte_Path(const Pose2d& point,const Pose2d& target,const Trajectory& trajectory,
|
|
|
- std::vector<Pose2d>& poses,int point_num)
|
|
|
-{
|
|
|
- double gradient=0;
|
|
|
- if(fabs((target-point).x())==0)
|
|
|
- gradient=200.0;
|
|
|
- else
|
|
|
- gradient= (target-point).y()/(target-point).x();
|
|
|
-
|
|
|
- double theta=gradient2theta(gradient,target.x()-point.x()>=0);
|
|
|
-
|
|
|
- if(trajectory.size() < point_num)
|
|
|
- return false;
|
|
|
-
|
|
|
- poses.clear();
|
|
|
-
|
|
|
- //double theta=point.theta();
|
|
|
-
|
|
|
- //if(fabs(theta)>M_PI/2.0) theta+=M_PI;
|
|
|
- for (int i = 0; i < trajectory.size(); i++)
|
|
|
- {
|
|
|
-
|
|
|
- // 平移加反向旋转到小车坐标系
|
|
|
- Pose2d offset=trajectory[i]-point;
|
|
|
- double x = offset.x();
|
|
|
- double y = offset.y();
|
|
|
- double trans_x = x * cos(-theta) - y * sin(-theta);
|
|
|
- double trans_y = x * sin(-theta) + y * cos(-theta);
|
|
|
- if (trans_x>=0 && poses.size() < point_num)
|
|
|
- {
|
|
|
-
|
|
|
- // 旋转到原坐标系
|
|
|
- float nx=trans_x * cos(theta) - trans_y * sin(theta);
|
|
|
- float ny=trans_x*sin(theta)+trans_y*cos(theta);
|
|
|
- Pose2d pose(nx+point.x(),ny+point.y(),trajectory[i].theta());
|
|
|
-
|
|
|
- poses.push_back(pose);
|
|
|
-
|
|
|
- }
|
|
|
- }
|
|
|
- return true;
|
|
|
-}
|
|
|
-
|
|
|
-Eigen::VectorXd LoadedMPC::fit_path(const std::vector<Pose2d>& trajectory)
|
|
|
-{
|
|
|
- int order=3;
|
|
|
- assert(order >= 1 && order <= trajectory.size() - 1);
|
|
|
- Eigen::MatrixXd A(trajectory.size(), order + 1);
|
|
|
- Eigen::VectorXd yvals(trajectory.size());
|
|
|
- for (int i = 0; i < trajectory.size(); i++) {
|
|
|
- A(i, 0) = 1.0;
|
|
|
- yvals[i]=trajectory[i].y();
|
|
|
- }
|
|
|
-
|
|
|
- for (int j = 0; j < trajectory.size(); j++) {
|
|
|
- for (int i = 0; i < order; i++) {
|
|
|
- A(j, i + 1) = A(j, i) * trajectory[j].x();
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- auto Q = A.householderQr();
|
|
|
- auto result = Q.solve(yvals);
|
|
|
- return result;
|
|
|
-}
|